diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 5998bff70..ba3c57e6d 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -30,10 +30,10 @@ part of linting.
 In most systems, both packages can be easily installed using `pip`.
 BEFORE proceeding to a pull request, the minimal requirement is that you run
 
-::
-
-  $ make lint
-  $ make pretty
+```
+$ make lint
+$ make pretty
+```
 
 This will ensure the formatting and linting requirement are applied in the whole
 directory tree. Please resolve any warnings or errors that may appear. Your
@@ -42,6 +42,6 @@ commit will not pass the CI tests otherwise.
 For a more flexible setup, we also provide the script `i-pi-style`, for
 which instructions can be obtained by typing
 
-::
-
-  $ i-pi-style -h
+```
+$ i-pi-style -h
+```
\ No newline at end of file
diff --git a/drivers/README b/drivers/README.md
similarity index 89%
rename from drivers/README
rename to drivers/README.md
index dad690afe..0c8f9f04a 100644
--- a/drivers/README
+++ b/drivers/README.md
@@ -21,7 +21,7 @@ py
 --
 
 A Python version of the driver. Some simple potentials, and interfaces
-with codes that have a Python API can be found in the `forcefields/` 
+with codes that have a Python API can be found in the `pes/` 
 folder. 
 
 
diff --git a/ipi/engine/beads.py b/ipi/engine/beads.py
index 6249a1147..92b88c4d4 100644
--- a/ipi/engine/beads.py
+++ b/ipi/engine/beads.py
@@ -166,7 +166,7 @@ def resize(self, natoms, nbeads):
             dependencies=[b._kstress for b in self._blist],
         )
 
-    def copy(self, nbeads=-1):
+    def clone(self, nbeads=-1):
         """Creates a new beads object with newP <= P beads from the original.
 
         Returns:
@@ -174,7 +174,7 @@ def copy(self, nbeads=-1):
         """
 
         if nbeads > self.nbeads:
-            raise ValueError("Cannot copy to an object with larger number of beads")
+            raise ValueError("Cannot clone to an object with larger number of beads")
         elif nbeads == -1:
             nbeads = self.nbeads
 
diff --git a/ipi/engine/cell.py b/ipi/engine/cell.py
index 372923afa..ca69a20d2 100644
--- a/ipi/engine/cell.py
+++ b/ipi/engine/cell.py
@@ -51,7 +51,7 @@ def __init__(self, h=None):
         )
         self._V = depend_value(name="V", func=self.get_volume, dependencies=[self._h])
 
-    def copy(self):
+    def clone(self):
         return Cell(dstrip(self.h).copy())
 
     def get_ih(self):
diff --git a/ipi/engine/forces.py b/ipi/engine/forces.py
index cb0a99992..61bee0c2a 100644
--- a/ipi/engine/forces.py
+++ b/ipi/engine/forces.py
@@ -24,6 +24,7 @@
 from ipi.utils.depend import *
 from ipi.utils.nmtransform import nm_rescale
 from ipi.engine.beads import Beads
+from ipi.engine.cell import Cell
 
 
 __all__ = ["Forces", "ForceComponent"]
@@ -1006,57 +1007,80 @@ def make_rpc(rpc, beads):
             self._pots.add_dependency(fc._weight)
             self._virs.add_dependency(fc._weight)
 
-    def copy(self, beads=None, cell=None):
-        """Returns a copy of this force object that can be used to compute forces,
-        e.g. for use in internal loops of geometry optimizers, or for property
+    def clone(self, beads, cell):
+        """Duplicates the force object, so that it can be used to compute forces
+        using the same forcefields, but with a separate beads and cell object
+        (that can e.g. be created by cloning existing beads and cell objects).
+        This is useful whenever one wants to compute the same forces and energies
+        without touching the physical system, e.g. for use in internal loops of
+        geometry optimizers, when attempting Monte Carlo moves or for property
         calculation.
 
         Args:
