GROMACS 2025/2026 patches: follow-up of #1506 (constraint virial with MTS, shells in reruns and with ENERGY, rerun stop, docs)

CHANGES/v2.10.md

@@ -104,7 +104,7 @@ Before switching to version 2.10, users are invited to carefully read the follow
     to what the manual stated. See [issue 567](https://github.com/plumed/plumed2/issues/567).
   - Fixed the GROMACS 2025 and 2026 patches so that PLUMED is told before the force calculation
     that GROMACS writes a checkpoint, as with the GROMACS 2024 patch. Actions that write their
-    state on checkpoints, such as [OPES_METAD](OPES_METAD.md) with STATE_WFILE and no STATE_WSTRIDE,
+    state on checkpoints, such as \ref OPES_METAD with STATE_WFILE and no STATE_WSTRIDE,
     wrote it one step after each GROMACS checkpoint, and never at the last step. With more than one
     MPI rank per simulation, only the main rank wrote it, and the run deadlocked at the first checkpoint.
     See [#1430](https://github.com/plumed/plumed2/pull/1430).
@@ -119,11 +119,11 @@ Before switching to version 2.10, users are invited to carefully read the follow
   - Replica exchange with more than one rank per replica (e.g. GROMACS -replex or -hrex) hung after an exchange
     attempt on which no action needed the atoms: the exchange of positions started for the exchange was never
     completed. With small systems, the following exchanges could instead receive the positions of an earlier step.
-  - GROMACS 2025 and 2026 patches: PLUMED can now stop the simulation, e.g. with COMMITTOR.
-    The request reaches all ranks and the simulations that share state, and mdrun stops at the next
-    step at which GROMACS communicates signals (every nstglobalcomm steps, and at least every 200 steps
-    between simulations that share state; with the modular simulator, at the next neighbor-search step
-    after that), writing a checkpoint there.
+  - GROMACS 2025 and 2026 patches: PLUMED can now stop the simulation, e.g. with \ref COMMITTOR.
+    The request reaches all ranks and the simulations that share state, and mdrun runs one more step
+    after the next step at which GROMACS communicates signals (every nstglobalcomm steps, and between
+    simulations that share state every 200 steps or more; with the modular simulator, it stops at the next
+    neighbor-search step after that), writing a checkpoint at that last step.
   - GROMACS 2025 and 2026 patches: when the PLUMED input needs the potential energy, the forces
     that PLUMED applied on virtual sites were lost. They are now spread on the constructing atoms,
     on the steps when PLUMED uses virtual sites.
@@ -133,6 +133,11 @@ Before switching to version 2.10, users are invited to carefully read the follow
   - GROMACS 2025 and 2026 patches: errors raised on purpose by the PLUMED module, e.g. for replica
     exchange with a bias, are no longer reported as internal errors.
   - GROMACS 2025 and 2026 patches: a request of PLUMED to stop (e.g. \ref COMMITTOR) was ignored in reruns.
+    The rerun now stops after the frame on which PLUMED asks, as with the GROMACS 2020-2024 patches.
   - GROMACS 2025 and 2026 patches: with shell (polarizable) particles, PLUMED also ran on each shell relaxation
     step, as step 1: e.g. \ref PRINT printed extra lines. It now only calculates its forces again on those steps,
-    without updating.
+    without updating, in MD and in reruns. When the PLUMED input needs the energy, PLUMED updates with the energy
+    before the shells are relaxed, and mdrun prints a warning.
+  - GROMACS 2025 and 2026 patches: with multiple time stepping and a PLUMED input that needs the energy, the
+    constraint virial lacked the forces of PLUMED on the slow steps: they are now also added to the force buffer
+    that GROMACS uses for it.

patches/gromacs-2025.0.config

@@ -18,9 +18,10 @@ integration in the next GROMACS version, including:
   the native GROMACS module does not compute or pass the potential energy.
 - checkpoints: PLUMED is told before the force calculation that GROMACS writes a checkpoint,
   so that it writes its own checkpoint data on the same step.
-- stop requests (e.g. COMMITTOR): GROMACS stops at the next step at which it communicates signals,
+- stop requests (e.g. COMMITTOR): GROMACS stops one step after the next step at which it communicates signals,
   writing a checkpoint.
-- polarizable models with shells: PLUMED calculates once per step.
+- polarizable models with shells: PLUMED updates once per step (with ENERGY, before the shells
+  are relaxed: mdrun warns about it).
 
 Patching must be done in the gromacs root directory  _before_ the cmake command is invoked.
 

patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp

@@ -191,7 +191,7 @@ try
     long int         lstep = prepareStepIsCalled_ ? currentStep_ : forceProviderInput.step_;
     // With shell particles, relax_shell_flexcon() calls do_force() again on the same step (passing
     // step 1): PLUMED then calculates its forces again, without updating (depositing, printing...).
-    // Only in MD, where prepareStep() gives the step: e.g. reruns can have frames with the same step
+    // "Again" means since the last prepareStep(), which MD and reruns call before each step
     repeatedCall_   = prepareStepIsCalled_ && (calculatedStep_ == lstep);
     calculatedStep_ = lstep;
     plumed_->cmd("setStepLong", &lstep);
@@ -298,6 +298,7 @@ try
     plumed_->cmd("setStepLong", &lstep);
     prepareStepIsCalled_ = true;
     currentStep_         = step;
+    calculatedStep_.reset();
     if (plumedAPIversion_ > 3)
     {
         // Set on every step: PLUMED resets it only on the steps when some action is active

patches/gromacs-2025.0.diff/src/gromacs/mdlib/sim_util.cpp

@@ -2637,10 +2637,17 @@ void do_force(FILE*                         fplog,
     /* PLUMED */
     if (fr->forceProviders != nullptr)
     {
-        ArrayRef<RVec> forceForPlumed = force.unpaddedArrayRef();
+        ArrayRef<RVec>    forceForPlumed = force.unpaddedArrayRef();
+        std::vector<RVec> forceBeforePlumed;
         if (simulationWork.useMts && stepWork.computeSlowForces)
         {
+            /* The update uses the combined MTS buffer, and the constraint virial the normal one,
+               which is complete when the virial is computed: PLUMED's forces are added to it below */
             forceForPlumed = forceView->forceMtsCombined();
+            if (stepWork.computeVirial)
+            {
+                forceBeforePlumed.assign(forceForPlumed.begin(), forceForPlumed.begin() + mdatoms->homenr);
+            }
         }
         const bool plumedForcesOnVirtualSites = fr->forceProviders->applyAfterPotentialEnergy(
                 stepWork.computeEnergy, &enerd->term[F_EPOT], forceForPlumed, vir_force);
@@ -2662,6 +2669,14 @@ void do_force(FILE*                         fplog,
                                 wcycle);
             m_add(vir_force, virial, vir_force);
         }
+        if (!forceBeforePlumed.empty())
+        {
+            ArrayRef<RVec> forceNormal = force.unpaddedArrayRef();
+            for (int i = 0; i < mdatoms->homenr; i++)
+            {
+                forceNormal[i] += forceForPlumed[i] - forceBeforePlumed[i];
+            }
+        }
     }
     /* END PLUMED */
 

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

@@ -870,6 +870,9 @@ void gmx::LegacySimulator::do_md()
     }
 
     bool usedMdGpuGraphLastStep = false;
+    /* PLUMED */
+    bool plumedWarnedAboutShells = false;
+    /* END PLUMED */
     /* and stop now if we should */
     bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
     while (!bLastStep)
@@ -1155,6 +1158,19 @@ void gmx::LegacySimulator::do_md()
                           step,
                           ir->mtsLevels.back().stepFactor);
             }
