From b93354b568b1289eeb5864856d2f46244c612eb5 Mon Sep 17 00:00:00 2001
From: Michele Ceriotti <michele.ceriotti@gmail.com>
Date: Fri, 6 Dec 2024 06:26:51 -0500
Subject: [PATCH] Added an OPES example, and fixed several issues with
 conserved quantity and metad/eval syncing

---
 .../opes_metadynamics_zundel/h5o2+.xyz        |  1 +
 .../h5o2.dms4B.coeff.com.dat                  |  1 +
 .../h5o2.pes4B.coeff.dat                      |  1 +
 .../opes_metadynamics_zundel/input.xml        | 77 +++++++++++++++++++
 .../plumed/plumed.dat                         | 21 +++++
 .../pimd_metadynamics_zundel/input.xml        |  6 +-
 ipi/engine/forcefields.py                     | 10 +++
 ipi/engine/smotion/metad.py                   |  8 +-
 8 files changed, 121 insertions(+), 4 deletions(-)
 create mode 120000 examples/features/metadynamics/opes_metadynamics_zundel/h5o2+.xyz
 create mode 120000 examples/features/metadynamics/opes_metadynamics_zundel/h5o2.dms4B.coeff.com.dat
 create mode 120000 examples/features/metadynamics/opes_metadynamics_zundel/h5o2.pes4B.coeff.dat
 create mode 100644 examples/features/metadynamics/opes_metadynamics_zundel/input.xml
 create mode 100644 examples/features/metadynamics/opes_metadynamics_zundel/plumed/plumed.dat

diff --git a/examples/features/metadynamics/opes_metadynamics_zundel/h5o2+.xyz b/examples/features/metadynamics/opes_metadynamics_zundel/h5o2+.xyz
new file mode 120000
index 000000000..12c06742d
--- /dev/null
+++ b/examples/features/metadynamics/opes_metadynamics_zundel/h5o2+.xyz
@@ -0,0 +1 @@
+../../../init_files/h5o2+.xyz
\ No newline at end of file
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
new file mode 120000
index 000000000..35797e5cc
--- /dev/null
+++ b/examples/features/metadynamics/opes_metadynamics_zundel/h5o2.dms4B.coeff.com.dat
@@ -0,0 +1 @@
+../../../../drivers/f90/pes/h5o2.dms4B.coeff.com.dat
\ No newline at end of file
diff --git a/examples/features/metadynamics/opes_metadynamics_zundel/h5o2.pes4B.coeff.dat b/examples/features/metadynamics/opes_metadynamics_zundel/h5o2.pes4B.coeff.dat
new file mode 120000
index 000000000..d5173136a
--- /dev/null
+++ b/examples/features/metadynamics/opes_metadynamics_zundel/h5o2.pes4B.coeff.dat
@@ -0,0 +1 @@
+../../../../drivers/f90/pes/h5o2.pes4B.coeff.dat
\ No newline at end of file
diff --git a/examples/features/metadynamics/opes_metadynamics_zundel/input.xml b/examples/features/metadynamics/opes_metadynamics_zundel/input.xml
new file mode 100644
index 000000000..4f849a516
--- /dev/null
+++ 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
new file mode 100644
index 000000000..76974a4d5
--- /dev/null
+++ 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
index 663df51c6..aa3dac6e3 100644
--- a/examples/features/metadynamics/pimd_metadynamics_zundel/input.xml
+++ b/examples/features/metadynamics/pimd_metadynamics_zundel/input.xml
@@ -9,13 +9,13 @@
    <ffplumed name="plumed">
          <file mode="xyz">./h5o2+.xyz</file>
          <plumed_dat> plumed/plumed.dat </plumed_dat>
-         <plumed_extras> [doo, co1.lessthan, co2.lessthan, mtd.bias ] </plumed_extras>
+         <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/ipi/engine/forcefields.py b/ipi/engine/forcefields.py
index 491216670..7aeb4f557 100644
--- a/ipi/engine/forcefields.py
+++ b/ipi/engine/forcefields.py
@@ -892,6 +892,16 @@ def mtd_update(self, pos, cell):
         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")
diff --git a/ipi/engine/smotion/metad.py b/ipi/engine/smotion/metad.py
index c48d2782f..b4a36297a 100644
--- a/ipi/engine/smotion/metad.py
+++ b/ipi/engine/smotion/metad.py
@@ -84,10 +84,16 @@ 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!)
 
+                # 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.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.