Implement fixatoms_dof tag and variables (i-pi#394)

* Implement fixatoms_dof tag and variables * lint * lint with proper black version * bugifix * bug fix * bugfix * bugfix * update some fixatoms to fixatoms_dof still remaining * more bugfix * bugfix * lint * bugfix * forgot about hesstools, now it is updated * sort and flatten fixatoms_dof * add regtest * upload ref_vib.out missing * imporve help-string fixdof * add examples inside geop * Removing unnecessary debug print statement --------- Co-authored-by: Mariana Rossi <[email protected]>
lab-cosmo · Nov 12, 2024 · 4ae8223 · 4ae8223
1 parent ffce0d3
commit 4ae8223
Show file tree

Hide file tree

Showing 41 changed files with 839 additions and 243 deletions.
diff --git a/demos/te-pigs/1_dataset_curation/get_pigs_dataset.py b/demos/te-pigs/1_dataset_curation/get_pigs_dataset.py
@@ -11,7 +11,6 @@
     iread(centroid_force_filename),
     iread(physical_force_filename),
 ):
-
     atoms = a_pos.copy()
 
     atoms.arrays["centroid_force"] = a_cforce.arrays["f_centroid"]

diff --git a/examples/features/geometry_optimization/SD_pswater_with_fixatoms/init.xyz b/examples/features/geometry_optimization/SD_pswater_with_fixatoms/init.xyz
@@ -0,0 +1,5 @@
+3
+# CELL(abcABC):   35.23300    35.23300    35.23300    90.00000    90.00000    90.00000  Traj: positions{angstrom} Step:          11  Bead:       0 
+       O  0.00000e+00 -2.00000e-02  0.00000e+00
+       H  7.50000e-01  5.00000e-01  0.00000e+00
+       H -7.50000e-01  5.00000e-01  0.00000e+00
diff --git a/examples/features/geometry_optimization/SD_pswater_with_fixatoms/input.xml b/examples/features/geometry_optimization/SD_pswater_with_fixatoms/input.xml
@@ -0,0 +1,31 @@
+<simulation mode='static' verbosity='high'>
+  <output prefix='simulation'>
+    <properties stride='1' filename='out'>  [ step, potential ] </properties>
+    <trajectory filename='pos' stride='1'> positions </trajectory>
+  </output>
+  <total_steps> 1000 </total_steps>
+  <prng>
+    <seed> 32342 </seed>
+  </prng>
+  <ffsocket name='pswater' mode='unix' pbc='false'>
+    <address> h2o-geop </address>
+  </ffsocket>
+  <system>
+    <initialize nbeads='1'>
+      <file mode='xyz'> init.xyz </file>
+    </initialize>
+    <forces>
+      <force forcefield='pswater'> </force>
+    </forces>
+    <motion mode='minimize'>
+       <fixatoms>[0]</fixatoms>
+      <optimizer mode='sd'>
+        <tolerances>
+          <energy> 1e-5 </energy>
+          <force> 1e-5 </force>
+          <position> 1e-5 </position>
+        </tolerances>
+      </optimizer>
+    </motion>
+  </system>
+</simulation>
diff --git a/examples/features/geometry_optimization/SD_pswater_with_fixatoms/run.sh b/examples/features/geometry_optimization/SD_pswater_with_fixatoms/run.sh
@@ -0,0 +1,16 @@
+ipi=i-pi
+driver="i-pi-driver -m pswater -u -a h2o-geop"
+sleep_time=4
+
+${ipi} input.xml > log.i-pi & 
+echo "# i-PI is running"
+
+echo "# Waiting for ${sleep_time} (s) before executing driver"
+sleep ${sleep_time}
+
+${driver} > /dev/null & 
+echo "# Driver is running"
+
+wait
+
+echo "# Simulation complete"
diff --git a/examples/features/geometry_optimization/SD_pswater_with_fixdofs/init.xyz b/examples/features/geometry_optimization/SD_pswater_with_fixdofs/init.xyz
@@ -0,0 +1,5 @@
+3
+# CELL(abcABC):   35.23300    35.23300    35.23300    90.00000    90.00000    90.00000  Traj: positions{angstrom} Step:          11  Bead:       0 
+       O  0.00000e+00 -2.00000e-02  0.00000e+00
+       H  7.50000e-01  5.00000e-01  0.00000e+00
+       H -7.50000e-01  5.00000e-01  0.00000e+00
diff --git a/examples/features/geometry_optimization/SD_pswater_with_fixdofs/input.xml b/examples/features/geometry_optimization/SD_pswater_with_fixdofs/input.xml
@@ -0,0 +1,31 @@
+<simulation mode='static' verbosity='high'>
+  <output prefix='simulation'>
+    <properties stride='1' filename='out'>  [ step, potential ] </properties>
+    <trajectory filename='pos' stride='1'> positions </trajectory>
+  </output>
+  <total_steps> 1000 </total_steps>
+  <prng>
+    <seed> 32342 </seed>
+  </prng>
+  <ffsocket name='pswater' mode='unix' pbc='false'>
+    <address> h2o-geop </address>
+  </ffsocket>
+  <system>
+    <initialize nbeads='1'>
+      <file mode='xyz'> init.xyz </file>
+    </initialize>
+    <forces>
+      <force forcefield='pswater'> </force>
+    </forces>
+    <motion mode='minimize'>
+       <fixatoms_dof>[0,1,2]</fixatoms_dof>
+      <optimizer mode='sd'>
+        <tolerances>
+          <energy> 1e-5 </energy>
+          <force> 1e-5 </force>
+          <position> 1e-5 </position>
+        </tolerances>
+      </optimizer>
+    </motion>
+  </system>
+</simulation>
diff --git a/examples/features/geometry_optimization/SD_pswater_with_fixdofs/run.sh b/examples/features/geometry_optimization/SD_pswater_with_fixdofs/run.sh
@@ -0,0 +1,16 @@
+ipi=i-pi
+driver="i-pi-driver -m pswater -u -a h2o-geop"
+sleep_time=4
+
+${ipi} input.xml > log.i-pi & 
+echo "# i-PI is running"
+
+echo "# Waiting for ${sleep_time} (s) before executing driver"
+sleep ${sleep_time}
+
+${driver} > /dev/null & 
+echo "# Driver is running"
+
+wait
+
+echo "# Simulation complete"
diff --git a/ipi/engine/motion/al6xxx_kmc.py b/ipi/engine/motion/al6xxx_kmc.py
@@ -68,7 +68,7 @@ def __init__(
         thermostat=None,
         barostat=None,
         fixcom=False,
-        fixatoms=None,
+        fixatoms_dof=None,
         nmts=None,
         max_cache_len=1000,
     ):
