From 65591c1fca975b58bc98fb390ab2f2d07f9bf9aa Mon Sep 17 00:00:00 2001 From: Vahideh Alizadeh <82591913+V-Alizade@users.noreply.github.com> Date: Tue, 24 Sep 2024 22:42:23 +0200 Subject: [PATCH 1/8] update documentation (#377) Co-authored-by: Michele Ceriotti --- CONTRIBUTING.md | 14 +++++++------- drivers/{README => README.md} | 2 +- 2 files changed, 8 insertions(+), 8 deletions(-) rename drivers/{README => README.md} (89%) 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. From c301d1aef496fd6f291117ae459c93eaa97a2be2 Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Tue, 24 Sep 2024 16:29:37 +0200 Subject: [PATCH 2/8] Rename force, bead and cell copy to clone Small dev friendliness change - it's confusing that .copy does not really generate a copy, and can break the depend mechanism badly if beads and cell are not initialized separately. So, this PR renames to clone() to stress this does more, and improves the docstring to be less of a jerk to pretty much every developer but MC. --- ipi/engine/beads.py | 4 +-- ipi/engine/cell.py | 2 +- ipi/engine/forces.py | 46 +++++++++++++++++------------ ipi/engine/motion/al6xxx_kmc.py | 2 +- ipi/engine/motion/atomswap.py | 6 ++-- ipi/engine/motion/geop.py | 12 ++++---- ipi/engine/motion/instanton.py | 12 ++++---- ipi/engine/motion/neb.py | 12 ++++---- ipi/engine/motion/phonons.py | 6 ++-- ipi/engine/motion/planetary.py | 4 +-- ipi/engine/motion/scphonons.py | 6 ++-- ipi/engine/motion/stringmep.py | 20 ++++++------- ipi/engine/motion/vscf.py | 6 ++-- ipi/engine/properties.py | 18 +++++------ ipi/utils/instools.py | 6 ++-- tools/py/Instanton_interpolation.py | 2 +- tools/py/Instanton_postproc.py | 2 +- tools/py/neb_interpolate.py | 2 +- 18 files changed, 88 insertions(+), 80 deletions(-) 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..6be377fef 100644 --- a/ipi/engine/forces.py +++ b/ipi/engine/forces.py @@ -1006,28 +1006,36 @@ 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 ) @@ -1252,9 +1260,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 +1289,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 +1318,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..c201aac53 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.""" @@ -245,12 +245,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 From 7dc8369ad02137cfe3e3f725ab58d6c2c5bac3f6 Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Tue, 24 Sep 2024 16:31:47 +0200 Subject: [PATCH 3/8] First error fixed (Cell was not defined) --- ipi/engine/forces.py | 1 + 1 file changed, 1 insertion(+) diff --git a/ipi/engine/forces.py b/ipi/engine/forces.py index 6be377fef..a1e8a4e8f 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"] From 8e537a4cdea1ee931aac53dd9c48adaacac40fd8 Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Tue, 24 Sep 2024 20:50:52 +0200 Subject: [PATCH 4/8] Cleaner transfer_force implementation that also exposes separately the saving and loading of the force state --- ipi/engine/forces.py | 80 +++++++++++++++++++++++++++++++++----------- 1 file changed, 60 insertions(+), 20 deletions(-) diff --git a/ipi/engine/forces.py b/ipi/engine/forces.py index a1e8a4e8f..47ba87490 100644 --- a/ipi/engine/forces.py +++ b/ipi/engine/forces.py @@ -1042,30 +1042,43 @@ def clone(self, beads, cell): ) 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 + """ - 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(mself.beads._q._value), + bead_forces, + ) + ) + return force_data + + def load_state(self, state): + """ + Loads a previously-saved state of a force object. See + `transfere_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. @@ -1077,13 +1090,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)) ): From eecf9961273d0684beeb3577bb0fcff6b55d6764 Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Tue, 24 Sep 2024 20:53:36 +0200 Subject: [PATCH 5/8] oops! silly mistake should check before pushing --- ipi/engine/forces.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ipi/engine/forces.py b/ipi/engine/forces.py index 47ba87490..269c57b17 100644 --- a/ipi/engine/forces.py +++ b/ipi/engine/forces.py @@ -1122,7 +1122,7 @@ def transfer_forces(self, refforce): force.load_state(old_state) """ - self.load_state(refforce.dump_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)) From 5721dd043f8399058e03d003708fbf102c7639b3 Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Tue, 24 Sep 2024 21:10:20 +0200 Subject: [PATCH 6/8] More errors... --- ipi/engine/forces.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ipi/engine/forces.py b/ipi/engine/forces.py index 269c57b17..22f424e10 100644 --- a/ipi/engine/forces.py +++ b/ipi/engine/forces.py @@ -1058,7 +1058,7 @@ def dump_state(self): force_data.append( ( # will also need bead positions, see _load_state - deepcopy(mself.beads._q._value), + deepcopy(dstrip(mself.beads.q)), bead_forces, ) ) From 9b5c3ed9a89a72a91c0a63a734cd9e92d13e774e Mon Sep 17 00:00:00 2001 From: Mariana Rossi Date: Wed, 25 Sep 2024 15:49:56 +0200 Subject: [PATCH 7/8] Explanations and typos --- .../features/geometry_optimization/BFGS_pswater/input.xml | 2 +- ipi/engine/forces.py | 8 +++++--- ipi/engine/motion/neb.py | 3 +++ 3 files changed, 9 insertions(+), 4 deletions(-) diff --git a/examples/features/geometry_optimization/BFGS_pswater/input.xml b/examples/features/geometry_optimization/BFGS_pswater/input.xml index c4979cc42..2f1260f3f 100644 --- a/examples/features/geometry_optimization/BFGS_pswater/input.xml +++ b/examples/features/geometry_optimization/BFGS_pswater/input.xml @@ -1,4 +1,4 @@ - + [ step, potential{electronvolt}] x_centroid diff --git a/ipi/engine/forces.py b/ipi/engine/forces.py index 22f424e10..61bee0c2a 100644 --- a/ipi/engine/forces.py +++ b/ipi/engine/forces.py @@ -1044,7 +1044,9 @@ def clone(self, beads, cell): def dump_state(self): """ - Stores a (deep) copy of the internal state of a force object + Stores a (deep) copy of the internal state of a force object. + See `transfer_forces` for an explanation of the usage and + logic. """ force_data = [] @@ -1067,7 +1069,7 @@ def dump_state(self): def load_state(self, state): """ Loads a previously-saved state of a force object. See - `transfere_forces` for an explanation of the usage and + `transfer_forces` for an explanation of the usage and logic. """ @@ -1135,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 diff --git a/ipi/engine/motion/neb.py b/ipi/engine/motion/neb.py index c201aac53..4b1c5068e 100644 --- a/ipi/engine/motion/neb.py +++ b/ipi/engine/motion/neb.py @@ -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( From d645307507e0cb9cf861ee696fc375b1b92c2ed1 Mon Sep 17 00:00:00 2001 From: Mariana Rossi Date: Wed, 25 Sep 2024 15:52:48 +0200 Subject: [PATCH 8/8] Reverting input mod --- examples/features/geometry_optimization/BFGS_pswater/input.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/features/geometry_optimization/BFGS_pswater/input.xml b/examples/features/geometry_optimization/BFGS_pswater/input.xml index 2f1260f3f..c4979cc42 100644 --- a/examples/features/geometry_optimization/BFGS_pswater/input.xml +++ b/examples/features/geometry_optimization/BFGS_pswater/input.xml @@ -1,4 +1,4 @@ - + [ step, potential{electronvolt}] x_centroid