+            if (shellfc != nullptr && !plumedWarnedAboutShells)
+            {
+                GMX_LOG(mdLog_.warning)
+                        .asParagraph()
+                        .appendText(
+                                "With shell particles or flexible constraints, PLUMED updates (e.g. "
+                                "prints, deposits its bias, checks COMMITTOR) on the first force "
+                                "calculation of each step, before the shells are relaxed. This PLUMED "
+                                "input needs the potential energy: the energy that it uses for these "
+                                "updates, like the collective variables that depend on the shells, is "
+                                "the one before the relaxation.");
+                plumedWarnedAboutShells = true;
+            }
             force_flags |= GMX_FORCE_ENERGY | GMX_FORCE_VIRIAL;
         }
         /* END PLUMED */

patches/gromacs-2025.0.diff/src/gromacs/mdrun/rerun.cpp

@@ -111,6 +111,7 @@
 #include "gromacs/mdtypes/forcebuffers.h"
 #include "gromacs/mdtypes/forcerec.h"
 #include "gromacs/mdtypes/group.h"
+#include "gromacs/mdtypes/iforceprovider.h" /* PLUMED */
 #include "gromacs/mdtypes/inputrec.h"
 #include "gromacs/mdtypes/interaction_const.h"
 #include "gromacs/mdtypes/md_enums.h"
@@ -645,6 +646,15 @@ void gmx::LegacySimulator::do_rerun()
                                                        runScheduleWork_->domainWork,
                                                        runScheduleWork_->simulationWork);
 