-           beads: Optionally, bind this to a different beads object than the one
-              this Forces is currently bound
-           cell: Optionally, bind this to a different cell object
+           beads: Beads object that the clone should be bound to
+           cell: Cell object that the clone should be bound to
 
         Returns: The copy of the Forces object
         """
 
         if not self.bound:
             raise ValueError("Cannot copy a forces object that has not yet been bound.")
+        if not type(beads) is Beads:
+            raise ValueError(
+                "The 'beads' argument must be provided and be a Beads object."
+            )
+        if not type(cell) is Cell:
+            raise ValueError(
+                "The 'cell' argument must be provided and be a Cell object."
+            )
+
         nforce = Forces()
         nbeads = beads
-        if nbeads is None:
-            nbeads = self.beads
         ncell = cell
-        if cell is None:
-            ncell = self.cell
         nforce.bind(
             nbeads, ncell, self.fcomp, self.ff, self.open_paths, self.output_maker
         )
         return nforce
 
-    def transfer_forces(self, refforce):
-        """Low-level function copying over the value of a second force object,
-        triggering updates but un-tainting this force depends themselves.
+    def dump_state(self):
+        """
+        Stores a (deep) copy of the internal state of a force object.
+        See `transfer_forces` for an explanation of the usage and
+        logic.
+        """
 
-        We have noted that in some corner cases it is necessary to copy only
-        the values of updated forces rather than the full depend object, in order to
-        avoid triggering a repeated call to the client code that is potentially
-        very costly. This happens routinely in geometry relaxation routines, for example.
+        force_data = []
+        for k in range(len(self.mforces)):
+            mself = self.mforces[k]
+            # access the actual data values
+
+            bead_forces = []
+            for b in range(mself.nbeads):
+                bead_forces.append(deepcopy(mself._forces[b]._ufvx._value))
+            force_data.append(
+                (
+                    # will also need bead positions, see _load_state
+                    deepcopy(dstrip(mself.beads.q)),
+                    bead_forces,
+                )
+            )
+        return force_data
+
+    def load_state(self, state):
+        """
+        Loads a previously-saved state of a force object. See
+        `transfer_forces` for an explanation of the usage and
+        logic.
         """
 
-        if len(self.mforces) != len(refforce.mforces):
+        if len(state) != len(self.mforces):
             raise ValueError(
-                "Cannot copy forces between objects with different numbers of components"
+                "Cannot load force state from an force with a different numbers of components"
             )
 
         for k in range(len(self.mforces)):
-            mreff = refforce.mforces[k]
             mself = self.mforces[k]
-            if mreff.nbeads != mself.nbeads:
-                raise ValueError(
-                    "Cannot copy forces between objects with different numbers of beads for the "
-                    + str(k)
-                    + "th component"
-                )
+            q, fb = state[k]
 
             # this is VERY subtle. beads in this force component are
             # obtained as a contraction, and so are computed automatically.
@@ -1068,13 +1092,40 @@ def transfer_forces(self, refforce):
             # the value of the contracted bead, so that it's marked as NOT
             # tainted - it should not be as it's an internal of the force and
             # therefore get copied
