Improve the PLUMED interface (i-pi#411)

* Major update to how the plumed module collects hills work This avoids a lot of unnecessary computations, and one can switch off the calculation of total work to save even more. Note that this change makes calculations somewhat less robust, as the assumption mtd_update is called right after having updated the bias i-pi-side is even more important. However this was an assumption all along, so overall this is just being more consequential about it. * Added an OPES example, and fixed several issues with conserved quantity and metad/eval syncing * Added a note that PLUMED 2.10 is needed to fetch colvar values
m-stack-org · Dec 27, 2024 · 1646182 · 1646182
1 parent 9f2c98d
commit 1646182
Show file tree

Hide file tree

Showing 18 changed files with 209 additions and 61 deletions.
diff --git a/demos/ensemble-deltamu/README.md b/demos/ensemble-deltamu/README.md
@@ -130,7 +130,7 @@ The relevant section in the `input.xml` file is
 ```
   <ffplumed name="plumed">
       <file mode="xyz">init.xyz</file>
-      <plumeddat>plumed.dat </plumeddat>
+      <plumed_dat>plumed.dat </plumed_dat>
   </ffplumed>
 ```
 

diff --git a/demos/ensemble-deltamu/input_n2p2.xml b/demos/ensemble-deltamu/input_n2p2.xml
@@ -26,7 +26,7 @@
   </ffcommittee>
   <ffplumed name="plumed">                                                    
       <file mode="xyz">init_ase.xyz</file>
-      <plumeddat>plumed.dat </plumeddat>                              
+      <plumed_dat>plumed.dat </plumed_dat>                              
   </ffplumed> 
    <system>
       <initialize nbeads="1">

diff --git a/demos/ensemble-deltamu/input_rascaline.xml b/demos/ensemble-deltamu/input_rascaline.xml
@@ -17,7 +17,7 @@
   </ffcommittee>
   <ffplumed name="plumed">                                                    
       <file mode="xyz">init_ase.xyz</file>
-      <plumeddat>plumed.dat </plumeddat>
+      <plumed_dat>plumed.dat </plumed_dat>
   </ffplumed> 
    <system>
       <initialize nbeads="1">

diff --git a/examples/features/metadynamics/opes_metadynamics_zundel/h5o2+.xyz b/examples/features/metadynamics/opes_metadynamics_zundel/h5o2+.xyz
@@ -0,0 +1 @@
+../../../init_files/h5o2+.xyz
diff --git a/examples/features/metadynamics/opes_metadynamics_zundel/h5o2.dms4B.coeff.com.dat b/examples/features/metadynamics/opes_metadynamics_zundel/h5o2.dms4B.coeff.com.dat
@@ -0,0 +1 @@
+../../../../drivers/f90/pes/h5o2.dms4B.coeff.com.dat
diff --git a/examples/features/metadynamics/opes_metadynamics_zundel/h5o2.pes4B.coeff.dat b/examples/features/metadynamics/opes_metadynamics_zundel/h5o2.pes4B.coeff.dat
@@ -0,0 +1 @@
+../../../../drivers/f90/pes/h5o2.pes4B.coeff.dat
diff --git a/examples/features/metadynamics/opes_metadynamics_zundel/input.xml b/examples/features/metadynamics/opes_metadynamics_zundel/input.xml
@@ -0,0 +1,77 @@
+<simulation>
+    <ffsocket mode='unix' name='driver'>
+         <latency>  1.00000000e-02</latency>
+         <slots>4</slots>
+         <port>20614</port>
+         <timeout>  6.00000000e+02</timeout>
+         <address>zundel</address>
+   </ffsocket>
+   <ffplumed name="plumed">
+         <file mode="xyz">./h5o2+.xyz</file>
+         <plumed_dat> plumed/plumed.dat </plumed_dat>
+         <plumed_extras> [doo, dc, opes.bias ] </plumed_extras>
+   </ffplumed>
+   <total_steps>400</total_steps>
+   <output prefix="data">
+      <trajectory stride="40" filename="pos" cell_units="angstrom">positions{angstrom}</trajectory>
+      <trajectory stride="20" filename="xc" format="xyz">x_centroid{angstrom}</trajectory>
+      <trajectory stride="20" filename="colvar" bead="0" extra_type="doo,dc,opes.bias"> extras_bias </trajectory>
+      <properties stride="2"> [ step, time, conserved, temperature{kelvin}, kinetic_cv, 
+            potential, kinetic_cv(H), kinetic_cv(O), ensemble_bias ] </properties>
+   </output>
+   <prng>
+      <seed>18885</seed>
+   </prng>
+   <system>
+      <forces> 
+          <force forcefield="driver"></force>
+      </forces>
+      <initialize nbeads="8">
+        <file mode="xyz">./h5o2+.xyz</file>
+        <cell>
+         [   25.29166, 0, 0, 0, 25.29166, 0, 0, 0, 25.29166 ]
+        </cell>
+      </initialize>
+      <ensemble>
+         <temperature units="kelvin"> 300.0 </temperature>
+         <bias>
+           <force forcefield="plumed" nbeads="1">
+             <interpolate_extras> [ doo, dc, opes.bias ] </interpolate_extras>
+           </force>
+         </bias>
+      </ensemble>
+      <motion mode="dynamics">
+        <dynamics mode="nvt">
+         <timestep units="femtosecond"> 0.5 </timestep>
+<!--
+# Generated at http://cosmo-epfl.github.io/gle4md
+# Please cite:
+# M. Ceriotti, G. Bussi and M. Parrinello, J. Chem. Theory Comput. 6, 1170 (2010) 
+# M. Ceriotti, G. Bussi and M. Parrinello, Phys. Rev. Lett. 102, 020601 (2009)
+# Smart-sampling GLE. Enforces efficient sampling, focussing the effort on the slowest mode 
+# accessible by the simulation. Generated from the parameter file
+# library/smart/smart-0.5_6-2.a, 
+# and shifted so that they are effective to sample optimally 
+# a time scale of t_opt=1 picoseconds,
+# and do as well as possible upt to a cutoff frequency of 
+# νmax=100 THz [3336 cm^-1]
+-->
+  <thermostat mode='gle'>
+    <A shape='(7,7)'>
+      [   8.191023526179e-4,    8.328506066524e-3,    1.657771834013e-3,    9.736989925341e-4,    2.841803794895e-4,   -3.176846864198e-5,   -2.967010478210e-4, 
+  -8.389856546341e-4,    2.405526974742e-2,   -1.507872374848e-2,    2.589784240185e-3,    1.516783633362e-3,   -5.958833418565e-4,    4.198422349789e-4, 
+   7.798710586406e-4,    1.507872374848e-2,    8.569039501219e-3,    6.001000899602e-3,    1.062029383877e-3,    1.093939147968e-3,   -2.661575532976e-3, 
+  -9.676783161546e-4,   -2.589784240185e-3,   -6.001000899602e-3,    2.680459336535e-5,   -5.214694469742e-5,    4.231304910751e-4,   -2.104894919743e-5, 
+  -2.841997149166e-4,   -1.516783633362e-3,   -1.062029383877e-3,    5.214694469742e-5,    1.433903506353e-9,   -4.241574212449e-5,    7.910178912362e-5, 
+   3.333208286893e-5,    5.958833418565e-4,   -1.093939147968e-3,   -4.231304910751e-4,    4.241574212449e-5,    2.385554468441e-8,   -3.139255482869e-5, 
+   2.967533789056e-4,   -4.198422349789e-4,    2.661575532976e-3,    2.104894919743e-5,   -7.910178912362e-5,    3.139255482869e-5,   2.432567259684e-11
+     ]
+    </A>
+  </thermostat>
+        </dynamics>
+      </motion>
+  </system>
+  <smotion mode="metad">
+     <metad> <metaff> [ plumed ] </metaff> </metad>
+  </smotion>
+</simulation>
diff --git a/examples/features/metadynamics/opes_metadynamics_zundel/plumed/plumed.dat b/examples/features/metadynamics/opes_metadynamics_zundel/plumed/plumed.dat
@@ -0,0 +1,21 @@
+# default units are LENGTH=nm ENERGY=kJ/mol TIME=ps
+doo: DISTANCE ATOMS=1,2 
+co1: DISTANCES GROUPA=1 GROUPB=3-7 LESS_THAN={RATIONAL R_0=0.14}
+co2: DISTANCES GROUPA=2 GROUPB=3-7 LESS_THAN={RATIONAL R_0=0.14}
+dc: COMBINE ARG=co1.lessthan,co2.lessthan COEFFICIENTS=1,-1 PERIODIC=NO
+OPES_METAD ...
+LABEL=opes
+ARG=doo,dc
+PACE=50
+BARRIER=50
+TEMP=300
+FILE=plumed/KERNELS
+STATE_RFILE=plumed/STATE.latest
+STATE_WFILE=plumed/STATE
+STATE_WSTRIDE=1*10000
+STORE_STATES
+... OPES_METAD
+uwall: UPPER_WALLS ARG=doo AT=0.4 KAPPA=250
+
+PRINT ARG=doo,co1.*,co2.*,opes.*,uwall.* STRIDE=10 FILE=plumed/COLVAR 
+FLUSH STRIDE=1
diff --git a/examples/features/metadynamics/pimd_metadynamics_zundel/input.xml b/examples/features/metadynamics/pimd_metadynamics_zundel/input.xml
@@ -8,14 +8,14 @@
    </ffsocket>
    <ffplumed name="plumed">
          <file mode="xyz">./h5o2+.xyz</file>
-         <plumeddat> plumed/plumed.dat </plumeddat>
-         <plumed_extras> [doo, co1.lessthan, co2.lessthan, mtd.bias ] </plumed_extras>
+         <plumed_dat> plumed/plumed.dat </plumed_dat>
+         <plumed_extras> [doo, dc, mtd.bias ] </plumed_extras>
    </ffplumed>
    <total_steps>400</total_steps>
    <output prefix="data">
       <trajectory stride="40" filename="pos" cell_units="angstrom">positions{angstrom}</trajectory>
       <trajectory stride="20" filename="xc" format="xyz">x_centroid{angstrom}</trajectory>
-      <trajectory stride="20" filename="colvar" bead="0" extra_type="doo,co1.lessthan,co2.lessthan,mtd.bias"> extras_bias </trajectory>
+      <trajectory stride="20" filename="colvar" bead="0" extra_type="doo,dc,mtd.bias"> extras_bias </trajectory>
       <properties stride="2"> [ step, time, conserved, temperature{kelvin}, kinetic_cv, 
             potential, kinetic_cv(H), kinetic_cv(O), ensemble_bias ] </properties>
    </output>
@@ -36,7 +36,7 @@
          <temperature units="kelvin"> 300.0 </temperature>
          <bias>
            <force forcefield="plumed" nbeads="1">
-             <interpolate_extras> [ doo, mtd.bias ] </interpolate_extras>
+             <interpolate_extras> [ doo, dc, mtd.bias ] </interpolate_extras>
            </force>
          </bias>
       </ensemble>

diff --git a/examples/features/metadynamics/remd_metadynamics_zundel/input.xml b/examples/features/metadynamics/remd_metadynamics_zundel/input.xml
@@ -8,15 +8,15 @@
    </ffsocket>
    <ffplumed name="plumed-200K">
          <file mode="xyz">./h5o2+.xyz</file>
-         <plumeddat> plumed_200K/plumed.dat </plumeddat>
+         <plumed_dat> plumed_200K/plumed.dat </plumed_dat>
    </ffplumed>
    <ffplumed name="plumed-250K">
          <file mode="xyz">./h5o2+.xyz</file>
-         <plumeddat> plumed_250K/plumed.dat </plumeddat>
+         <plumed_dat> plumed_250K/plumed.dat </plumed_dat>
    </ffplumed>
    <ffplumed name="plumed-320K">
          <file mode="xyz">./h5o2+.xyz</file>
-         <plumeddat> plumed_320K/plumed.dat </plumeddat>
+         <plumed_dat> plumed_320K/plumed.dat </plumed_dat>
    </ffplumed>
    <total_steps>20000</total_steps>
    <output prefix="data">

diff --git a/examples/features/metadynamics/wte_metadynamics_zundel/input.xml b/examples/features/metadynamics/wte_metadynamics_zundel/input.xml
@@ -8,7 +8,7 @@
    </ffsocket>
    <ffplumed name="plumed">
          <file mode="xyz">./h5o2+.xyz</file>
-         <plumeddat> plumed/plumed.dat </plumeddat>
+         <plumed_dat> plumed/plumed.dat </plumed_dat>
    </ffplumed>
    <total_steps>400</total_steps>
    <output prefix="data">

diff --git a/ipi/engine/forcefields.py b/ipi/engine/forcefields.py
@@ -734,8 +734,9 @@ def __init__(
         dopbc=False,
         threaded=False,
         init_file="",
-        plumeddat="",
-        plumedstep=0,
+        compute_work=True,
+        plumed_dat="",
+        plumed_step=0,
         plumed_extras=[],
     ):
         """Initialises FFPlumed.
@@ -756,9 +757,10 @@ def __init__(
             latency, offset, name, pars, dopbc=False, threaded=threaded
         )
         self.plumed = plumed.Plumed()
-        self.plumeddat = plumeddat
-        self.plumedstep = plumedstep
+        self.plumed_dat = plumed_dat
+        self.plumed_step = plumed_step
         self.plumed_extras = plumed_extras
+        self.compute_work = compute_work
         self.init_file = init_file
 
         if self.init_file.mode == "xyz":
@@ -772,7 +774,7 @@ def __init__(
         self.natoms = myatoms.natoms
         self.plumed.cmd("setRealPrecision", 8)  # i-PI uses double precision
         self.plumed.cmd("setMDEngine", "i-pi")
-        self.plumed.cmd("setPlumedDat", self.plumeddat)
+        self.plumed.cmd("setPlumedDat", self.plumed_dat)
         self.plumed.cmd("setNatoms", self.natoms)
         timeunit = 2.4188843e-05  # atomic time to ps
         self.plumed.cmd("setMDTimeUnits", timeunit)
@@ -787,7 +789,7 @@ def __init__(
             "setMDLengthUnits", 0.052917721
         )  # Pass a pointer to the conversion factor between the length unit used in your code and nm
         self.plumedrestart = False
-        if self.plumedstep > 0:
+        if self.plumed_step > 0:
             # we are restarting, signal that PLUMED should continue
             self.plumedrestart = True
             self.plumed.cmd("setRestart", 1)
@@ -810,6 +812,12 @@ def __init__(
         self.masses = dstrip(myatoms.m)
         self.lastq = np.zeros(3 * self.natoms)
         self.system_force = None  # reference to physical force calculator
+        softexit.register_function(self.softexit)
+
+    def softexit(self):
+        """Takes care of cleaning up upon softexit"""
+
+        self.plumed.finalize()
 
     def evaluate(self, r):
         """A wrapper function to call the PLUMED evaluation routines
@@ -829,7 +837,7 @@ def evaluate(self, r):
         # linking with the current value in simulations is non-trivial, as masses
         # are not expected to be the force evaluator's business, and charges are not
         # i-PI's business.
-        self.plumed.cmd("setStep", self.plumedstep)
+        self.plumed.cmd("setStep", self.plumed_step)
         self.plumed.cmd("setCharges", self.charges)
         self.plumed.cmd("setMasses", self.masses)
 
@@ -871,27 +879,41 @@ def mtd_update(self, pos, cell):
         """Makes updates to the potential that only need to be triggered
         upon completion of a time step."""
 
-        self.plumedstep += 1
-        f = np.zeros((self.natoms, 3))
-        vir = np.zeros((3, 3))
+        # NB - this assumes this is called at the end of a step,
+        # when the bias has already been computed to integrate MD
+        # unexpected behavior will happen if it's called when the
+        # bias force is not "freshly computed"
 
-        self.plumed.cmd("setStep", self.plumedstep)
-        self.plumed.cmd("setCharges", self.charges)
-        self.plumed.cmd("setMasses", self.masses)
-        rpos = pos.reshape((-1, 3))
-        self.plumed.cmd("setPositions", rpos)
-        self.plumed.cmd("setBox", cell.T.copy())
-        if self.system_force is not None:
-            f[:] = dstrip(self.system_force.f).reshape((-1, 3))
-            vir[:] = -dstrip(self.system_force.vir)
-            self.plumed.cmd("setEnergy", dstrip(self.system_force.pot))
-        self.plumed.cmd("setForces", f)
-        self.plumed.cmd("setVirial", vir)
-        self.plumed.cmd("prepareCalc")
-        self.plumed.cmd("performCalcNoUpdate")
+        self.plumed_step += 1
+
+        bias_before = np.zeros(1, float)
+        bias_after = np.zeros(1, float)
+
+        if self.compute_work:
+            self.plumed.cmd("getBias", bias_before)
+
+        # Checks that the update is called on the right position.
+        # this should be the case for most workflows - if this error
+        # is triggered and your input makes sense, the right thing to
+        # do is to perform a full plumed-side update (which will have a cost,
+        # so see if you can avoid it)
+        if np.linalg.norm(self.lastq - pos) > 1e-10:
+            raise ValueError(
+                "Metadynamics update is performed using an incorrect position"
+            )
+
+        # sets the step and does the actual update
+        self.plumed.cmd("setStep", self.plumed_step)
         self.plumed.cmd("update")
 
-        return True
+        # recompute the bias so we can compute the work
+        if self.compute_work:
+            self.plumed.cmd("performCalcNoForces")
+            self.plumed.cmd("getBias", bias_after)
+
+        work = (bias_before - bias_after).item()
+
+        return work
 
 
 class FFYaff(FFEval):

diff --git a/ipi/engine/smotion/metad.py b/ipi/engine/smotion/metad.py
@@ -71,6 +71,7 @@ def step(self, step=None):
                 k = bc.ffield
                 if k not in self.metaff:
                     continue  # only does metad for the indicated forcefield
+
                 f = s.ensemble.bias.ff[k]
                 if not hasattr(
                     f, "mtd_update"
@@ -83,16 +84,25 @@ def step(self, step=None):
                 if s.ensemble.bweights[ik] == 0:
                     continue  # do not put metad bias on biases with zero weights (useful to do remd+metad!)
 
-                meta_pot_before = s.ensemble.bias.pot
-
-                fmtd = f.mtd_update(pos=s.beads.qc, cell=s.cell.h)
-                if fmtd:  # if metadyn has updated, then we must recompute forces.
-                    # hacky but cannot think of a better way: we must manually taint *just* that component
+                # MTD is hardcoded to be applied on the centroid variable.
+                # this is the "right" thing if you need to compute kinetics based
+                # on the resulting FES
+                mtd_work = f.mtd_update(pos=s.beads.qc, cell=s.cell.h)
+
+                # updates the conserved quantity with the change in bias so that
+                # we remove the shift due to added hills
+                s.ensemble.eens += (
+                    mtd_work * s.beads.nbeads
+                )  # apply ring polymer contraction!
+
+                if mtd_work != 0:
+                    # hacky but cannot think of a better way: we must manually taint *just* that component.
+                    # we also use the fact that the bias force from a hill is zero when it's added so we
+                    # don't need changes to the forces, only to the bias
                     for fc in s.ensemble.bias.mforces:
                         if fc.ffield == k:
                             for fb in fc._forces:
-                                fb._ufvx.taint()
-                    meta_pot_after = s.ensemble.bias.pot
-                    # updates the conserved quantity with the change in bias so that
-                    # we remove the shift due to added hills
-                    s.ensemble.eens += meta_pot_before - meta_pot_after
+                                # this open-heart surgery on a depend object is fugly
+                                # but can't se a better way
+                                fb._ufvx._value[0] -= mtd_work
+                                fb._ufvx.taint(taintme=False)
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1 @@
		../../../../drivers/f90/pes/h5o2.dms4B.coeff.com.dat
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1 @@
		../../../../drivers/f90/pes/h5o2.pes4B.coeff.dat