@@ -165,7 +165,7 @@ def __init__(
 
         # the KMC step is variable and so it cannot be stored as proper timing
         self._dt = depend_value(name="dt", value=0.0)
-        self.fixatoms = np.asarray([])
+        self.fixatoms_dof = np.asarray([])
         self.fixcom = True
         self.optimizer = [None] * self.neval
         # geop should not trigger exit if there is early convergence, but just carry on.
@@ -174,7 +174,7 @@ def __init__(
         for i in range(self.neval):
             # geometry optimizer should not have *any* hystory dependence
             self.optimizer[i] = GeopMotion(
-                fixcom=fixcom, fixatoms=fixatoms, **optimizer
+                fixcom=fixcom, fixatoms_dof=fixatoms_dof, **optimizer
             )  # mode="cg", ls_options={"tolerance": 1, "iter": 20,  "step": 1e-3, "adaptive": 0.0}, tolerances={"energy": 1e-7, "force": 1e-2, "position": 1e-4}, ) #!TODO: set the geop parameters properly
 
         # dictionary of previous energy evaluations - kind of tricky to use this with the omaker thingie

diff --git a/ipi/engine/motion/alchemy.py b/ipi/engine/motion/alchemy.py
@@ -30,7 +30,7 @@ class AlchemyMC(Motion):
     """
 
     def __init__(
-        self, fixcom=False, fixatoms=None, mode=None, names=[], nxc=1, ealc=None
+        self, fixcom=False, fixatoms_dof=None, mode=None, names=[], nxc=1, ealc=None
     ):
         """Initialises a "alchemical exchange" motion object.
 
@@ -40,7 +40,7 @@ def __init__(
 
         """
 
-        super(AlchemyMC, self).__init__(fixcom=fixcom, fixatoms=fixatoms)
+        super(AlchemyMC, self).__init__(fixcom=fixcom, fixatoms_dof=fixatoms_dof)
 
         self.names = names
         self.nxc = nxc

diff --git a/ipi/engine/motion/atomswap.py b/ipi/engine/motion/atomswap.py
@@ -25,7 +25,7 @@ class AtomSwap(Motion):
     """
 
     def __init__(
-        self, fixcom=False, fixatoms=None, mode=None, names=[], nxc=1, ealc=None
+        self, fixcom=False, fixatoms_dof=None, mode=None, names=[], nxc=1, ealc=None
     ):
         """Initialises a "alchemical exchange" motion object.
 
@@ -35,7 +35,7 @@ def __init__(
 
         """
 
-        super(AtomSwap, self).__init__(fixcom=fixcom, fixatoms=fixatoms)
+        super(AtomSwap, self).__init__(fixcom=fixcom, fixatoms_dof=fixatoms_dof)
 
         self.names = names
         self.nxc = nxc

diff --git a/ipi/engine/motion/dynamics.py b/ipi/engine/motion/dynamics.py
@@ -53,7 +53,7 @@ def __init__(
         thermostat=None,
         barostat=None,
         fixcom=False,
-        fixatoms=None,
+        fixatoms_dof=None,
         nmts=None,
         efield=None,
         bec=None,
@@ -66,7 +66,7 @@ def __init__(
                 motion will be constrained or not. Defaults to False.
         """
 
-        super(Dynamics, self).__init__(fixcom=fixcom, fixatoms=fixatoms)
+        super(Dynamics, self).__init__(fixcom=fixcom, fixatoms_dof=fixatoms_dof)
 
         # initialize time step. this is the main time step that covers a full time step
         self._dt = depend_value(name="dt", value=timestep)
@@ -76,7 +76,7 @@ def __init__(
         else:
             if (
                 thermostat.__class__.__name__ is ("ThermoPILE_G" or "ThermoNMGLEG ")
-            ) and (len(fixatoms) > 0):
+            ) and (len(fixatoms_dof) > 0):
                 softexit.trigger(
                     status="bad",
                     message="!! Sorry, fixed atoms and global thermostat on the centroid not supported. Use a local thermostat. !!",
@@ -115,21 +115,17 @@ def __init__(
 
         # constraints
         self.fixcom = fixcom
-        if fixatoms is None:
-            self.fixatoms = np.zeros(0, int)
-            self.fixatoms3 = np.zeros(0, int)
+        if fixatoms_dof is None:
+            self.fixatoms_dof = np.zeros(0, int)
         else:
-            self.fixatoms = fixatoms
-            self.fixatoms3 = np.array(
-                [[3 * i, 3 * i + 1, 3 * i + 2] for i in self.fixatoms]
-            ).flatten()
+            self.fixatoms_dof = fixatoms_dof
 
     def get_fixdof(self):
         """Calculate the number of fixed degrees of freedom, required for
         temperature and pressure calculations.
         """
 
-        fixdof = len(self.fixatoms) * 3 * self.beads.nbeads
+        fixdof = len(self.fixatoms_dof) * self.beads.nbeads
         if self.fixcom:
             fixdof += 3
         return fixdof
@@ -279,8 +275,7 @@ def bind(self, motion):
         self.thermostat = motion.thermostat
         self.barostat = motion.barostat
         self.fixcom = motion.fixcom
-        self.fixatoms = motion.fixatoms
-        self.fixatoms3 = motion.fixatoms3
+        self.fixatoms_dof = motion.fixatoms_dof
         self.enstype = motion.enstype
 
         # no need to dpipe these are really just references
@@ -368,13 +363,14 @@ def pconstraints(self):
 
             self.ensemble.eens += np.sum(vcom**2) * 0.5 * Mnb  # COM kinetic energy
 
-        if len(self.fixatoms) > 0:
+        if len(self.fixatoms_dof) > 0:
             m3 = dstrip(beads.m3)
             p = dstrip(beads.p)
-            fixatoms3 = self.fixatoms3
 
-            self.ensemble.eens += 0.5 * np.sum(p[:, fixatoms3] ** 2 / m3[:, fixatoms3])
-            beads.p[:, fixatoms3] = 0.0
+            self.ensemble.eens += 0.5 * np.sum(
+                p[:, self.fixatoms_dof] ** 2 / m3[:, self.fixatoms_dof]
+            )
+            beads.p[:, self.fixatoms_dof] = 0.0
 
 
 dproperties(