+        /* PLUMED */
+        /* As in md.cpp, PLUMED gets the step before the force calculation: with shell particles, it
+           then only calculates again, without updating, in the rest of the shell relaxation */
+        if (fr_->forceProviders != nullptr)
+        {
+            fr_->forceProviders->prepareStep(step, false);
+        }
+        /* END PLUMED */
+
         if (shellfc)
         {
             /* Now is the time to relax the shells */
@@ -892,8 +902,9 @@ void gmx::LegacySimulator::do_rerun()
         if (MAIN(cr_))
         {
             /* read next frame from input trajectory */
-            /* PLUMED: keep a stop request, e.g. from PLUMED (COMMITTOR) through the stop handler */
-            isLastStep = isLastStep || !read_next_frame(oenv_, status, &rerun_fr);
+            /* PLUMED: stop after this frame if asked, e.g. by PLUMED (COMMITTOR) through the stop handler */
+            isLastStep = isLastStep || stopHandler->stoppingAfterCurrentStep(bNS)
+                         || !read_next_frame(oenv_, status, &rerun_fr);
         }
 
         if (PAR(cr_))

patches/gromacs-2026.0.config

@@ -17,9 +17,10 @@ This patch provides extra features to be implemented in the official PLUMED inte
   the native GROMACS module does not compute or pass the potential energy
 - checkpoints: PLUMED is told before the force calculation that GROMACS writes a checkpoint,
   so that it writes its own checkpoint data on the same step
-- stop requests (e.g. COMMITTOR): GROMACS stops at the next step at which it communicates signals,
+- stop requests (e.g. COMMITTOR): GROMACS stops one step after the next step at which it communicates signals,
   writing a checkpoint
-- polarizable models with shells: PLUMED calculates once per step
+- polarizable models with shells: PLUMED updates once per step (with ENERGY, before the shells
+  are relaxed: mdrun warns about it)
 
 Patching must be done in the gromacs root directory  _before_ the cmake command is invoked.
 

patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp

@@ -209,7 +209,7 @@ try
     long int lstep = prepareStepIsCalled_ ? currentStep_ : forceProviderInput.step_;
     // With shell particles, relax_shell_flexcon() calls do_force() again on the same step (passing
     // step 1): PLUMED then calculates its forces again, without updating (depositing, printing...).
-    // Only in MD, where prepareStep() gives the step: e.g. reruns can have frames with the same step
+    // "Again" means since the last prepareStep(), which MD and reruns call before each step
     repeatedCall_   = prepareStepIsCalled_ && (calculatedStep_ == lstep);
     calculatedStep_ = lstep;
     plumed_->cmd("setStepLong", &lstep);
@@ -342,6 +342,7 @@ try
     plumed_->cmd("setStepLong", &lstep);
     prepareStepIsCalled_ = true;
     currentStep_         = step;
+    calculatedStep_.reset();
     if (plumedAPIversion_ > 3)
     {
         // Set on every step: PLUMED resets it only on the steps when some action is active

patches/gromacs-2026.0.diff/src/gromacs/mdlib/sim_util.cpp

@@ -2735,10 +2735,17 @@ void do_force(FILE*                         fplog,
     /* PLUMED */
     if (fr->forceProviders != nullptr)
     {
-        ArrayRef<RVec> forceForPlumed = force.unpaddedArrayRef();
+        ArrayRef<RVec>    forceForPlumed = force.unpaddedArrayRef();
+        std::vector<RVec> forceBeforePlumed;
         if (simulationWork.useMts && stepWork.computeSlowForces)
         {
+            /* The update uses the combined MTS buffer, and the constraint virial the normal one,
+               which is complete when the virial is computed: PLUMED's forces are added to it below */
             forceForPlumed = forceView->forceMtsCombined();
+            if (stepWork.computeVirial)
+            {
+                forceBeforePlumed.assign(forceForPlumed.begin(), forceForPlumed.begin() + mdatoms->homenr);
+            }
         }
         const bool plumedForcesOnVirtualSites = fr->forceProviders->applyAfterPotentialEnergy(
                 stepWork.computeEnergy, &enerd->term[InteractionFunction::PotentialEnergy], forceForPlumed, vir_force);
@@ -2760,6 +2767,14 @@ void do_force(FILE*                         fplog,
                                 wcycle);
             m_add(vir_force, virial, vir_force);
         }
+        if (!forceBeforePlumed.empty())
+        {
+            ArrayRef<RVec> forceNormal = force.unpaddedArrayRef();
+            for (int i = 0; i < mdatoms->homenr; i++)
+            {
+                forceNormal[i] += forceForPlumed[i] - forceBeforePlumed[i];
+            }
+        }
     }
     /* END PLUMED */
 

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

@@ -863,6 +863,9 @@ void gmx::LegacySimulator::do_md()
     }
 
     bool usedMdGpuGraphLastStep = false;
+    /* PLUMED */
+    bool plumedWarnedAboutShells = false;
+    /* END PLUMED */
     /* and stop now if we should */
     bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
     while (!bLastStep)
@@ -1141,6 +1144,19 @@ void gmx::LegacySimulator::do_md()
                           step,
                           ir->mtsLevels.back().stepFactor);
             }