-            mself.beads._q.set(mreff.beads.q, manual=False)
-            for b in range(mself.nbeads):
-                mself._forces[b]._ufvx.set(
-                    deepcopy(mreff._forces[b]._ufvx._value), manual=False
+            mself.beads._q.set(q, manual=False)
+            if len(fb) != mself.nbeads:
+                raise ValueError(
+                    "Cannot copy forces between objects with different numbers of beads for the "
+                    + str(k)
+                    + "th component"
                 )
+
+            for b in range(mself.nbeads):
+                mself._forces[b]._ufvx.set(fb[b], manual=False)
                 mself._forces[b]._ufvx.taint(taintme=False)
 
+    def transfer_forces(self, refforce):
+        """Low-level function copying over the value of a second force object,
+        triggering updates but un-tainting this force depends themselves.
+
+        We have noted that in some corner cases it is necessary to copy only
+        the values of updated forces rather than the full depend object, in order to
+        avoid triggering a repeated call to the client code that is potentially
+        very costly. This happens routinely in geometry relaxation routines,
+        for example.
+
+        The transfer can also be implemented in a delayed way, by using
+        separately set_state and dump_state. For instance, one can
+
+        old_state = force.dump_state()
+
+        do-something-that-changes-force
+
+        force.load_state(old_state)
+        """
+
+        self.load_state(refforce.dump_state())
+
     def transfer_forces_manual(
         self, new_q, new_v, new_forces, new_x=None, vir=np.zeros((3, 3))
     ):
@@ -1086,7 +1137,7 @@ def transfer_forces_manual(
           - new_f list of length equal to number of force type, containing the beads forces
           - new_x list of length equal to number of force type, containing the beads extras
         """
-        msg = "Unconsistent dimensions inside transfer_forces_manual"
+        msg = "Inconsistent dimensions inside transfer_forces_manual"
         assert len(self.mforces) == len(new_q), msg
         assert len(self.mforces) == len(new_v), msg
         assert len(self.mforces) == len(new_forces), msg
@@ -1252,9 +1303,9 @@ def forcesvirs_4th_order(self, index):
             if self.alpha == 0:
                 # we use an aux force evaluator with half the number of beads.
                 if self.dforces is None:
-                    self.dbeads = self.beads.copy(self.nbeads // 2)
-                    self.dcell = self.cell.copy()
-                    self.dforces = self.copy(self.dbeads, self.dcell)
+                    self.dbeads = self.beads.clone(self.nbeads // 2)
+                    self.dcell = self.cell.clone()
+                    self.dforces = self.clone(self.dbeads, self.dcell)
 
                 self.dcell.h = self.cell.h
 
@@ -1281,9 +1332,9 @@ def forcesvirs_4th_order(self, index):
             else:
                 # we use an aux force evaluator with the same number of beads.
                 if self.dforces is None:
-                    self.dbeads = self.beads.copy()
-                    self.dcell = self.cell.copy()
-                    self.dforces = self.copy(self.dbeads, self.dcell)
+                    self.dbeads = self.beads.clone()
+                    self.dcell = self.cell.clone()
+                    self.dforces = self.clone(self.dbeads, self.dcell)
 
                 self.dcell.h = self.cell.h
 
@@ -1310,9 +1361,9 @@ def forcesvirs_4th_order(self, index):
         if self.mforces[index].epsilon < 0.0:
             # we use an aux force evaluator with the same number of beads.
             if self.dforces is None:
-                self.dbeads = self.beads.copy()
-                self.dcell = self.cell.copy()
-                self.dforces = self.copy(self.dbeads, self.dcell)
+                self.dbeads = self.beads.clone()
+                self.dcell = self.cell.clone()
+                self.dforces = self.clone(self.dbeads, self.dcell)
 
             self.dcell.h = self.cell.h
 
diff --git a/ipi/engine/motion/al6xxx_kmc.py b/ipi/engine/motion/al6xxx_kmc.py
index 10ab09570..34f73dbf2 100755
--- a/ipi/engine/motion/al6xxx_kmc.py
+++ b/ipi/engine/motion/al6xxx_kmc.py
@@ -290,7 +290,7 @@ def bind(self, ens, beads, nm, cell, bforce, prng, omaker):
         self.dens = [None] * self.neval
         self.dbias = [None] * self.neval
         for i in range(self.neval):
-            self.dbeads[i] = beads.copy()
+            self.dbeads[i] = beads.clone()
             self.dforces[i] = bforce.copy(self.dbeads[i], self.dcell)
             self.dnm[i] = nm.copy()
             self.dens[i] = ens.copy()
diff --git a/ipi/engine/motion/atomswap.py b/ipi/engine/motion/atomswap.py
index 5ddcf0a54..99f62ab82 100644
--- a/ipi/engine/motion/atomswap.py
+++ b/ipi/engine/motion/atomswap.py
@@ -64,9 +64,9 @@ def bind(self, ens, beads, cell, bforce, nm, prng, omaker):
 
         super(AtomSwap, self).bind(ens, beads, cell, bforce, nm, prng, omaker)
         self.ensemble.add_econs(self._ealc)
-        self.dbeads = self.beads.copy()
-        self.dcell = self.cell.copy()
-        self.dforces = self.forces.copy(self.dbeads, self.dcell)
+        self.dbeads = self.beads.clone()
+        self.dcell = self.cell.clone()
+        self.dforces = self.forces.clone(self.dbeads, self.dcell)
 
     def AXlist(self, atomtype):
         """This compile a list of atoms ready for exchanges."""
diff --git a/ipi/engine/motion/geop.py b/ipi/engine/motion/geop.py
index 0e068a2e6..9847b868e 100755
--- a/ipi/engine/motion/geop.py
+++ b/ipi/engine/motion/geop.py
@@ -195,9 +195,9 @@ def __init__(self):
         self.fcount = 0
 
     def bind(self, dumop):
-        self.dbeads = dumop.beads.copy()
-        self.dcell = dumop.cell.copy()
-        self.dforces = dumop.forces.copy(self.dbeads, self.dcell)
+        self.dbeads = dumop.beads.clone()
+        self.dcell = dumop.cell.clone()
+        self.dforces = dumop.forces.clone(self.dbeads, self.dcell)
 
         self.fixatoms_mask = np.ones(
             3 * dumop.beads.natoms, dtype=bool
@@ -249,9 +249,9 @@ def __init__(self):
         pass
 
     def bind(self, dumop):
-        self.dbeads = dumop.beads.copy()
-        self.dcell = dumop.cell.copy()
-        self.dforces = dumop.forces.copy(self.dbeads, self.dcell)
+        self.dbeads = dumop.beads.clone()
+        self.dcell = dumop.cell.clone()
+        self.dforces = dumop.forces.clone(self.dbeads, self.dcell)
 
         self.fixatoms_mask = np.ones(
             3 * dumop.beads.natoms, dtype=bool
diff --git a/ipi/engine/motion/instanton.py b/ipi/engine/motion/instanton.py
index 23f74467e..d7eb2e935 100644
--- a/ipi/engine/motion/instanton.py
+++ b/ipi/engine/motion/instanton.py
@@ -268,9 +268,9 @@ def __init__(self):
         pass
 
     def bind(self, mapper):
-        self.dbeads = mapper.beads.copy()
-        self.dcell = mapper.cell.copy()
-        self.dforces = mapper.forces.copy(self.dbeads, self.dcell)
+        self.dbeads = mapper.beads.clone()
+        self.dcell = mapper.cell.clone()
+        self.dforces = mapper.forces.clone(self.dbeads, self.dcell)
 
         # self.nm = mapper.nm
         # self.rp_factor = mapper.rp_factor
@@ -355,8 +355,8 @@ def interpolation(self, full_q, full_mspath, get_all_info=False):
         reduced_b.m[:] = self.dbeads.m
         reduced_b.names[:] = self.dbeads.names
 
-        reduced_cell = self.dcell.copy()
-        reduced_forces = self.dforces.copy(reduced_b, reduced_cell)
+        reduced_cell = self.dcell.clone()
+        reduced_forces = self.dforces.clone(reduced_b, reduced_cell)
 
         # Evaluate energy and forces (and maybe friction)
         rpots = reduced_forces.pots  # reduced energy
@@ -690,7 +690,7 @@ def bind(self, mapper):
         self.temp = mapper.temp
         self.fix = mapper.fix
         self.coef = mapper.coef
-        self.dbeads = mapper.beads.copy()
+        self.dbeads = mapper.beads.clone()
         # self.nm = mapper.nm
         # self.rp_factor = mapper.rp_factor
         if self.dbeads.nbeads > 1:
diff --git a/ipi/engine/motion/neb.py b/ipi/engine/motion/neb.py
index c96ec8e99..4b1c5068e 100644
--- a/ipi/engine/motion/neb.py
+++ b/ipi/engine/motion/neb.py
@@ -51,9 +51,9 @@ def bind(self, ens):
         # In principle, there is no need in dforces within the Mapper,
         # BUT dbeads are needed to calculate tangents for the endpoints,
         # and dforces are needed outside the Mapper to construct the "main" forces.
-        self.dbeads = ens.beads.copy()
-        self.dcell = ens.cell.copy()
-        self.dforces = ens.forces.copy(self.dbeads, self.dcell)
+        self.dbeads = ens.beads.clone()
+        self.dcell = ens.cell.clone()
+        self.dforces = ens.forces.clone(self.dbeads, self.dcell)
         self.fixatoms = ens.fixatoms.copy()
 
         # Mask to exclude fixed atoms from 3N-arrays
@@ -66,7 +66,7 @@ def bind(self, ens):
         # Create reduced bead and force object
         self.rbeads = Beads(ens.beads.natoms, ens.beads.nbeads - 2)
         self.rbeads.q[:] = ens.beads.q[1:-1]
-        self.rforces = ens.forces.copy(self.rbeads, self.dcell)
+        self.rforces = ens.forces.clone(self.rbeads, self.dcell)
 
     def __call__(self, x):
         """Returns the potential for all beads and the gradient."""
@@ -91,6 +91,9 @@ def __call__(self, x):
             self.allpots = dstrip(self.dforces.pots).copy()
             # We want to be greedy about force calls,
             # so we transfer from full beads to the reduced ones.
+            # MR note, 2024: The call below is overly convoluted and
+            # likely unnecessary. It should be possible to dump values
+            # now, or use usual transfer_forces. Needs deep revision.
             tmp_f = self.dforces.f.copy()[1:-1]
             tmp_v = self.allpots.copy()[1:-1]
             self.rforces.transfer_forces_manual(
@@ -245,12 +248,12 @@ def bind(self, ens):
 
         # Reduced Beads object is needed to calculate only required beads.
         self.rbeads = Beads(ens.beads.natoms, 1)
-        self.rcell = ens.cell.copy()
+        self.rcell = ens.cell.clone()
         # Coords of the bead before and after the climbing one.
         self.q_prev = np.zeros(3 * (ens.beads.natoms - len(ens.fixatoms)))
         self.q_next = np.zeros(3 * (ens.beads.natoms - len(ens.fixatoms)))
         # Make reduced forces dependent on reduced beads
-        self.rforces = ens.forces.copy(self.rbeads, self.rcell)
+        self.rforces = ens.forces.clone(self.rbeads, self.rcell)
 
     def __call__(self, x):
         """Returns climbing force for a climbing image."""
diff --git a/ipi/engine/motion/phonons.py b/ipi/engine/motion/phonons.py
index 6de0b2ccb..5e1cd35e6 100644
--- a/ipi/engine/motion/phonons.py
+++ b/ipi/engine/motion/phonons.py
@@ -98,9 +98,9 @@ def bind(self, ens, beads, nm, cell, bforce, prng, omaker):
         self.m = dstrip(self.beads.m)
         self.phcalc.bind(self)
 
-        self.dbeads = self.beads.copy()
-        self.dcell = self.cell.copy()
-        self.dforces = self.forces.copy(self.dbeads, self.dcell)
+        self.dbeads = self.beads.clone()
+        self.dcell = self.cell.clone()
+        self.dforces = self.forces.clone(self.dbeads, self.dcell)
 
     def step(self, step=None):
         """Executes one step of phonon computation."""
diff --git a/ipi/engine/motion/planetary.py b/ipi/engine/motion/planetary.py
index a6bf735fc..58ad1a553 100644
--- a/ipi/engine/motion/planetary.py
+++ b/ipi/engine/motion/planetary.py
@@ -115,8 +115,8 @@ def bind(self, ens, beads, nm, cell, bforce, prng, omaker):
         self.basenm = nm
 
         # copies of all of the helper classes that are needed to bind the ccdyn object
-        self.dbeads = beads.copy(nbeads=self.nbeads)
-        self.dcell = cell.copy()
+        self.dbeads = beads.clone(nbeads=self.nbeads)
+        self.dcell = cell.clone()
         self.dforces = bforce.copy(self.dbeads, self.dcell)
 
         # options for NM propagation - hardcoded frequencies unless using a GLE thermo
diff --git a/ipi/engine/motion/scphonons.py b/ipi/engine/motion/scphonons.py
index ce2817a86..9c6cdad00 100644
--- a/ipi/engine/motion/scphonons.py
+++ b/ipi/engine/motion/scphonons.py
@@ -132,9 +132,9 @@ def bind(self, ens, beads, nm, cell, bforce, prng, omaker):
 
         # Creates duplicate classes to simplify computation of forces.
         self.dof = 3 * self.beads.natoms
-        self.dbeads = self.beads.copy(nbeads=self.nparallel)
-        self.dcell = self.cell.copy()
-        self.dforces = self.forces.copy(self.dbeads, self.dcell)
+        self.dbeads = self.beads.clone(nbeads=self.nparallel)
+        self.dcell = self.cell.clone()
+        self.dforces = self.forces.clone(self.dbeads, self.dcell)
 
         # Sets temperature.
         self.temp = self.ensemble.temp
diff --git a/ipi/engine/motion/stringmep.py b/ipi/engine/motion/stringmep.py
index 256bc5675..3a1fc9311 100644
--- a/ipi/engine/motion/stringmep.py
+++ b/ipi/engine/motion/stringmep.py
@@ -61,9 +61,9 @@ def bind(self, ens):
         # In principle, there is no need in dforces within the Mapper,
         # BUT dbeads are needed to calculate tangents for the endpoints,
         # and dforces are needed outside the Mapper to construct the "main" forces.
-        self.dbeads = ens.beads.copy()
-        self.dcell = ens.cell.copy()
-        self.dforces = ens.forces.copy(self.dbeads, self.dcell)
+        self.dbeads = ens.beads.clone()
+        self.dcell = ens.cell.clone()
+        self.dforces = ens.forces.clone(self.dbeads, self.dcell)
         self.fixatoms = ens.fixatoms.copy()
 
         # Mask to exclude fixed atoms from 3N-arrays
@@ -76,7 +76,7 @@ def bind(self, ens):
         # Create reduced bead and force object
         self.rbeads = Beads(ens.beads.natoms, ens.beads.nbeads - 2)
         self.rbeads.q[:] = ens.beads.q[1:-1]
-        self.rforces = ens.forces.copy(self.rbeads, self.dcell)
+        self.rforces = ens.forces.clone(self.rbeads, self.dcell)
 
     def __call__(self, x):
         """Returns the potential for all beads and the gradient."""
@@ -130,9 +130,9 @@ def bind(self, ens):
         # In principle, there is no need in dforces within the Mapper,
         # BUT dbeads are needed to calculate tangents for the endpoints,
         # and dforces are needed outside the Mapper to construct the "main" forces.
-        self.dbeads = ens.beads.copy()
-        self.dcell = ens.cell.copy()
-        self.dforces = ens.forces.copy(self.dbeads, self.dcell)
+        self.dbeads = ens.beads.clone()
+        self.dcell = ens.cell.clone()
+        self.dforces = ens.forces.clone(self.dbeads, self.dcell)
         self.fixatoms = ens.fixatoms.copy()
 
         # Mask to exclude fixed atoms from 3N-arrays
@@ -145,7 +145,7 @@ def bind(self, ens):
         # Create reduced bead and force object
         self.rbeads = Beads(ens.beads.natoms, ens.beads.nbeads - 2)
         self.rbeads.q[:] = ens.beads.q[1:-1]
-        self.rforces = ens.forces.copy(self.rbeads, self.dcell)
+        self.rforces = ens.forces.clone(self.rbeads, self.dcell)
 
     def __call__(self, x):
         """Returns the potential for all beads and the gradient."""
@@ -199,12 +199,12 @@ def bind(self, ens):
 
         # Reduced Beads object is needed to calculate only required beads.
         self.rbeads = Beads(ens.beads.natoms, 1)
-        self.rcell = ens.cell.copy()
+        self.rcell = ens.cell.clone()
         # Coords of the bead before and after the climbing one.
         self.q_prev = np.zeros(3 * (ens.beads.natoms - len(ens.fixatoms)))
         self.q_next = np.zeros(3 * (ens.beads.natoms - len(ens.fixatoms)))
         # Make reduced forces dependent on reduced beads
-        self.rforces = ens.forces.copy(self.rbeads, self.rcell)
+        self.rforces = ens.forces.clone(self.rbeads, self.rcell)
 
     def __call__(self, x):
         """Returns climbing force for a climbing image."""
diff --git a/ipi/engine/motion/vscf.py b/ipi/engine/motion/vscf.py
index c55f674a8..effe76c27 100644
--- a/ipi/engine/motion/vscf.py
+++ b/ipi/engine/motion/vscf.py
@@ -145,9 +145,9 @@ def bind(self, ens, beads, nm, cell, bforce, prng, omaker):
         self.m = dstrip(self.beads.m)
         self.calc.bind(self)
 
-        self.dbeads = self.beads.copy(nbeads=self.nparallel)
-        self.dcell = self.cell.copy()
-        self.dforces = self.forces.copy(self.dbeads, self.dcell)
+        self.dbeads = self.beads.clone(nbeads=self.nparallel)
+        self.dcell = self.cell.clone()
+        self.dforces = self.forces.clone(self.dbeads, self.dcell)
 
     def step(self, step=None):
         """Executes one step of phonon computation."""
diff --git a/ipi/engine/properties.py b/ipi/engine/properties.py
index 6357f974f..a3ab7edc7 100644
--- a/ipi/engine/properties.py
+++ b/ipi/engine/properties.py
@@ -957,9 +957,9 @@ def bind(self, system):
         # dummy beads and forcefield objects so that we can use scaled and
         # displaced path estimators without changing the simulation bead
         # coordinates
-        self.dbeads = system.beads.copy()
-        self.dcell = system.cell.copy()
-        self.dforces = system.forces.copy(self.dbeads, self.dcell)
+        self.dbeads = system.beads.clone()
+        self.dcell = system.cell.clone()
+        self.dforces = system.forces.clone(self.dbeads, self.dcell)
         self.fqref = None
         self._threadlock = (
             system._propertylock
@@ -2989,15 +2989,15 @@ def bind(self, system):
         # dummy beads and forcefield objects so that we can use scaled and
         # displaced path estimators without changing the simulation bead
         # coordinates
-        self.dbeads = system.beads.copy()
-        self.dcell = system.cell.copy()
-        self.dforces = self.system.forces.copy(self.dbeads, self.dcell)
+        self.dbeads = system.beads.clone()
+        self.dcell = system.cell.clone()
+        self.dforces = self.system.forces.clone(self.dbeads, self.dcell)
         self._threadlock = system._propertylock
 
         if system.beads.nbeads >= 2:
-            self.scdbeads = system.beads.copy(system.beads.nbeads // 2)
-            self.scdcell = system.cell.copy()
-            self.scdforces = self.system.forces.copy(self.scdbeads, self.scdcell)
+            self.scdbeads = system.beads.clone(system.beads.nbeads // 2)
+            self.scdcell = system.cell.clone()
+            self.scdforces = self.system.forces.clone(self.scdbeads, self.scdcell)
 
     def get_akcv(self):
         """Calculates the contribution to the kinetic energy due to each degree
diff --git a/ipi/utils/instools.py b/ipi/utils/instools.py
index a7db3e532..387ba154a 100644
--- a/ipi/utils/instools.py
+++ b/ipi/utils/instools.py
@@ -328,9 +328,9 @@ def __init__(self, fixatoms, beads, nbeads=None):
         self.mask2 = np.arange(3 * self.natoms * self.nbeads)[mask2]
 
         self.fixbeads = Beads(beads.natoms - len(fixatoms), beads.nbeads)
-        self.fixbeads.q[:] = self.get_active_vector(beads.copy().q, 1)
-        self.fixbeads.m[:] = self.get_active_vector(beads.copy().m, 0)
-        self.fixbeads.names[:] = self.get_active_vector(beads.copy().names, 0)
+        self.fixbeads.q[:] = self.get_active_vector(beads.clone().q, 1)
+        self.fixbeads.m[:] = self.get_active_vector(beads.clone().m, 0)
+        self.fixbeads.names[:] = self.get_active_vector(beads.clone().names, 0)
 
         mask3a = np.ones(9 * self.natoms, dtype=bool)
         for i in range(9):
diff --git a/tools/py/Instanton_interpolation.py b/tools/py/Instanton_interpolation.py
index 1cc9b4d61..befe7fd25 100755
--- a/tools/py/Instanton_interpolation.py
+++ b/tools/py/Instanton_interpolation.py
@@ -128,7 +128,7 @@
             print("We can't find {}".format(chk))
             sys.exit()
         cell = simulation.syslist[0].cell
-        beads = simulation.syslist[0].motion.beads.copy()
+        beads = simulation.syslist[0].motion.beads.clone()
         natoms = simulation.syslist[0].motion.beads.natoms
         nbeads = beads.nbeads
         q = beads.q
diff --git a/tools/py/Instanton_postproc.py b/tools/py/Instanton_postproc.py
index b52164e21..5dd91b384 100755
--- a/tools/py/Instanton_postproc.py
+++ b/tools/py/Instanton_postproc.py
@@ -263,7 +263,7 @@ def save_frequencies(d, nzeros, filename="freq.dat"):
 )
 
 
-beads = simulation.syslist[0].motion.beads.copy()
+beads = simulation.syslist[0].motion.beads.clone()
 m = simulation.syslist[0].motion.beads.m.copy()
 nbeads = simulation.syslist[0].motion.beads.nbeads
 natoms = simulation.syslist[0].motion.beads.natoms
diff --git a/tools/py/neb_interpolate.py b/tools/py/neb_interpolate.py
index e11ba1346..aaaa9a45b 100755
--- a/tools/py/neb_interpolate.py
+++ b/tools/py/neb_interpolate.py
@@ -358,7 +358,7 @@ def spline_resample(q, nbeads_old, nbeads_new, k=3):
             print("Error: cannot find {}.".format(input_chk))
             sys.exit()
         cell = simulation.syslist[0].cell
-        beads = simulation.syslist[0].motion.beads.copy()
+        beads = simulation.syslist[0].motion.beads.clone()
         natoms = simulation.syslist[0].motion.beads.natoms
         masses = simulation.syslist[0].motion.beads.m  # KF: not sure about this line
         nbeads_old = beads.nbeads