GROMACS 2025/2026 patches and kernel: fix the virial with an ENERGY bias, a replica exchange deadlock, COMMITTOR, virtual sites, MTS and shells

CHANGES/v2.10.md

@@ -113,3 +113,26 @@ Before switching to version 2.10, users are invited to carefully read the follow
   - Fixed the GROMACS 2025 patch so that `-multidir` simulations without `-plumed` do not have to run
     in lockstep, which made them deadlock when they had different numbers of steps. As in the
     GROMACS 2026 patch, the multi-simulation record is now used only when PLUMED is active.
+  - Fixed the virial with a bias on ENERGY (e.g. OPES multithermal or the well-tempered ensemble) and
+    domain decomposition: only the first rank rescaled its share of the virial of the MD code, so the
+    pressure was wrong, by orders of magnitude in some cases. As in PLUMED 2.9, it is now rescaled on all ranks.
+  - 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: 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.
+  - GROMACS 2025 and 2026 patches: with multiple time stepping, the potential energy passed to PLUMED
+    on the fast steps lacked the slow terms (e.g. the PME mesh, by 650-1250 kJ/mol in a water test).
+    mdrun now stops with an error if PLUMED needs the energy on such a step.
+  - 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.
+  - 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.

patches/gromacs-2025.0.config

@@ -18,10 +18,13 @@ 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,
+  writing a checkpoint.
+- polarizable models with shells: PLUMED calculates once per step.
 
 Patching must be done in the gromacs root directory  _before_ the cmake command is invoked.
 
-The GROMACS integration has some improvements with the compatibility wiht the multi-threading implementation
+The GROMACS integration has some improvements with the compatibility with the multi-threading implementation
 but it is still preferable to configure gromacs as
 
 cmake -DGMX_THREAD_MPI=OFF and add -DGMX_MPI=ON if you want to use MPI.

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

@@ -55,7 +55,7 @@
 #include "plumedOptions.h"
 #include "plumedforceprovider.h"
 
-class gmx_multisim_t;
+struct gmx_multisim_t;
 
 namespace gmx
 {

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

@@ -44,6 +44,9 @@
 #include "gromacs/math/vec.h"
 #include "gromacs/mdrunutility/handlerestart.h"
 #include "gromacs/mdrunutility/mdmodulesnotifiers.h"
+#include "gromacs/topology/atoms.h"
+#include "gromacs/topology/mtop_atomloops.h"
+#include "gromacs/topology/mtop_util.h"
 #include "gromacs/topology/topology.h"
 
 namespace gmx
@@ -52,6 +55,16 @@ namespace gmx
 void PlumedOptionProvider::setTopology(const gmx_mtop_t& mtop)
 {
     opts_.natoms_ = mtop.natoms;
+    opts_.isVirtualSite_.clear();
+    if (gmx_mtop_particletype_count(mtop)[ParticleType::VSite] > 0)
+    {
+        opts_.isVirtualSite_.resize(mtop.natoms, false);
+        for (const AtomProxy atomP : AtomRange(mtop))
+        {
+            opts_.isVirtualSite_[atomP.globalAtomNumber()] =
+                    (atomP.atom().ptype == ParticleType::VSite);
+        }
+    }
 }
 
 void PlumedOptionProvider::setEnsembleTemperature(const EnsembleTemperature& temp)

patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h

@@ -44,6 +44,7 @@
 
 #include <optional>
 #include <string>
+#include <vector>
 
 #include "gromacs/utility/real.h"
 
@@ -67,6 +68,8 @@ struct PlumedOptions
     StartingBehavior    startingBehavior_{};
     bool                active_{ false };
     bool                replex_{ false };
+    //! Whether each atom is a virtual site, empty if there are none
+    std::vector<bool>   isVirtualSite_;
 };
 
 class PlumedOptionProvider

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

@@ -41,6 +41,8 @@
 
 #include "plumedforceprovider.h"
 
+#include <algorithm>
+
 #include "gromacs/domdec/domdec.h"
 #include "gromacs/domdec/domdec_struct.h"
 #include "gromacs/gmxlib/network.h"
@@ -74,7 +76,9 @@ namespace gmx
 
 PlumedForceProvider::~PlumedForceProvider() = default;
 PlumedForceProvider::PlumedForceProvider(const PlumedOptions& options)
-try : plumed_(std::make_unique<PLMD::Plumed>()),replex_(options.replex_)
+try : plumed_(std::make_unique<PLMD::Plumed>()),
+      replex_(options.replex_),
+      isVirtualSite_(options.isVirtualSite_)
 {
     // I prefer to pass a struct with data because it stops the coupling
     // at the implementation and not at the function signature:
@@ -147,6 +151,11 @@ try : plumed_(std::make_unique<PLMD::Plumed>()),replex_(options.replex_)
     plumed_->cmd("setTimestep", &options.simulationTimeStep_);
     plumed_->cmd("init", nullptr);
 }
