Merge branch 'main' into xtb

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
+```

drivers/README → 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. 
 
 

ipi/engine/beads.py

@@ -166,15 +166,15 @@ 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:
            A Beads object with the first newP q, p, m and names arrays as the original.
         """
 
         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
 

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):

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
 

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()

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."""

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

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:

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."""