diff --git a/python/metatomic_ase/setup.py b/python/metatomic_ase/setup.py
index 3d57b23f..a45fe0ab 100644
--- a/python/metatomic_ase/setup.py
+++ b/python/metatomic_ase/setup.py
@@ -122,7 +122,7 @@ def create_version_number(version):
         # No dependency on ASE itself until this package is no longer a direct
         # dependency of metatomic-torch
         # "ase >=3.22.0",
-        "vesin >=0.5.2,<0.6",
+        "vesin >=0.5.5,<0.6",
     ]
 
     # when packaging a sdist for release, we should never use local dependencies
diff --git a/python/metatomic_ase/src/metatomic_ase/_calculator.py b/python/metatomic_ase/src/metatomic_ase/_calculator.py
index 2247a054..b6f52b5c 100644
--- a/python/metatomic_ase/src/metatomic_ase/_calculator.py
+++ b/python/metatomic_ase/src/metatomic_ase/_calculator.py
@@ -21,7 +21,7 @@
     pick_output,
 )
 
-from ._neighbors import _compute_requested_neighbors
+from ._neighbors import AllNeighborsCalculator
 
 
 import ase  # isort: skip
@@ -317,6 +317,11 @@ def __init__(
                     "be positive"
                 )
 
+        self._nl_calculators = AllNeighborsCalculator(
+            requested_options=self._model.requested_neighbor_lists(),
+            check_consistency=check_consistency,
+        )
+
         # We do our own check to verify if a property is implemented in `calculate()`,
         # so we pretend to be able to compute all properties ASE knows about.
         self.implemented_properties = ALL_ASE_PROPERTIES