+catch (const GromacsException&)
+{
+    // Keep the type and message of the errors raised on purpose, e.g. NotImplementedError
+    throw;
+}
 catch (const std::exception& ex)
 {
     GMX_THROW(InternalError(
@@ -161,6 +170,11 @@ try
         checkpointRequested_ = 1;
     }
 }
+catch (const GromacsException&)
+{
+    // Keep the type and message of the errors raised on purpose, e.g. NotImplementedError
+    throw;
+}
 catch (const std::exception& ex)
 {
     GMX_THROW(InternalError(
@@ -174,7 +188,12 @@ try
     // setup: these instructions in the original patch are BEFORE do_force()
     // now this is called within do_force(), but this does not impact the results
     const t_commrec* cr    = &(forceProviderInput.cr_);
-    long int         lstep = forceProviderInput.step_;
+    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
+    repeatedCall_   = prepareStepIsCalled_ && (calculatedStep_ == lstep);
+    calculatedStep_ = lstep;
     plumed_->cmd("setStepLong", &lstep);
     if (haveDDAtomOrdering(*cr))
     {
@@ -201,7 +220,10 @@ try
     }
     preparedStep_.reset();
 
-    plumedWantsToStop_ = 0;
+    if (!repeatedCall_)
+    {
+        plumedWantsToStop_ = 0;
+    }
     plumed_->cmd("setStopFlag", &plumedWantsToStop_);
 
     if (!prepareStepIsCalled_)
@@ -229,19 +251,42 @@ try
     clear_mat(plumedVir_);
     plumed_->cmd("setVirial", &plumedVir_[0][0]);
 
-    plumed_->cmd("performCalc", nullptr);
+    plumed_->cmd(repeatedCall_ ? "performCalcNoUpdate" : "performCalc", nullptr);
 
     checkReplicaExchangeBias();
 
     msmul(plumedVir_, 0.5, plumedVir_);
     forceProviderOutput->forceWithVirial_.addVirialContribution(plumedVir_);
 }
+catch (const GromacsException&)
+{
+    // Keep the type and message of the errors raised on purpose, e.g. NotImplementedError
+    throw;
+}
 catch (const std::exception& ex)
 {
     GMX_THROW(InternalError(
             std::string("An error occurred while PLUMED was calculating the forces\n:") + ex.what()));
 }
 
+bool PlumedForceProvider::usesVirtualSites()
+{
+    if (isVirtualSite_.empty())
+    {
+        return false;
+    }
+    /* PLUMED applies forces only on the atoms it uses. Their list is the same on all ranks,
+       so that all take the same decision. It is empty when PLUMED only biases the energy. */
+    int numAtoms = 0;
+    plumed_->cmd("createFullList", &numAtoms);
+    const int* atoms = nullptr;
+    plumed_->cmd("getFullList", &atoms);
+    const bool found =
+            std::any_of(atoms, atoms + numAtoms, [this](int atom) { return isVirtualSite_[atom]; });
+    plumed_->cmd("clearFullList", nullptr);
+    return found;
+}
+
 bool PlumedForceProvider::prepareStep(int64_t step, bool isCheckpointingStep)
 try
 {
@@ -252,6 +297,7 @@ try
     long int lstep = step;
     plumed_->cmd("setStepLong", &lstep);
     prepareStepIsCalled_ = true;
+    currentStep_         = step;
     if (plumedAPIversion_ > 3)
     {
         // Set on every step: PLUMED resets it only on the steps when some action is active
@@ -264,12 +310,22 @@ try
     preparedStep_ = step;
     return plumedNeedsEnergy_ != 0;
 }
+catch (const GromacsException&)
+{
+    // Keep the type and message of the errors raised on purpose, e.g. NotImplementedError
+    throw;
+}
 catch (const std::exception& ex)
 {
     GMX_THROW(InternalError(
             std::string("An error occurred while PLUMED was preparing the step\n:") + ex.what()));
 }
 
+bool PlumedForceProvider::requestsStop() const
+{
+    return plumedWantsToStop_ != 0;
+}
+
 void PlumedForceProvider::checkReplicaExchangeBias()
 {
     if (replex_)
@@ -284,15 +340,15 @@ void PlumedForceProvider::checkReplicaExchangeBias()
     }
 }
 
-void PlumedForceProvider::applyAfterPotentialEnergy(bool          energyWasComputed,
+bool PlumedForceProvider::applyAfterPotentialEnergy(bool          energyWasComputed,
                                                     real*         potentialEnergy,
                                                     ArrayRef<RVec> force,
                                                     tensor        virial)
 try
 {
     if (!plumedNeedsEnergy_)
     {
-        return;
+        return false;
     }
     if (!energyWasComputed || potentialEnergy == nullptr)
     {
@@ -312,11 +368,17 @@ try
     // Match the legacy GROMACS patch: pass 2*virial into PLUMED, then replace.
     msmul(virial, 2.0, plumedVir_);
     plumed_->cmd("setVirial", &plumedVir_[0][0]);
-    plumed_->cmd("performCalc", nullptr);
+    plumed_->cmd(repeatedCall_ ? "performCalcNoUpdate" : "performCalc", nullptr);
 
     checkReplicaExchangeBias();
 
     msmul(plumedVir_, 0.5, virial);
+    return usesVirtualSites();
+}
+catch (const GromacsException&)
+{
+    // Keep the type and message of the errors raised on purpose, e.g. NotImplementedError
+    throw;
 }
 catch (const std::exception& ex)
 {

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

@@ -43,6 +43,7 @@
 
 #include <memory>
 #include <optional>
+#include <vector>
 
 #include "gromacs/mdtypes/iforceprovider.h"
 
@@ -81,17 +82,22 @@ class PlumedForceProvider final : public IForceProvider
     void calculateForces(const ForceProviderInput& forceProviderInput,
                          ForceProviderOutput*      forceProviderOutput) override;
     bool prepareStep(int64_t step, bool isCheckpointingStep) override;
-    void applyAfterPotentialEnergy(bool energyWasComputed,
+    bool requestsStop() const override;
+    bool applyAfterPotentialEnergy(bool energyWasComputed,
                                    real* potentialEnergy,
                                    ArrayRef<RVec> force,
                                    tensor virial) override;
 
 private:
     void checkReplicaExchangeBias();
+    //! Whether the atoms that PLUMED used on this step include virtual sites
+    bool usesVirtualSites();
 
     std::unique_ptr<PLMD::Plumed> plumed_;
     int                           plumedAPIversion_;
     bool                          replex_;
+    //! Whether each atom is a virtual site, empty if there are none
+    std::vector<bool>             isVirtualSite_;
     int                           plumedNeedsEnergy_ = 0;
     int                           plumedWantsToStop_ = 0;
     matrix                        plumedVir_         = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } };
@@ -101,6 +107,12 @@ class PlumedForceProvider final : public IForceProvider
     bool prepareStepIsCalled_ = false;
     //! Set on the main rank by writeCheckpointData(), passed to PLUMED by calculateForces().
     int checkpointRequested_ = 0;
+    //! Step for which prepareStep() was last called
+    int64_t currentStep_ = 0;
+    //! Last step that PLUMED calculated
+    std::optional<int64_t> calculatedStep_;
+    //! Whether PLUMED already calculated this MD step (shell relaxation): it is then not updated again
+    bool repeatedCall_ = false;
 };
 
 } // namespace gmx

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

@@ -2642,8 +2642,26 @@ void do_force(FILE*                         fplog,
         {
             forceForPlumed = forceView->forceMtsCombined();
         }
-        fr->forceProviders->applyAfterPotentialEnergy(
+        const bool plumedForcesOnVirtualSites = fr->forceProviders->applyAfterPotentialEnergy(
                 stepWork.computeEnergy, &enerd->term[F_EPOT], forceForPlumed, vir_force);
+        if (plumedForcesOnVirtualSites && vsite)
+        {
+            /* The forces on the virtual sites were spread above and are now zero, except
+               those PLUMED just applied: spread these too, as in postProcessForces() */
+            const VirtualSitesHandler::VirialHandling virialHandling =
+                    (stepWork.computeVirial ? VirtualSitesHandler::VirialHandling::NonLinear
+                                            : VirtualSitesHandler::VirialHandling::None);
+            matrix virial = { { 0 } };
+            vsite->spreadForces(x.unpaddedArrayRef(),
+                                forceForPlumed,
+                                virialHandling,
+                                {},
+                                virial,
+                                nrnb,
+                                box,
+                                wcycle);
+            m_add(vir_force, virial, vir_force);
+        }
     }
     /* END PLUMED */
 

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

@@ -1143,6 +1143,18 @@ void gmx::LegacySimulator::do_md()
         {
             /* setupStepWorkload() snapshots force_flags below, so PLUMED has to be
                asked for the potential energy before that, not from within do_force(). */
+            if (simulationWork.useMts && !do_per_step(step, ir->mtsLevels.back().stepFactor))
+            {
+                gmx_fatal(FARGS,
+                          "PLUMED needs the potential energy at step %" PRId64
+                          ", but with multiple time stepping the slow forces, and so the "
+                          "potential energy, are computed only every %d steps "
+                          "(mts-level2-factor). Use a PLUMED input that needs the energy only "
+                          "on those steps (e.g. no bias on ENERGY), or do not use multiple time "
+                          "stepping.",
+                          step,
+                          ir->mtsLevels.back().stepFactor);
+            }
             force_flags |= GMX_FORCE_ENERGY | GMX_FORCE_VIRIAL;
         }
         /* END PLUMED */