+            if (shellfc != nullptr && !plumedWarnedAboutShells)
+            {
+                GMX_LOG(mdLog_.warning)
+                        .asParagraph()
+                        .appendText(
+                                "With shell particles or flexible constraints, PLUMED updates (e.g. "
+                                "prints, deposits its bias, checks COMMITTOR) on the first force "
+                                "calculation of each step, before the shells are relaxed. This PLUMED "
+                                "input needs the potential energy: the energy that it uses for these "
+                                "updates, like the collective variables that depend on the shells, is "
+                                "the one before the relaxation.");
+                plumedWarnedAboutShells = true;
+            }
             force_flags |= GMX_FORCE_ENERGY | GMX_FORCE_VIRIAL;
         }
         /* END PLUMED */

patches/gromacs-2026.0.diff/src/gromacs/mdrun/rerun.cpp

@@ -109,6 +109,7 @@
 #include "gromacs/mdtypes/forcebuffers.h"
 #include "gromacs/mdtypes/forcerec.h"
 #include "gromacs/mdtypes/group.h"
+#include "gromacs/mdtypes/iforceprovider.h" /* PLUMED */
 #include "gromacs/mdtypes/inputrec.h"
 #include "gromacs/mdtypes/interaction_const.h"
 #include "gromacs/mdtypes/md_enums.h"
@@ -646,6 +647,15 @@ void gmx::LegacySimulator::do_rerun()
                                                        runScheduleWork_->domainWork,
                                                        runScheduleWork_->simulationWork);
 
+        /* PLUMED */
+        /* As in md.cpp, PLUMED gets the step before the force calculation: with shell particles, it
+           then only calculates again, without updating, in the rest of the shell relaxation */
+        if (fr_->forceProviders != nullptr)
+        {
+            fr_->forceProviders->prepareStep(step, false);
+        }
+        /* END PLUMED */
+
         if (shellfc)
         {
             /* Now is the time to relax the shells */
@@ -891,8 +901,9 @@ void gmx::LegacySimulator::do_rerun()
         if (isMainRank)
         {
             /* read next frame from input trajectory */
-            /* PLUMED: keep a stop request, e.g. from PLUMED (COMMITTOR) through the stop handler */
-            isLastStep = isLastStep || !read_next_frame(oenv_, status, &rerun_fr);
+            /* PLUMED: stop after this frame if asked, e.g. by PLUMED (COMMITTOR) through the stop handler */
+            isLastStep = isLastStep || stopHandler->stoppingAfterCurrentStep(step)
+                         || !read_next_frame(oenv_, status, &rerun_fr);
         }
 
         if (cr_->commMyGroup.isParallel())

user-doc/Introduction.md

@@ -60,15 +60,22 @@ and the patches for the older versions of GROMACS have some more:
 
 | Feature | GROMACS 2025 or later, without patches | patches for GROMACS 2025 and 2026 | patches for GROMACS 2022 to 2024 |
 |:--|:--:|:--:|:--:|
-| biases on atoms, with domain decomposition and multiple time stepping | yes | yes | yes |
-| \ref ENERGY collective variable and biases on the energy | no | yes | yes |
-| multiple walkers (e.g. `WALKERS_MPI`) | no | yes | yes |
-| PLUMED checkpoint data written at the steps of the GROMACS checkpoints | no | yes | yes |
+| biases on atoms, with domain decomposition and multiple time stepping | yes (1) | yes | yes |
+| \ref ENERGY collective variable and biases on the energy | no | yes (2,3) | yes (3) |
+| multiple simulations that communicate through MPI (e.g. `WALKERS_MPI`) | no | yes | yes |
+| PLUMED checkpoint data written at the steps of the GROMACS checkpoints | no | yes (2) | yes |
 | PLUMED stopping the simulation (e.g. \ref COMMITTOR) | no | yes | yes |
-| polarizable models with shells | no | yes | no |
+| polarizable models with shells | no | yes (4) | no |
 | replica exchange with a PLUMED bias (`-replex`, `-hrex`) | no | no | yes |
 | lambda dynamics (expanded ensemble) | no | no | yes |
 
+Notes on the table:
+
+1. GROMACS 2026.0 to 2026.3 without patches stop at startup with more than one MPI rank.
+2. Not with the modular simulator of GROMACS, the default for md-vv: set `GMX_DISABLE_MODULAR_SIMULATOR=1`.
+3. With multiple time stepping, only on the steps when the slow forces are computed.
+4. With \ref ENERGY, PLUMED updates (e.g. deposits its bias) with the energy before the shells are relaxed.
+
 The interface with GROMACS is tested on the [PLUMED testcenter](https://plumed-testcenter.github.io/browse.html).
 
 Alternatively, one