@@ -398,11 +403,7 @@ def run_model(
             systems.append(system)
 
         # Compute the neighbors lists requested by the model
-        input_systems = _compute_requested_neighbors(
-            systems=systems,
-            requested_options=self._model.requested_neighbor_lists(),
-            check_consistency=self.parameters["check_consistency"],
-        )
+        input_systems = self._nl_calculators.compute(systems=systems)
 
         available_outputs = self._model.capabilities().outputs
         for key in outputs:
@@ -538,11 +539,7 @@ def calculate(
         with record_function("MetatomicCalculator::compute_neighbors"):
             # convert from ase.Atoms to metatomic.torch.System
             system = System(types, positions, cell, pbc)
-            input_system = _compute_requested_neighbors(
-                systems=[system],
-                requested_options=self._model.requested_neighbor_lists(),
-                check_consistency=self.parameters["check_consistency"],
-            )[0]
+            input_system = self._nl_calculators.compute(systems=[system])[0]
 
         with record_function("MetatomicCalculator::get_model_inputs"):
             for name, option in self._model.requested_inputs().items():
@@ -721,11 +718,7 @@ def compute_energy(
             systems.append(system)
 
         # Compute the neighbors lists requested by the model
-        input_systems = _compute_requested_neighbors(
-            systems=systems,
-            requested_options=self._model.requested_neighbor_lists(),
-            check_consistency=self.parameters["check_consistency"],
-        )
+        input_systems = self._nl_calculators.compute(systems=systems)
 
         predictions = self._model(
             systems=input_systems,
diff --git a/python/metatomic_ase/src/metatomic_ase/_neighbors.py b/python/metatomic_ase/src/metatomic_ase/_neighbors.py
index ddbf2d95..eac31412 100644
--- a/python/metatomic_ase/src/metatomic_ase/_neighbors.py
+++ b/python/metatomic_ase/src/metatomic_ase/_neighbors.py
@@ -31,57 +31,68 @@
     HAS_NVALCHEMIOPS = False
 
 
-def _compute_requested_neighbors(
-    systems: List[System],
-    requested_options: List[NeighborListOptions],
-    check_consistency=False,
-) -> List[System]:
-    """
-    Compute all neighbor lists requested by ``model`` and store them inside the systems.
-    """
-    can_use_nvalchemi = HAS_NVALCHEMIOPS and all(
-        system.device.type == "cuda" for system in systems
-    )
-
-    if can_use_nvalchemi:
-        full_nl_options = []
-        half_nl_options = []
-        for options in requested_options:
-            if options.full_list:
-                full_nl_options.append(options)
-            else:
-                half_nl_options.append(options)
-
-        # Do the full neighbor lists with nvalchemi, and the rest with vesin
-        systems = _compute_requested_neighbors_nvalchemi(
-            systems=systems,
-            requested_options=full_nl_options,
-        )
-        systems = _compute_requested_neighbors_vesin(
-            systems=systems,
-            requested_options=half_nl_options,
-            check_consistency=check_consistency,
+class AllNeighborsCalculator:
+    def __init__(
+        self,
+        requested_options: List[NeighborListOptions],
+        check_consistency=False,
+    ):
+        self.check_consistency = check_consistency
+        self._full_nl_options = [
+            options for options in requested_options if options.full_list
+        ]
+        self._full_vesin_calculators = [
+            vesin.metatomic.NeighborList(
+                options=options,
+                length_unit="angstrom",
+                check_consistency=check_consistency,
+            )
+            for options in requested_options
+            if options.full_list
+        ]
+        self._half_vesin_calculators = [
+            vesin.metatomic.NeighborList(
+                options=options,
+                length_unit="angstrom",
+                check_consistency=check_consistency,
+            )
+            for options in requested_options
+            if not options.full_list
+        ]
+
+    def compute(self, systems: List[System]) -> List[System]:
+        assert isinstance(systems, list)
+        assert isinstance(systems[0], torch.ScriptObject)
+
+        can_use_nvalchemi = HAS_NVALCHEMIOPS and all(
+            system.device.type == "cuda" for system in systems
         )
-    else:
+
+        if can_use_nvalchemi:
+            # Do the full neighbor lists with nvalchemi
+            systems = _compute_requested_neighbors_nvalchemi(
+                systems=systems,
+                requested_options=self._full_nl_options,
+            )
+        else:
+            systems = _compute_requested_neighbors_vesin(
+                systems=systems,
+                calculators=self._full_vesin_calculators,
+            )
+
+        # always compute the half neighbor lists with vesin
         systems = _compute_requested_neighbors_vesin(
             systems=systems,
-            requested_options=requested_options,
-            check_consistency=check_consistency,
+            calculators=self._half_vesin_calculators,
         )
 
-    return systems
+        return systems
 
 
 def _compute_requested_neighbors_vesin(
     systems: List[System],
-    requested_options: List[NeighborListOptions],
-    check_consistency=False,
+    calculators: List[vesin.metatomic.NeighborList],
 ) -> List[System]:
-    """
-    Compute all neighbor lists requested by ``model`` and store them inside the systems,
-    using vesin.
-    """
-
     system_devices = []
     moved_systems = []
     for system in systems:
@@ -91,12 +102,13 @@ def _compute_requested_neighbors_vesin(
         else:
             moved_systems.append(system)
 
-    vesin.metatomic.compute_requested_neighbors_from_options(
-        systems=moved_systems,
-        system_length_unit="angstrom",
-        options=requested_options,
-        check_consistency=check_consistency,
-    )
+    for calculator in calculators:
+        calculator.add_neighbor_list(
+            systems=moved_systems,
+            # if we have more than one system, we can no keep the data as a reference
+            # to memory allocated in the calculator and we need to make a copy
+            copy=len(systems) > 1,
+        )
 
     systems = []
     for system, device in zip(moved_systems, system_devices, strict=True):
@@ -142,6 +154,7 @@ def _compute_requested_neighbors_nvalchemi(systems, requested_options):
                         "cell_shift_c",
                     ],
                     values=torch.hstack([P, S]),
+                    assume_unique=True,
                 ),
                 components=[
                     Labels("xyz", torch.tensor([[0], [1], [2]], device=system.device))
diff --git a/python/metatomic_torchsim/CHANGELOG.md b/python/metatomic_torchsim/CHANGELOG.md
index bf24c999..561dbc5a 100644
--- a/python/metatomic_torchsim/CHANGELOG.md
+++ b/python/metatomic_torchsim/CHANGELOG.md
@@ -16,6 +16,13 @@ follows [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased](https://github.com/metatensor/metatomic/)
 
+## [Version 0.1.2](https://github.com/metatensor/metatomic/releases/tag/metatomic-torchsim-v0.1.2) - 2026-04-22
+
+### Changed
+
+- Removed the upper-version pin on `torch-sim-atomistic` to make updating the
+  code in there that re-exports this package easier.
+
 ## [Version 0.1.1](https://github.com/metatensor/metatomic/releases/tag/metatomic-torchsim-v0.1.1) - 2026-04-01
 
 ### Fixed
diff --git a/python/metatomic_torchsim/metatomic_torchsim/_model.py b/python/metatomic_torchsim/metatomic_torchsim/_model.py
index b42fc6d9..2318f90b 100644
--- a/python/metatomic_torchsim/metatomic_torchsim/_model.py
+++ b/python/metatomic_torchsim/metatomic_torchsim/_model.py
@@ -28,7 +28,7 @@
     pick_output,
 )
 
-from ._neighbors import _compute_requested_neighbors
+from ._neighbors import AllNeighborsCalculator
 
 
 try:
@@ -249,7 +249,6 @@ def __init__(
                     "be positive"
                 )
 
-        self._requested_neighbor_lists = self._model.requested_neighbor_lists()
         self._requested_inputs = self._model.requested_inputs()
         if len(self._requested_inputs) != 0:
             raise ValueError(
@@ -283,6 +282,11 @@ def __init__(
             outputs=run_outputs,
         )
 
+        self._nl_calculators = AllNeighborsCalculator(
+            requested_options=self._model.requested_neighbor_lists(),
+            check_consistency=check_consistency,
+        )
+
         self.additional_outputs: Dict[str, TensorMap] = {}
         """
         Additional outputs computed by :py:meth:`forward` are stored here.
@@ -355,11 +359,7 @@ def forward(self, state: "ts.SimState") -> Dict[str, torch.Tensor]:
             )
 
         # Compute neighbor lists
-        systems = _compute_requested_neighbors(
-            systems=systems,
-            requested_options=self._requested_neighbor_lists,
-            check_consistency=self._check_consistency,
-        )
+        systems = self._nl_calculators.compute(systems=systems)
 
         # Run the model (evaluation options precomputed in __init__)
         model_outputs = self._model(
diff --git a/python/metatomic_torchsim/metatomic_torchsim/_neighbors.py b/python/metatomic_torchsim/metatomic_torchsim/_neighbors.py
index e343fb14..eac31412 100644
--- a/python/metatomic_torchsim/metatomic_torchsim/_neighbors.py
+++ b/python/metatomic_torchsim/metatomic_torchsim/_neighbors.py
@@ -24,7 +24,6 @@
             "If you encounter errors, please update to this version range.",
             stacklevel=1,
         )
-
     from nvalchemiops.torch.neighbors import neighbor_list as nvalchemi_neighbor_list
 
     HAS_NVALCHEMIOPS = True
@@ -32,53 +31,68 @@
     HAS_NVALCHEMIOPS = False
 
 
-def _compute_requested_neighbors(
-    systems: List[System],
-    requested_options: List[NeighborListOptions],
-    check_consistency: bool = False,
-) -> List[System]:
-    """Compute all neighbor lists requested by the model and store them in the systems.
+class AllNeighborsCalculator:
+    def __init__(
+        self,
+        requested_options: List[NeighborListOptions],
+        check_consistency=False,
+    ):
+        self.check_consistency = check_consistency
+        self._full_nl_options = [
+            options for options in requested_options if options.full_list
+        ]
+        self._full_vesin_calculators = [
+            vesin.metatomic.NeighborList(
+                options=options,
+                length_unit="angstrom",
+                check_consistency=check_consistency,
+            )
+            for options in requested_options
+            if options.full_list
+        ]
+        self._half_vesin_calculators = [
+            vesin.metatomic.NeighborList(
+                options=options,
+                length_unit="angstrom",
+                check_consistency=check_consistency,
+            )
+            for options in requested_options
+            if not options.full_list
+        ]
 
-    Uses nvalchemiops for full neighbor lists on CUDA when available, vesin otherwise.
-    """
-    can_use_nvalchemi = HAS_NVALCHEMIOPS and all(
-        system.device.type == "cuda" for system in systems
-    )
-
-    if can_use_nvalchemi:
-        full_nl_options = []
-        half_nl_options = []
-        for options in requested_options:
-            if options.full_list:
-                full_nl_options.append(options)
-            else:
-                half_nl_options.append(options)
-
-        systems = _compute_requested_neighbors_nvalchemi(
-            systems=systems,
-            requested_options=full_nl_options,
-        )
-        systems = _compute_requested_neighbors_vesin(
-            systems=systems,
-            requested_options=half_nl_options,
-            check_consistency=check_consistency,
+    def compute(self, systems: List[System]) -> List[System]:
+        assert isinstance(systems, list)
+        assert isinstance(systems[0], torch.ScriptObject)
+
+        can_use_nvalchemi = HAS_NVALCHEMIOPS and all(
+            system.device.type == "cuda" for system in systems
         )
-    else:
+
+        if can_use_nvalchemi:
+            # Do the full neighbor lists with nvalchemi
+            systems = _compute_requested_neighbors_nvalchemi(
+                systems=systems,
+                requested_options=self._full_nl_options,
+            )
+        else:
+            systems = _compute_requested_neighbors_vesin(
+                systems=systems,
+                calculators=self._full_vesin_calculators,
+            )
+
+        # always compute the half neighbor lists with vesin
         systems = _compute_requested_neighbors_vesin(
             systems=systems,
-            requested_options=requested_options,
-            check_consistency=check_consistency,
+            calculators=self._half_vesin_calculators,
         )
 
-    return systems
+        return systems
 
 
 def _compute_requested_neighbors_vesin(
     systems: List[System],
-    requested_options: List[NeighborListOptions],
-    check_consistency: bool = False,
+    calculators: List[vesin.metatomic.NeighborList],
 ) -> List[System]:
-    """Compute neighbor lists using vesin."""
     system_devices = []
     moved_systems = []
     for system in systems:
@@ -88,12 +102,13 @@ def _compute_requested_neighbors_vesin(
         else:
             moved_systems.append(system)
 
-    vesin.metatomic.compute_requested_neighbors_from_options(
-        systems=moved_systems,
-        system_length_unit="angstrom",
-        options=requested_options,
-        check_consistency=check_consistency,
-    )
+    for calculator in calculators:
+        calculator.add_neighbor_list(
+            systems=moved_systems,
+            # if we have more than one system, we can no keep the data as a reference
+            # to memory allocated in the calculator and we need to make a copy
+            copy=len(systems) > 1,
+        )
 
     systems = []
     for system, device in zip(moved_systems, system_devices, strict=True):
@@ -102,11 +117,13 @@ def _compute_requested_neighbors_vesin(
     return systems
 
 
-def _compute_requested_neighbors_nvalchemi(
-    systems: List[System],
-    requested_options: List[NeighborListOptions],
-) -> List[System]:
-    """Compute full neighbor lists on CUDA using nvalchemiops."""
+def _compute_requested_neighbors_nvalchemi(systems, requested_options):
+    """
+    Compute all neighbor lists requested by ``model`` and store them inside the systems,
+    using nvalchemiops. This function should only be called if all systems are on CUDA
+    and all neighbor list options require a full neighbor list.
+    """
+
     for options in requested_options:
         assert options.full_list
         for system in systems:
@@ -137,16 +154,13 @@ def _compute_requested_neighbors_nvalchemi(
                         "cell_shift_c",
                     ],
                     values=torch.hstack([P, S]),
+                    assume_unique=True,
                 ),
                 components=[
-                    Labels(
-                        "xyz",
-                        torch.tensor([[0], [1], [2]], device=system.device),
-                    )
+                    Labels("xyz", torch.tensor([[0], [1], [2]], device=system.device))
                 ],
                 properties=Labels(
-                    "distance",
-                    torch.tensor([[0]], device=system.device),
+                    "distance", torch.tensor([[0]], device=system.device)
                 ),
             )
             system.add_neighbor_list(options, neighbors)
diff --git a/python/metatomic_torchsim/setup.py b/python/metatomic_torchsim/setup.py
index fd63256e..4bd6ab05 100644
--- a/python/metatomic_torchsim/setup.py
+++ b/python/metatomic_torchsim/setup.py
@@ -10,7 +10,7 @@
 ROOT = os.path.realpath(os.path.dirname(__file__))
 METATOMIC_TORCH = os.path.realpath(os.path.join(ROOT, "..", "metatomic_torch"))
 
-METATOMIC_TORCHSIM_VERSION = "0.1.1"
+METATOMIC_TORCHSIM_VERSION = "0.1.2"
 
 
 class sdist_generate_data(sdist):
@@ -105,7 +105,10 @@ def create_version_number(version):
     with open(os.path.join(ROOT, "AUTHORS")) as fd:
         authors = fd.read().splitlines()
 
-    install_requires = ["torch-sim-atomistic >=0.5,<0.6"]
+    install_requires = [
+        "torch-sim-atomistic >=0.5",
+        "vesin >=0.5.5,<0.6",
+    ]
 
     # when packaging a sdist for release, we should never use local dependencies
     METATOMIC_NO_LOCAL_DEPS = os.environ.get("METATOMIC_NO_LOCAL_DEPS", "0") == "1"
diff --git a/tox.ini b/tox.ini
index 7774419c..8c623b50 100644
--- a/tox.ini
+++ b/tox.ini
@@ -127,7 +127,7 @@ deps =
     torch=={env:METATOMIC_TESTS_TORCH_VERSION:2.11}.*
     numpy {env:METATOMIC_TESTS_NUMPY_VERSION_PIN}
 
-    vesin
+    vesin >=0.5.5,<0.6
     ase
 
 changedir = python/metatomic_torch
@@ -151,7 +151,7 @@ deps =
     torch=={env:METATOMIC_TESTS_TORCH_VERSION:2.11}.*
     numpy {env:METATOMIC_TESTS_NUMPY_VERSION_PIN}
 
-    vesin
+    vesin >=0.5.5,<0.6
     ase
 
 commands =
@@ -176,7 +176,7 @@ deps =
     numpy {env:METATOMIC_TESTS_NUMPY_VERSION_PIN}
 
     ase
-    vesin
+    vesin >=0.5.5,<0.6
     nvalchemi-toolkit-ops >=0.3.0,<0.4
 
     # for symmetrized calculator
@@ -208,7 +208,7 @@ deps =
     torch=={env:METATOMIC_TESTS_TORCH_VERSION:2.11}.*
     numpy {env:METATOMIC_TESTS_NUMPY_VERSION_PIN}
 
-    vesin
+    vesin >=0.5.5,<0.6
     ase
     torch-sim-atomistic
     nvalchemi-toolkit-ops >=0.3.0,<0.4
@@ -280,7 +280,7 @@ deps =
 
     # required for examples
     ase
-    vesin
+    vesin >=0.5.5,<0.6
     chemiscope
 
 commands =