Fix ENERGY CV in gromacs2024 (#1444)

Author: invemichele-peptone

CHANGES/v2.9.md

@@ -79,3 +79,9 @@ Changes from version 2.8 which are relevant for users:
   - Added FMT keyword in \ref METAINFERENCE
   - Several fixes in regression tests and in CI
   - Fixed an issue where the checkpoint step remains true forever, (see \issue{1428})
+
+## Version 2.9.6 (in progress)
+  - Fixed the GROMACS 2024 patch so that the potential energy is computed on every step
+    where PLUMED needs it (ENERGY CV, OPES multithermal, well-tempered ensemble). See \issue{1205}.
+  - Documentation: \ref ENERGY does include the GROMACS long tail (`DispCorr`) corrections,
+    contrary to what the manual stated. See \issue{567}.

patches/gromacs-2024.3.diff/src/gromacs/mdrun/md.cpp

@@ -1326,6 +1326,14 @@ void gmx::LegacySimulator::do_md()
         force_flags = (GMX_FORCE_STATECHANGED | ((inputrecDynamicBox(ir)) ? GMX_FORCE_DYNAMICBOX : 0)
                        | GMX_FORCE_ALLFORCES | (bCalcVir ? GMX_FORCE_VIRIAL : 0)
                        | (bCalcEner ? GMX_FORCE_ENERGY : 0) | (computeDHDL ? GMX_FORCE_DHDL : 0));
+        /* PLUMED */
+        if (plumedswitch && plumedNeedsEnergy)
+        {
+            /* GROMACS >= 2024 snapshots force_flags into stepWork below.
+               Keep computing energy/virial if PLUMED needed them last step. */
+            force_flags |= GMX_FORCE_ENERGY | GMX_FORCE_VIRIAL;
+        }
+        /* END PLUMED */
         if (simulationWork.useMts && !do_per_step(step, ir->nstfout))
         {
             // TODO: merge this with stepWork.useOnlyMtsCombinedForceBuffer
@@ -1372,7 +1380,7 @@ void gmx::LegacySimulator::do_md()
         {
             // Reset graph on search step (due to changing neighbour list etc)
             // or virial step (due to changing shifts and box).
-            if (bNS || bCalcVir)
+            if (bNS || bCalcVir || plumedNeedsEnergy)
             {
                 fr_->mdGraph[MdGraphEvenOrOddStep::EvenStep]->reset();
                 fr_->mdGraph[MdGraphEvenOrOddStep::OddStep]->reset();
@@ -1381,7 +1389,7 @@ void gmx::LegacySimulator::do_md()
             {
                 mdGraph->setUsedGraphLastStep(usedMdGpuGraphLastStep);
                 bool canUseMdGpuGraphThisStep =
-                        !bNS && !bCalcVir && !doTemperatureScaling && !doParrinelloRahman && !bGStat
+                        !bNS && !bCalcVir && !plumedNeedsEnergy && !doTemperatureScaling && !doParrinelloRahman && !bGStat
                         && !needHalfStepKineticEnergy && !do_per_step(step, ir->nstxout)
                         && !do_per_step(step, ir->nstxout_compressed)
                         && !do_per_step(step, ir->nstvout) && !do_per_step(step, ir->nstfout)
@@ -1471,7 +1479,21 @@ void gmx::LegacySimulator::do_md()
                   if(pversion>3) plumed_cmd(plumedmain,"doCheckPoint",&checkp);
                   plumed_cmd(plumedmain,"setForces",&f.view().force()[0][0]);
                   plumed_cmd(plumedmain,"isEnergyNeeded",&plumedNeedsEnergy);
-                  if(plumedNeedsEnergy) force_flags |= GMX_FORCE_ENERGY | GMX_FORCE_VIRIAL;
+                  if(plumedNeedsEnergy) {
+                    /* GROMACS >= 2024 ignores force_flags at the do_force call:
+                       stepWork was already built above, so rebuild it here. */
+                    force_flags |= GMX_FORCE_ENERGY | GMX_FORCE_VIRIAL;
+                    runScheduleWork_->stepWork = setupStepWorkload(
+                            legacyForceFlags | GMX_FORCE_ENERGY | GMX_FORCE_VIRIAL,
+                            ir->mtsLevels,
+                            step,
+                            runScheduleWork_->domainWork,
+                            simulationWork);
+                    if(!runScheduleWork_->stepWork.computeEnergy)
+                    {
+                      gmx_fatal(FARGS, "PLUMED needs the potential energy but GROMACS is not computing it this step");
+                    }
+                  }
                   clear_mat(plumed_vir);
                   plumed_cmd(plumedmain,"setVirial",&plumed_vir[0][0]);
                 }

src/colvar/Energy.cpp

@@ -38,10 +38,9 @@ it is available, and when also replica exchange is available,
 metadynamics applied to ENERGY can be used to decrease the
 number of required replicas.
 
-\bug This \ref ENERGY does not include long tail corrections.
-Thus when using e.g. LAMMPS `"pair_modify tail yes"` or GROMACS `"DispCorr Ener"` (or `"DispCorr EnerPres"`),
-the potential energy from \ref ENERGY will be slightly different form the one of the MD code.
-You should still be able to use \ref ENERGY and then reweight your simulation with the correct MD energy value.
+\warning Some MD engines communicate to PLUMED the energy without the long tail corrections,
+leading to a small mismatch. If this is the case, use the energy from the MD engine for
+reweighting. See \issue{567}.
 
 \bug Acceptance for replica exchange when \ref ENERGY is biased
 is computed correctly only if all the replicas have the same