From 6a39b78d063b6cc51a6a7a7b37101a047e2d30b4 Mon Sep 17 00:00:00 2001 From: Michele Invernizzi Date: Wed, 14 Oct 2026 16:44:21 +0200 Subject: [PATCH 01/10] Rescale the virial on all domains when ENERGY is biased On the domains other than the first, the Box input was flagged noforce, so that the virial computed by PLUMED is added only once. But the same flag also stopped ENERGY from rescaling that rank's share of the virial of the MD code, which is spread over all ranks. The pressure was then wrong with any bias on ENERGY and domain decomposition: in a GROMACS test, the virial at step 0 was -725.7 on 2 ranks instead of 5.457. PLUMED 2.9 rescaled the virial on all ranks. A new flag, forces_added_elsewhere, now only stops adding PLUMED's forces on those ranks. rescale_force skips a buffer the MD code did not pass, as PLUMED 2.9 did for the virial. Co-Authored-By: Claude Opus 5.5 --- CHANGES/v2.10.md | 3 +++ src/core/ActionToPutData.cpp | 3 ++- src/core/ActionToPutData.h | 3 +++ src/core/DataPassingObject.cpp | 4 ++++ src/core/DomainDecomposition.cpp | 4 +++- 5 files changed, 15 insertions(+), 2 deletions(-) diff --git a/CHANGES/v2.10.md b/CHANGES/v2.10.md index 78d6a4d29f..bfc0cecb7a 100644 --- a/CHANGES/v2.10.md +++ b/CHANGES/v2.10.md @@ -113,3 +113,6 @@ 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. diff --git a/src/core/ActionToPutData.cpp b/src/core/ActionToPutData.cpp index 3f9e1a00ff..6741929fa5 100644 --- a/src/core/ActionToPutData.cpp +++ b/src/core/ActionToPutData.cpp @@ -57,6 +57,7 @@ ActionToPutData::ActionToPutData(const ActionOptions&ao): Action(ao), ActionForInterface(ao), noforce(false), + forces_added_elsewhere(false), fixed(false), from_domains(false), resetable(false), @@ -226,7 +227,7 @@ void ActionToPutData::wait() { } void ActionToPutData::apply() { - if( getPntrToValue()->forcesWereAdded() && !noforce ) { + if( getPntrToValue()->forcesWereAdded() && !noforce && !forces_added_elsewhere ) { if( getName()=="ENERGY" || getDependencies().size()==0 ) { mydata->add_force( getPntrToValue() ); } diff --git a/src/core/ActionToPutData.h b/src/core/ActionToPutData.h index 6efe24cee4..53b98b8e56 100644 --- a/src/core/ActionToPutData.h +++ b/src/core/ActionToPutData.h @@ -37,6 +37,9 @@ class ActionToPutData : private: /// Are we not applying forces on this values bool noforce; +/// Are the forces on this value added by another rank? The forces that the MD code +/// holds in this value are still rescaled on this rank (e.g. its share of the virial) + bool forces_added_elsewhere; /// Is this quantity fixed bool fixed; /// Is this quantity passed from the domains diff --git a/src/core/DataPassingObject.cpp b/src/core/DataPassingObject.cpp index 6f44a68cc3..4a7844b918 100644 --- a/src/core/DataPassingObject.cpp +++ b/src/core/DataPassingObject.cpp @@ -268,6 +268,10 @@ void DataPassingObjectTyped::rescale_force( const unsigned& n, const double& } T* pp; getPointer( f, s, start, stride, pp ); + // The MD code might not pass this buffer on every rank (e.g. the virial) + if( !pp ) { + return; + } #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(pp,n)) for(unsigned i=0; i("Box"); - ap->noforce=true; + ap->forces_added_elsewhere=true; } } From c873aeebe261defc2935b512e9afae06ad249b52 Mon Sep 17 00:00:00 2001 From: Michele Invernizzi Date: Wed, 14 Oct 2026 16:44:21 +0200 Subject: [PATCH 02/10] GROMACS 2025/2026 patches: let PLUMED stop the simulation The patches passed setStopFlag to PLUMED but never read the flag, so COMMITTOR could not stop mdrun. The GROMACS 2024 patch set bLastStep after the force calculation, which with simulations that share state (e.g. WALKERS_MPI) stops one simulation and deadlocks the others. The force provider now reports the request with requestsStop(), and runner.cpp registers it as a stop condition with the StopHandlerBuilder. GROMACS then passes it to all ranks and to the simulations that share state, in the legacy and in the modular simulator, and mdrun stops at the next step at which GROMACS communicates signals (every nstglobalcomm steps, and at least every 200 steps between simulations), writing a checkpoint, with exit code 0. Only the main rank signals, since the flag is the same on all ranks. The modular simulator stops at the next neighbor-search step instead: when a stop is signalled for another step, it clears the task queue it is running from, and crashes (seen with md-vv). Tested on ala2 with COMMITTOR: the run stops 4-9 steps after committing with leap-frog on 1 and 2 ranks, 13-25 steps after with md-vv, and two walkers stop together. Before, all ran to the end. Co-Authored-By: Claude Opus 5.5 --- CHANGES/v2.10.md | 5 +++ .../plumed/plumedforceprovider.cpp | 5 +++ .../plumed/plumedforceprovider.h | 1 + .../src/gromacs/mdrun/runner.cpp | 32 +++++++++++++++++++ .../src/gromacs/mdtypes/iforceprovider.cpp | 12 +++++++ .../src/gromacs/mdtypes/iforceprovider.h | 11 +++++++ .../plumed/plumedforceprovider.cpp | 5 +++ .../plumed/plumedforceprovider.h | 2 ++ .../src/gromacs/mdrun/runner.cpp | 32 +++++++++++++++++++ .../src/gromacs/mdtypes/iforceprovider.cpp | 12 +++++++ .../src/gromacs/mdtypes/iforceprovider.h | 11 +++++++ 11 files changed, 128 insertions(+) diff --git a/CHANGES/v2.10.md b/CHANGES/v2.10.md index bfc0cecb7a..d3994244e9 100644 --- a/CHANGES/v2.10.md +++ b/CHANGES/v2.10.md @@ -116,3 +116,8 @@ Before switching to version 2.10, users are invited to carefully read the follow - 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. + - 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. diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp index 34b3ba7c22..db075fb0f8 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp @@ -270,6 +270,11 @@ catch (const std::exception& ex) 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_) diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h index 69dcc08a56..9741993169 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h @@ -81,6 +81,7 @@ class PlumedForceProvider final : public IForceProvider void calculateForces(const ForceProviderInput& forceProviderInput, ForceProviderOutput* forceProviderOutput) override; bool prepareStep(int64_t step, bool isCheckpointingStep) override; + bool requestsStop() const override; void applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef force, diff --git a/patches/gromacs-2025.0.diff/src/gromacs/mdrun/runner.cpp b/patches/gromacs-2025.0.diff/src/gromacs/mdrun/runner.cpp index 70f5bc976c..5f6c41313a 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/mdrun/runner.cpp +++ b/patches/gromacs-2025.0.diff/src/gromacs/mdrun/runner.cpp @@ -138,6 +138,7 @@ #include "gromacs/mdtypes/fcdata.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/group.h" +#include "gromacs/mdtypes/iforceprovider.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/interaction_const.h" #include "gromacs/mdtypes/locality.h" @@ -2313,6 +2314,37 @@ int Mdrunner::mdrunner() simulatorBuilder.add(SimulatorStateData( globalState.get(), localState, &observablesHistory, &enerd, &ekind)); simulatorBuilder.add(std::move(membedHolder)); + /* PLUMED */ + if (fr->forceProviders != nullptr && isSimulationMainRank) + { + /* PLUMED can ask to stop the simulation, e.g. COMMITTOR. As a stop condition, the + request reaches all ranks and the simulations that share state, in both the + legacy and the modular simulator. The simulation stops after the next step + with global communication. The request is the same on all ranks, so only + the main rank signals it, which also keeps the summed signal small. + The modular simulator stops at the next neighbor-search step: when asked to + stop at another step, it clears the task queue it is running, and crashes. */ + const StopSignal plumedStopSignal = useModularSimulator ? StopSignal::stopAtNextNSStep + : StopSignal::stopImmediately; + stopHandlerBuilder_->registerStopCondition( + [forceProviders = fr->forceProviders, + plumedStopSignal, + fplog, + reported = std::make_shared(false)]() + { + if (!forceProviders->requestsStop()) + { + return StopSignal::noSignal; + } + if (fplog != nullptr && !*reported) + { + fprintf(fplog, "\n\nPLUMED asked to stop the simulation\n\n"); + *reported = true; + } + return plumedStopSignal; + }); + } + /* END PLUMED */ simulatorBuilder.add(std::move(stopHandlerBuilder_)); simulatorBuilder.add(SimulatorConfig(mdrunOptions, startingBehavior, &runScheduleWork)); diff --git a/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.cpp b/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.cpp index 26fe70f5fd..f7a158a438 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.cpp +++ b/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.cpp @@ -109,6 +109,18 @@ bool ForceProviders::prepareStep(int64_t step, bool isCheckpointingStep) const return needed; } +bool ForceProviders::requestsStop() const +{ + for (const auto& provider : impl_->providers_) + { + if (provider.first->requestsStop()) + { + return true; + } + } + return false; +} + void ForceProviders::applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef force, diff --git a/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.h b/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.h index 56916fd4f3..68c32e5423 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.h +++ b/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.h @@ -193,6 +193,14 @@ class IForceProvider */ virtual bool prepareStep(int64_t /*step*/, bool /*isCheckpointingStep*/) { return false; } + /*! \brief + * Whether this provider asks to stop the simulation, after the forces it computed last. + * + * Used by PLUMED, e.g. when COMMITTOR finds the system in a basin. The MD loop gets the + * request through a stop condition registered with the StopHandlerBuilder. + */ + virtual bool requestsStop() const { return false; } + /*! \brief * Optional callback after the potential energy has been accumulated * and the total force buffer is complete. @@ -252,6 +260,9 @@ class ForceProviders //! Prepares all providers for \p step, and returns whether any needs the potential energy. bool prepareStep(int64_t step, bool isCheckpointingStep) const; + //! Whether any provider asks to stop the simulation. + bool requestsStop() const; + //! Forwards applyAfterPotentialEnergy() to all providers. void applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, diff --git a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp index 8251a0394d..73ed672371 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp +++ b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp @@ -268,6 +268,11 @@ catch (const std::exception& ex) std::string("An error occurred while PLUMED was calculating the forces\n:") + ex.what())); } +bool PlumedForceProvider::requestsStop() const +{ + return plumedWantsToStop_ != 0; +} + void PlumedForceProvider::checkReplicaExchangeBias() { if (replex_) diff --git a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h index 617b81f660..90a8bcd29a 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h +++ b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h @@ -91,6 +91,8 @@ class PlumedForceProvider final : public IForceProvider * \param[in] isCheckpointingStep whether GROMACS writes a checkpoint at \p step */ bool prepareStep(int64_t step, bool isCheckpointingStep) override; + //! Whether PLUMED asked to stop the simulation during the last calculation + bool requestsStop() const override; /*! \brief Passes the potential energy to PLUMED and applies energy biases * \param[in] energyWasComputed whether the potential energy is valid this step * \param[in,out] potentialEnergy the MD potential energy diff --git a/patches/gromacs-2026.0.diff/src/gromacs/mdrun/runner.cpp b/patches/gromacs-2026.0.diff/src/gromacs/mdrun/runner.cpp index a29da7e600..9c34449e8b 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/mdrun/runner.cpp +++ b/patches/gromacs-2026.0.diff/src/gromacs/mdrun/runner.cpp @@ -138,6 +138,7 @@ #include "gromacs/mdtypes/fcdata.h" #include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/group.h" +#include "gromacs/mdtypes/iforceprovider.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/interaction_const.h" #include "gromacs/mdtypes/locality.h" @@ -2351,6 +2352,37 @@ int Mdrunner::mdrunner() simulatorBuilder.add(SimulatorStateData( globalState.get(), localState, &observablesHistory, &enerd, &ekind)); simulatorBuilder.add(std::move(membedHolder)); + /* PLUMED */ + if (fr->forceProviders != nullptr && isSimulationMainRank) + { + /* PLUMED can ask to stop the simulation, e.g. COMMITTOR. As a stop condition, the + request reaches all ranks and the simulations that share state, in both the + legacy and the modular simulator. The simulation stops after the next step + with global communication. The request is the same on all ranks, so only + the main rank signals it, which also keeps the summed signal small. + The modular simulator stops at the next neighbor-search step: when asked to + stop at another step, it clears the task queue it is running, and crashes. */ + const StopSignal plumedStopSignal = useModularSimulator ? StopSignal::stopAtNextNSStep + : StopSignal::stopImmediately; + stopHandlerBuilder_->registerStopCondition( + [forceProviders = fr->forceProviders, + plumedStopSignal, + fplog, + reported = std::make_shared(false)]() + { + if (!forceProviders->requestsStop()) + { + return StopSignal::noSignal; + } + if (fplog != nullptr && !*reported) + { + fprintf(fplog, "\n\nPLUMED asked to stop the simulation\n\n"); + *reported = true; + } + return plumedStopSignal; + }); + } + /* END PLUMED */ simulatorBuilder.add(std::move(stopHandlerBuilder_)); simulatorBuilder.add(SimulatorConfig(mdrunOptions, startingBehavior, &runScheduleWork)); diff --git a/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.cpp b/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.cpp index 26fe70f5fd..f7a158a438 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.cpp +++ b/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.cpp @@ -109,6 +109,18 @@ bool ForceProviders::prepareStep(int64_t step, bool isCheckpointingStep) const return needed; } +bool ForceProviders::requestsStop() const +{ + for (const auto& provider : impl_->providers_) + { + if (provider.first->requestsStop()) + { + return true; + } + } + return false; +} + void ForceProviders::applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef force, diff --git a/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.h b/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.h index 3084fb526b..3bf75f50e0 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.h +++ b/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.h @@ -197,6 +197,14 @@ class IForceProvider */ virtual bool prepareStep(int64_t /*step*/, bool /*isCheckpointingStep*/) { return false; } + /*! \brief + * Whether this provider asks to stop the simulation, after the forces it computed last. + * + * Used by PLUMED, e.g. when COMMITTOR finds the system in a basin. The MD loop gets the + * request through a stop condition registered with the StopHandlerBuilder. + */ + virtual bool requestsStop() const { return false; } + /*! \brief * Optional callback after the potential energy has been accumulated * and the total force buffer is complete. @@ -256,6 +264,9 @@ class ForceProviders //! Prepares all providers for \p step, and returns whether any needs the potential energy. bool prepareStep(int64_t step, bool isCheckpointingStep) const; + //! Whether any provider asks to stop the simulation. + bool requestsStop() const; + //! Forwards applyAfterPotentialEnergy() to all providers. void applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, From 226beed855dd24a3e4c1a8e8bca99e154c1a7e3b Mon Sep 17 00:00:00 2001 From: Michele Invernizzi Date: Wed, 14 Oct 2026 16:44:21 +0200 Subject: [PATCH 03/10] GROMACS 2025/2026 patches: spread PLUMED forces on virtual sites in the energy path When PLUMED needs the potential energy, it runs at the end of do_force(), after GROMACS has spread the forces on the virtual sites and cleared them. The forces PLUMED applied on virtual sites were then lost. They are now spread like the other forces with an explicit virial in postProcessForces(), on the steps when the atoms that PLUMED uses include virtual sites: applyAfterPotentialEnergy() now returns whether they do. That list is the same on all ranks, so that all take the same decision. The extra spread costs about 45 us per step (2% of a GPU step for a 68k-atom TIP4P box), so inputs that only bias the energy, e.g. OPES multithermal, do not pay for it. Tested on TIP4P water with a restraint between two M sites: adding ENERGY to the input changed the forces on those waters by 6.9e3 kJ/mol/nm, now by 0.001 on 1 rank and by 0 on 2 ranks. With the restraint between two oxygens the forces do not change and there is no extra spread. Co-Authored-By: Claude Opus 5.5 --- CHANGES/v2.10.md | 3 ++ .../applied_forces/plumed/plumedOptions.cpp | 13 +++++++++ .../applied_forces/plumed/plumedOptions.h | 3 ++ .../plumed/plumedforceprovider.cpp | 29 +++++++++++++++++-- .../plumed/plumedforceprovider.h | 7 ++++- .../src/gromacs/mdlib/sim_util.cpp | 20 ++++++++++++- .../src/gromacs/mdtypes/iforceprovider.cpp | 8 +++-- .../src/gromacs/mdtypes/iforceprovider.h | 10 +++++-- .../applied_forces/plumed/plumedOptions.cpp | 13 +++++++++ .../applied_forces/plumed/plumedOptions.h | 3 ++ .../plumed/plumedforceprovider.cpp | 28 ++++++++++++++++-- .../plumed/plumedforceprovider.h | 8 ++++- .../src/gromacs/mdlib/sim_util.cpp | 20 ++++++++++++- .../src/gromacs/mdtypes/iforceprovider.cpp | 8 +++-- .../src/gromacs/mdtypes/iforceprovider.h | 10 +++++-- 15 files changed, 163 insertions(+), 20 deletions(-) diff --git a/CHANGES/v2.10.md b/CHANGES/v2.10.md index d3994244e9..85d961ef01 100644 --- a/CHANGES/v2.10.md +++ b/CHANGES/v2.10.md @@ -121,3 +121,6 @@ Before switching to version 2.10, users are invited to carefully read the follow 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. diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.cpp b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.cpp index 98a9ba2018..d49180c9c3 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.cpp +++ b/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) diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h index 310970a9a9..cabd1593b9 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h @@ -44,6 +44,7 @@ #include #include +#include #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 isVirtualSite_; }; class PlumedOptionProvider diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp index db075fb0f8..18f8f2b644 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp @@ -41,6 +41,8 @@ #include "plumedforceprovider.h" +#include + #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()),replex_(options.replex_) +try : plumed_(std::make_unique()), + 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: @@ -242,6 +246,24 @@ catch (const std::exception& ex) 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 { @@ -289,7 +311,7 @@ void PlumedForceProvider::checkReplicaExchangeBias() } } -void PlumedForceProvider::applyAfterPotentialEnergy(bool energyWasComputed, +bool PlumedForceProvider::applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef force, tensor virial) @@ -297,7 +319,7 @@ try { if (!plumedNeedsEnergy_) { - return; + return false; } if (!energyWasComputed || potentialEnergy == nullptr) { @@ -322,6 +344,7 @@ try checkReplicaExchangeBias(); msmul(plumedVir_, 0.5, virial); + return usesVirtualSites(); } catch (const std::exception& ex) { diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h index 9741993169..5667485e0e 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h @@ -43,6 +43,7 @@ #include #include +#include #include "gromacs/mdtypes/iforceprovider.h" @@ -82,17 +83,21 @@ class PlumedForceProvider final : public IForceProvider ForceProviderOutput* forceProviderOutput) override; bool prepareStep(int64_t step, bool isCheckpointingStep) override; bool requestsStop() const override; - void applyAfterPotentialEnergy(bool energyWasComputed, + bool applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef force, tensor virial) override; private: void checkReplicaExchangeBias(); + //! Whether the atoms that PLUMED used on this step include virtual sites + bool usesVirtualSites(); std::unique_ptr plumed_; int plumedAPIversion_; bool replex_; + //! Whether each atom is a virtual site, empty if there are none + std::vector isVirtualSite_; int plumedNeedsEnergy_ = 0; int plumedWantsToStop_ = 0; matrix plumedVir_ = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } }; diff --git a/patches/gromacs-2025.0.diff/src/gromacs/mdlib/sim_util.cpp b/patches/gromacs-2025.0.diff/src/gromacs/mdlib/sim_util.cpp index 235a91fb54..58a9f585cb 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/mdlib/sim_util.cpp +++ b/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 */ diff --git a/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.cpp b/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.cpp index f7a158a438..c53d8c4879 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.cpp +++ b/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.cpp @@ -121,13 +121,17 @@ bool ForceProviders::requestsStop() const return false; } -void ForceProviders::applyAfterPotentialEnergy(bool energyWasComputed, +bool ForceProviders::applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef force, tensor virial) const { + bool forcesOnVirtualSites = false; for (auto& provider : impl_->providers_) { - provider.first->applyAfterPotentialEnergy(energyWasComputed, potentialEnergy, force, virial); + forcesOnVirtualSites = provider.first->applyAfterPotentialEnergy( + energyWasComputed, potentialEnergy, force, virial) + || forcesOnVirtualSites; } + return forcesOnVirtualSites; } diff --git a/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.h b/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.h index 68c32e5423..e8adf8c1cb 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.h +++ b/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.h @@ -212,12 +212,15 @@ class IForceProvider * \param[in,out] potentialEnergy Pointer to the potential energy (F_EPOT) * \param[in,out] force Complete force array (home atoms) * \param[in,out] virial Force virial tensor + * \returns Whether the provider might have applied forces on virtual sites, + * which then have to be spread on the constructing atoms */ - virtual void applyAfterPotentialEnergy(bool /*energyWasComputed*/, + virtual bool applyAfterPotentialEnergy(bool /*energyWasComputed*/, real* /*potentialEnergy*/, ArrayRef /*force*/, tensor /*virial*/) { + return false; } protected: @@ -263,8 +266,9 @@ class ForceProviders //! Whether any provider asks to stop the simulation. bool requestsStop() const; - //! Forwards applyAfterPotentialEnergy() to all providers. - void applyAfterPotentialEnergy(bool energyWasComputed, + //! Forwards applyAfterPotentialEnergy() to all providers, and returns whether any might have + //! applied forces on virtual sites. + bool applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef force, tensor virial) const; diff --git a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.cpp b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.cpp index aaa6c97c58..cdd642ee50 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.cpp +++ b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.cpp @@ -43,6 +43,9 @@ #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" #include "gromacs/utility/mpicomm.h" #include "gromacs/utility/vec.h" @@ -53,6 +56,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) diff --git a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h index 46ef96b34c..c754414ffa 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h +++ b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h @@ -44,6 +44,7 @@ #include #include +#include #include "gromacs/utility/real.h" @@ -68,6 +69,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 isVirtualSite_; }; class PlumedOptionProvider diff --git a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp index 73ed672371..2258f1590b 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp +++ b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp @@ -41,6 +41,7 @@ #include "plumedforceprovider.h" +#include #include #include "gromacs/domdec/domdec.h" @@ -75,7 +76,9 @@ namespace gmx PlumedForceProvider::~PlumedForceProvider() = default; PlumedForceProvider::PlumedForceProvider(const PlumedOptions& options) -try : plumed_(std::make_unique()), replex_(options.replex_) +try : plumed_(std::make_unique()), + 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: @@ -287,6 +290,24 @@ void PlumedForceProvider::checkReplicaExchangeBias() } } +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 { @@ -315,7 +336,7 @@ catch (const std::exception& ex) std::string("An error occurred while PLUMED was preparing the step\n:") + ex.what())); } -void PlumedForceProvider::applyAfterPotentialEnergy(bool energyWasComputed, +bool PlumedForceProvider::applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef force, tensor virial) @@ -323,7 +344,7 @@ try { if (!plumedNeedsEnergy_) { - return; + return false; } if (!energyWasComputed || potentialEnergy == nullptr) { @@ -348,6 +369,7 @@ try checkReplicaExchangeBias(); msmul(plumedVir_, 0.5, virial); + return usesVirtualSites(); } catch (const std::exception& ex) { diff --git a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h index 90a8bcd29a..b05d6ffd6b 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h +++ b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h @@ -43,6 +43,7 @@ #include #include +#include #include "gromacs/mdtypes/iforceprovider.h" @@ -98,8 +99,9 @@ class PlumedForceProvider final : public IForceProvider * \param[in,out] potentialEnergy the MD potential energy * \param[in,out] force the complete force array * \param[in,out] virial the force virial + * \returns whether the atoms that PLUMED used include virtual sites */ - void applyAfterPotentialEnergy(bool energyWasComputed, + bool applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef force, tensor virial) override; @@ -107,10 +109,14 @@ class PlumedForceProvider final : public IForceProvider private: //! Throws if PLUMED computed a bias while replica exchange is active void checkReplicaExchangeBias(); + //! Whether the atoms that PLUMED used on this step include virtual sites + bool usesVirtualSites(); std::unique_ptr plumed_; int plumedAPIversion_; bool replex_; + //! Whether each atom is a virtual site, empty if there are none + std::vector isVirtualSite_; std::optional> globalAtomIndices_; int plumedNeedsEnergy_ = 0; diff --git a/patches/gromacs-2026.0.diff/src/gromacs/mdlib/sim_util.cpp b/patches/gromacs-2026.0.diff/src/gromacs/mdlib/sim_util.cpp index a089006882..36e9de6e67 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/mdlib/sim_util.cpp +++ b/patches/gromacs-2026.0.diff/src/gromacs/mdlib/sim_util.cpp @@ -2740,8 +2740,26 @@ void do_force(FILE* fplog, { forceForPlumed = forceView->forceMtsCombined(); } - fr->forceProviders->applyAfterPotentialEnergy( + const bool plumedForcesOnVirtualSites = fr->forceProviders->applyAfterPotentialEnergy( stepWork.computeEnergy, &enerd->term[InteractionFunction::PotentialEnergy], 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 */ diff --git a/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.cpp b/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.cpp index f7a158a438..c53d8c4879 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.cpp +++ b/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.cpp @@ -121,13 +121,17 @@ bool ForceProviders::requestsStop() const return false; } -void ForceProviders::applyAfterPotentialEnergy(bool energyWasComputed, +bool ForceProviders::applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef force, tensor virial) const { + bool forcesOnVirtualSites = false; for (auto& provider : impl_->providers_) { - provider.first->applyAfterPotentialEnergy(energyWasComputed, potentialEnergy, force, virial); + forcesOnVirtualSites = provider.first->applyAfterPotentialEnergy( + energyWasComputed, potentialEnergy, force, virial) + || forcesOnVirtualSites; } + return forcesOnVirtualSites; } diff --git a/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.h b/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.h index 3bf75f50e0..64a7fdb20f 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.h +++ b/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.h @@ -216,12 +216,15 @@ class IForceProvider * \param[in,out] potentialEnergy Pointer to the potential energy (F_EPOT) * \param[in,out] force Complete force array (home atoms) * \param[in,out] virial Force virial tensor + * \returns Whether the provider might have applied forces on virtual sites, + * which then have to be spread on the constructing atoms */ - virtual void applyAfterPotentialEnergy(bool /*energyWasComputed*/, + virtual bool applyAfterPotentialEnergy(bool /*energyWasComputed*/, real* /*potentialEnergy*/, ArrayRef /*force*/, tensor /*virial*/) { + return false; } protected: @@ -267,8 +270,9 @@ class ForceProviders //! Whether any provider asks to stop the simulation. bool requestsStop() const; - //! Forwards applyAfterPotentialEnergy() to all providers. - void applyAfterPotentialEnergy(bool energyWasComputed, + //! Forwards applyAfterPotentialEnergy() to all providers, and returns whether any might have + //! applied forces on virtual sites. + bool applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef force, tensor virial) const; From 8037947005eedaf83cdc4796ebfc7911494e4c9c Mon Sep 17 00:00:00 2001 From: Michele Invernizzi Date: Wed, 14 Oct 2026 16:44:21 +0200 Subject: [PATCH 04/10] GROMACS 2025/2026 patches: stop if PLUMED needs the energy on an MTS fast step With multiple time stepping, the slow forces and their energies are computed only every mts-level2-factor steps. If PLUMED needed the energy on another step, it got an energy without the slow terms: in a TIP4P test with mts-level2-factor 2, ENERGY was off by 650-1250 kJ/mol on the odd steps. mdrun now stops with an error, which says to need the energy only on the slow steps (e.g. no bias on ENERGY) or not to use multiple time stepping. Co-Authored-By: Claude Opus 5.5 --- CHANGES/v2.10.md | 3 +++ patches/gromacs-2025.0.diff/src/gromacs/mdrun/md.cpp | 12 ++++++++++++ patches/gromacs-2026.0.diff/src/gromacs/mdrun/md.cpp | 12 ++++++++++++ 3 files changed, 27 insertions(+) diff --git a/CHANGES/v2.10.md b/CHANGES/v2.10.md index 85d961ef01..0d6d59c1e7 100644 --- a/CHANGES/v2.10.md +++ b/CHANGES/v2.10.md @@ -124,3 +124,6 @@ Before switching to version 2.10, users are invited to carefully read the follow - 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. diff --git a/patches/gromacs-2025.0.diff/src/gromacs/mdrun/md.cpp b/patches/gromacs-2025.0.diff/src/gromacs/mdrun/md.cpp index 60612f6ff4..238081cdcf 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/mdrun/md.cpp +++ b/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 */ diff --git a/patches/gromacs-2026.0.diff/src/gromacs/mdrun/md.cpp b/patches/gromacs-2026.0.diff/src/gromacs/mdrun/md.cpp index ddd09af430..774a440ea5 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/mdrun/md.cpp +++ b/patches/gromacs-2026.0.diff/src/gromacs/mdrun/md.cpp @@ -1129,6 +1129,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 */ From 5097089c2428aa61ab7e3e9be1c49b387527ac75 Mon Sep 17 00:00:00 2001 From: Michele Invernizzi Date: Wed, 14 Oct 2026 16:44:21 +0200 Subject: [PATCH 05/10] GROMACS 2025/2026 patches: keep the type of the errors raised on purpose The force provider rethrew every exception as an InternalError, so GROMACS reported e.g. replica exchange with a PLUMED bias as "Internal error (bug)". GROMACS exceptions now pass unchanged, so this is reported as "Feature not implemented"; PLUMED errors are still wrapped as before. Co-Authored-By: Claude Opus 5.5 --- CHANGES/v2.10.md | 2 ++ .../plumed/plumedforceprovider.cpp | 25 +++++++++++++++++++ .../plumed/plumedforceprovider.cpp | 25 +++++++++++++++++++ 3 files changed, 52 insertions(+) diff --git a/CHANGES/v2.10.md b/CHANGES/v2.10.md index 0d6d59c1e7..3a493951f0 100644 --- a/CHANGES/v2.10.md +++ b/CHANGES/v2.10.md @@ -127,3 +127,5 @@ Before switching to version 2.10, users are invited to carefully read the follow - 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. diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp index 18f8f2b644..b310ab53ea 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp @@ -151,6 +151,11 @@ try : plumed_(std::make_unique()), 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( @@ -165,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( @@ -240,6 +250,11 @@ try 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( @@ -286,6 +301,11 @@ 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( @@ -346,6 +366,11 @@ try 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) { GMX_THROW(InternalError( diff --git a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp index 2258f1590b..e94722ff02 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp +++ b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp @@ -165,6 +165,11 @@ try : plumed_(std::make_unique()), 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( @@ -184,6 +189,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( @@ -265,6 +275,11 @@ try 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( @@ -330,6 +345,11 @@ 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( @@ -371,6 +391,11 @@ try 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) { GMX_THROW(InternalError( From b5d1521242626cbd583cdde737c4e8d803e82f4a Mon Sep 17 00:00:00 2001 From: Michele Invernizzi Date: Wed, 14 Oct 2026 16:44:21 +0200 Subject: [PATCH 06/10] GROMACS 2025 patch: declare gmx_multisim_t as a struct, as it is defined Co-Authored-By: Claude Opus 5.5 --- .../src/gromacs/applied_forces/plumed/plumedMDModule.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedMDModule.cpp b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedMDModule.cpp index d3a56f82f0..a49f579c17 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedMDModule.cpp +++ b/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 { From 040a5b446d4c065b6bc983908cfd2979ebb4119a Mon Sep 17 00:00:00 2001 From: Michele Invernizzi Date: Wed, 14 Oct 2026 16:44:21 +0200 Subject: [PATCH 07/10] GROMACS 2025/2026 patches: do not update PLUMED again in the shell relaxation With shell (polarizable) particles, relax_shell_flexcon() calls do_force() again for each shell relaxation iteration, passing step 1. PLUMED then calculated and updated on each of them, as step 1: e.g. PRINT printed extra lines at step 1, and actions updating at every step (e.g. PACE=1) were updated several times per step. PLUMED now uses the step from prepareStep(), and only recalculates its forces (performCalcNoUpdate) when called again on the same step. This is limited to MD, where prepareStep() gives the step: e.g. the frames of a rerun can all have the same step, and the conjugate gradient minimizer passes step -1. On a SW water dimer, COLVAR now has one line per step. Co-Authored-By: Claude Opus 5.5 --- CHANGES/v2.10.md | 3 +++ .../plumed/plumedforceprovider.cpp | 17 +++++++++++++---- .../applied_forces/plumed/plumedforceprovider.h | 6 ++++++ .../plumed/plumedforceprovider.cpp | 17 +++++++++++++---- .../applied_forces/plumed/plumedforceprovider.h | 6 ++++++ 5 files changed, 41 insertions(+), 8 deletions(-) diff --git a/CHANGES/v2.10.md b/CHANGES/v2.10.md index 3a493951f0..b0137909a6 100644 --- a/CHANGES/v2.10.md +++ b/CHANGES/v2.10.md @@ -129,3 +129,6 @@ Before switching to version 2.10, users are invited to carefully read the follow 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: 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. diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp index b310ab53ea..3fc8b236e4 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp @@ -188,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)) { @@ -215,7 +220,10 @@ try } preparedStep_.reset(); - plumedWantsToStop_ = 0; + if (!repeatedCall_) + { + plumedWantsToStop_ = 0; + } plumed_->cmd("setStopFlag", &plumedWantsToStop_); if (!prepareStepIsCalled_) @@ -243,7 +251,7 @@ try clear_mat(plumedVir_); plumed_->cmd("setVirial", &plumedVir_[0][0]); - plumed_->cmd("performCalc", nullptr); + plumed_->cmd(repeatedCall_ ? "performCalcNoUpdate" : "performCalc", nullptr); checkReplicaExchangeBias(); @@ -289,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 @@ -359,7 +368,7 @@ 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(); diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h index 5667485e0e..68f99759cf 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h @@ -107,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 calculatedStep_; + //! Whether PLUMED already calculated this MD step (shell relaxation): it is then not updated again + bool repeatedCall_ = false; }; } // namespace gmx diff --git a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp index e94722ff02..a405da34fa 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp +++ b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp @@ -206,7 +206,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 - 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 (globalAtomIndices_.has_value()) { @@ -233,7 +238,10 @@ try } preparedStep_.reset(); - plumedWantsToStop_ = 0; + if (!repeatedCall_) + { + plumedWantsToStop_ = 0; + } plumed_->cmd("setStopFlag", &plumedWantsToStop_); if (!prepareStepIsCalled_) @@ -268,7 +276,7 @@ try // in the original patch do_force() was called HERE // Do the work - plumed_->cmd("performCalc", nullptr); + plumed_->cmd(repeatedCall_ ? "performCalcNoUpdate" : "performCalc", nullptr); checkReplicaExchangeBias(); @@ -333,6 +341,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 @@ -384,7 +393,7 @@ try // Same virial convention as the legacy patch: hand PLUMED 2*virial, then replace it msmul(virial, 2.0, plumedVir_); plumed_->cmd("setVirial", &plumedVir_[0][0]); - plumed_->cmd("performCalc", nullptr); + plumed_->cmd(repeatedCall_ ? "performCalcNoUpdate" : "performCalc", nullptr); checkReplicaExchangeBias(); diff --git a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h index b05d6ffd6b..c58a5a7a93 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h +++ b/patches/gromacs-2026.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.h @@ -128,6 +128,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 calculatedStep_; + //! Whether PLUMED already calculated this MD step (shell relaxation): it is then not updated again + bool repeatedCall_ = false; }; } // namespace gmx From 7da918a4a2343216ba58b0721b1d47388cd2ac8a Mon Sep 17 00:00:00 2001 From: Michele Invernizzi Date: Wed, 14 Oct 2026 17:37:48 +0200 Subject: [PATCH 08/10] GROMACS 2025/2026 patches: honour stop requests in reruns In a rerun, GROMACS sets isLastStep from the stop handler at the start of each frame, but then overwrites it with the result of reading the next frame: a stop request, e.g. from COMMITTOR through the stop condition of the PLUMED module, was ignored and the rerun went through all the frames. rerun.cpp is now patched to keep it, as in the GROMACS 2020-2024 patches. Found by a regression test: a rerun with COMMITTOR committing at the first frame went through all 11 frames. Co-Authored-By: Claude Opus 5.5 --- CHANGES/v2.10.md | 1 + .../src/gromacs/mdrun/rerun.cpp | 943 ++++++++++++++++++ .../src/gromacs/mdrun/rerun.cpp.preplumed | 942 +++++++++++++++++ .../src/gromacs/mdrun/rerun.cpp | 942 +++++++++++++++++ .../src/gromacs/mdrun/rerun.cpp.preplumed | 941 +++++++++++++++++ 5 files changed, 3769 insertions(+) create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/mdrun/rerun.cpp create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/mdrun/rerun.cpp.preplumed create mode 100644 patches/gromacs-2026.0.diff/src/gromacs/mdrun/rerun.cpp create mode 100644 patches/gromacs-2026.0.diff/src/gromacs/mdrun/rerun.cpp.preplumed diff --git a/CHANGES/v2.10.md b/CHANGES/v2.10.md index b0137909a6..04d0a51440 100644 --- a/CHANGES/v2.10.md +++ b/CHANGES/v2.10.md @@ -129,6 +129,7 @@ Before switching to version 2.10, users are invited to carefully read the follow 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. diff --git a/patches/gromacs-2025.0.diff/src/gromacs/mdrun/rerun.cpp b/patches/gromacs-2025.0.diff/src/gromacs/mdrun/rerun.cpp new file mode 100644 index 0000000000..4113c5ec89 --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/mdrun/rerun.cpp @@ -0,0 +1,943 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2018- The GROMACS Authors + * and the project initiators Erik Lindahl, Berk Hess and David van der Spoel. + * Consult the AUTHORS/COPYING files and https://www.gromacs.org for details. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * https://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at https://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out https://www.gromacs.org. + */ +/*! \internal \file + * + * \brief Implements the loop for simulation reruns + * + * \author Pascal Merz + * \ingroup module_mdrun + */ +#include "gmxpre.h" + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "gromacs/applied_forces/awh/awh.h" +#include "gromacs/commandline/filenm.h" +#include "gromacs/compat/pointers.h" +#include "gromacs/domdec/collect.h" +#include "gromacs/domdec/dlbtiming.h" +#include "gromacs/domdec/domdec.h" +#include "gromacs/domdec/domdec_network.h" +#include "gromacs/domdec/domdec_struct.h" +#include "gromacs/domdec/localtopologychecker.h" +#include "gromacs/domdec/mdsetup.h" +#include "gromacs/domdec/partition.h" +#include "gromacs/essentialdynamics/edsam.h" +#include "gromacs/ewald/pme_load_balancing.h" +#include "gromacs/ewald/pme_pp.h" +#include "gromacs/fileio/enxio.h" +#include "gromacs/fileio/trxio.h" +#include "gromacs/gmxlib/network.h" +#include "gromacs/gmxlib/nrnb.h" +#include "gromacs/gpu_utils/gpu_utils.h" +#include "gromacs/listed_forces/listed_forces.h" +#include "gromacs/listed_forces/listed_forces_gpu.h" +#include "gromacs/math/functions.h" +#include "gromacs/math/utilities.h" +#include "gromacs/math/vec.h" +#include "gromacs/math/vectypes.h" +#include "gromacs/mdlib/checkpointhandler.h" +#include "gromacs/mdlib/constr.h" +#include "gromacs/mdlib/ebin.h" +#include "gromacs/mdlib/enerdata_utils.h" +#include "gromacs/mdlib/energyoutput.h" +#include "gromacs/mdlib/expanded.h" +#include "gromacs/mdlib/force.h" +#include "gromacs/mdlib/force_flags.h" +#include "gromacs/mdlib/forcerec.h" +#include "gromacs/mdlib/freeenergyparameters.h" +#include "gromacs/mdlib/md_support.h" +#include "gromacs/mdlib/mdatoms.h" +#include "gromacs/mdlib/mdoutf.h" +#include "gromacs/mdlib/membed.h" +#include "gromacs/mdlib/resethandler.h" +#include "gromacs/mdlib/sighandler.h" +#include "gromacs/mdlib/simulationsignal.h" +#include "gromacs/mdlib/stat.h" +#include "gromacs/mdlib/stophandler.h" +#include "gromacs/mdlib/tgroup.h" +#include "gromacs/mdlib/trajectory_writing.h" +#include "gromacs/mdlib/update.h" +#include "gromacs/mdlib/vcm.h" +#include "gromacs/mdlib/vsite.h" +#include "gromacs/mdrunutility/handlerestart.h" +#include "gromacs/mdrunutility/multisim.h" +#include "gromacs/mdrunutility/printtime.h" +#include "gromacs/mdtypes/awh_history.h" +#include "gromacs/mdtypes/awh_params.h" +#include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/df_history.h" +#include "gromacs/mdtypes/energyhistory.h" +#include "gromacs/mdtypes/forcebuffers.h" +#include "gromacs/mdtypes/forcerec.h" +#include "gromacs/mdtypes/group.h" +#include "gromacs/mdtypes/inputrec.h" +#include "gromacs/mdtypes/interaction_const.h" +#include "gromacs/mdtypes/md_enums.h" +#include "gromacs/mdtypes/mdatom.h" +#include "gromacs/mdtypes/mdrunoptions.h" +#include "gromacs/mdtypes/multipletimestepping.h" +#include "gromacs/mdtypes/observableshistory.h" +#include "gromacs/mdtypes/observablesreducer.h" +#include "gromacs/mdtypes/simulation_workload.h" +#include "gromacs/mdtypes/state.h" +#include "gromacs/mimic/utilities.h" +#include "gromacs/pbcutil/pbc.h" +#include "gromacs/pulling/output.h" +#include "gromacs/pulling/pull.h" +#include "gromacs/swap/swapcoords.h" +#include "gromacs/taskassignment/include/gromacs/taskassignment/decidesimulationworkload.h" +#include "gromacs/timing/wallcycle.h" +#include "gromacs/timing/walltime_accounting.h" +#include "gromacs/topology/atoms.h" +#include "gromacs/topology/idef.h" +#include "gromacs/topology/ifunc.h" +#include "gromacs/topology/mtop_util.h" +#include "gromacs/topology/topology.h" +#include "gromacs/trajectory/trajectoryframe.h" +#include "gromacs/utility/arrayref.h" +#include "gromacs/utility/basedefinitions.h" +#include "gromacs/utility/cstringutil.h" +#include "gromacs/utility/enumerationhelpers.h" +#include "gromacs/utility/exceptions.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/gmxassert.h" +#include "gromacs/utility/logger.h" +#include "gromacs/utility/real.h" + +#include "legacysimulator.h" +#include "replicaexchange.h" +#include "shellfc.h" + +struct gmx_edsam; +struct gmx_mdoutf; +struct gmx_shellfc_t; + +using gmx::SimulationSignaller; +using gmx::VirtualSitesHandler; + +/*! \brief Copy the state from \p rerunFrame to \p globalState and, if requested, construct vsites + * + * \param[in] rerunFrame The trajectory frame to compute energy/forces for + * \param[in,out] globalState The global state container + * \param[in] constructVsites When true, vsite coordinates are constructed + * \param[in] vsite Vsite setup, can be nullptr when \p constructVsites = false + */ +static void prepareRerunState(const t_trxframe& rerunFrame, + t_state* globalState, + bool constructVsites, + const VirtualSitesHandler* vsite) +{ + auto x = makeArrayRef(globalState->x); + auto rerunX = arrayRefFromArray(reinterpret_cast(rerunFrame.x), globalState->numAtoms()); + std::copy(rerunX.begin(), rerunX.end(), x.begin()); + copy_mat(rerunFrame.box, globalState->box); + + if (constructVsites) + { + GMX_ASSERT(vsite, "Need valid vsite for constructing vsites"); + + vsite->construct(globalState->x, globalState->v, globalState->box, gmx::VSiteOperation::PositionsAndVelocities); + } +} + +void gmx::LegacySimulator::do_rerun() +{ + // TODO Historically, the EM and MD "integrators" used different + // names for the t_inputrec *parameter, but these must have the + // same name, now that it's a member of a struct. We use this ir + // alias to avoid a large ripple of nearly useless changes. + // t_inputrec is being replaced by IMdpOptionsProvider, so this + // will go away eventually. + const t_inputrec* ir = inputRec_; + double t; + bool isLastStep = false; + bool doFreeEnergyPerturbation = false; + unsigned int force_flags; + tensor force_vir, shake_vir, total_vir, pres; + t_trxstatus* status = nullptr; + rvec mu_tot; + t_trxframe rerun_fr; + ForceBuffers f; + gmx_global_stat_t gstat; + gmx_shellfc_t* shellfc; + + double cycles; + + GMX_RELEASE_ASSERT(cr_->dd == nullptr || !ddUsesUpdateGroups(*cr_->dd), + "Update groups are not supported with rerun"); + + SimulationSignals signals; + // Most global communication stages don't propagate mdrun + // signals, and will use this object to achieve that. + SimulationSignaller nullSignaller(nullptr, nullptr, nullptr, false, false); + + GMX_LOG(mdLog_.info) + .asParagraph() + .appendText( + "Note that it is planned that the command gmx mdrun -rerun will " + "be available in a different form in a future version of GROMACS, " + "e.g. gmx rerun -f."); + + if (ir->efep != FreeEnergyPerturbationType::No + && (mdAtoms_->mdatoms()->nMassPerturbed > 0 || (constr_ && constr_->havePerturbedConstraints()))) + { + gmx_fatal(FARGS, + "Perturbed masses or constraints are not supported by rerun. " + "Either make a .tpr without mass and constraint perturbation, " + "or use GROMACS 2018.4, 2018.5 or later 2018 version."); + } + if (ir->bExpanded) + { + gmx_fatal(FARGS, "Expanded ensemble not supported by rerun."); + } + if (ir->bSimTemp) + { + gmx_fatal(FARGS, "Simulated tempering not supported by rerun."); + } + if (ir->bDoAwh) + { + gmx_fatal(FARGS, "AWH not supported by rerun."); + } + if (replExParams_.exchangeInterval > 0) + { + gmx_fatal(FARGS, "Replica exchange not supported by rerun."); + } + if (opt2bSet("-ei", nFile_, fnm_) || observablesHistory_->edsamHistory != nullptr) + { + gmx_fatal(FARGS, "Essential dynamics not supported by rerun."); + } + if (ir->bIMD) + { + gmx_fatal(FARGS, "Interactive MD not supported by rerun."); + } + if (isMultiSim(ms_)) + { + gmx_fatal(FARGS, "Multiple simulations not supported by rerun."); + } + if (std::any_of(ir->opts.annealing, + ir->opts.annealing + ir->opts.ngtc, + [](SimulatedAnnealing i) { return i != SimulatedAnnealing::No; })) + { + gmx_fatal(FARGS, "Simulated annealing not supported by rerun."); + } + + /* Rerun can't work if an output file name is the same as the input file name. + * If this is the case, the user will get an error telling them what the issue is. + */ + if (strcmp(opt2fn("-rerun", nFile_, fnm_), opt2fn("-o", nFile_, fnm_)) == 0 + || strcmp(opt2fn("-rerun", nFile_, fnm_), opt2fn("-x", nFile_, fnm_)) == 0) + { + gmx_fatal(FARGS, + "When using mdrun -rerun, the name of the input trajectory file " + "%s cannot be identical to the name of an output file (whether " + "given explicitly with -o or -x, or by default)", + opt2fn("-rerun", nFile_, fnm_)); + } + + /* Settings for rerun */ + { + // TODO: Avoid changing inputrec (#3854) + auto* nonConstInputrec = const_cast(inputRec_); + nonConstInputrec->nstlist = 1; + nonConstInputrec->nstcalcenergy = 1; + nonConstInputrec->nstxout_compressed = 0; + } + int nstglobalcomm = 1; + const bool bNS = true; + + ObservablesReducer observablesReducer = observablesReducerBuilder_->build(); + + const SimulationGroups* groups = &topGlobal_.groups; + if (ir->eI == IntegrationAlgorithm::Mimic) + { + auto* nonConstGlobalTopology = const_cast(&topGlobal_); + nonConstGlobalTopology->intermolecularExclusionGroup = genQmmmIndices(topGlobal_); + } + int* fep_state = MAIN(cr_) ? &stateGlobal_->fep_state : nullptr; + gmx::ArrayRef lambda = MAIN(cr_) ? stateGlobal_->lambda : gmx::ArrayRef(); + initialize_lambdas( + fpLog_, ir->efep, ir->bSimTemp, *ir->fepvals, ir->simtempvals->temperatures, ekind_, MAIN(cr_), fep_state, lambda); + const bool simulationsShareState = false; + gmx_mdoutf* outf = init_mdoutf(fpLog_, + nFile_, + fnm_, + mdrunOptions_, + cr_, + outputProvider_, + mdModulesNotifiers_, + ir, + topGlobal_, + oenv_, + wallCycleCounters_, + StartingBehavior::NewSimulation, + simulationsShareState, + ms_); + gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), + topGlobal_, + *ir, + pullWork_, + mdoutf_get_fp_dhdl(outf), + true, + StartingBehavior::NewSimulation, + simulationsShareState, + mdModulesNotifiers_); + + gstat = global_stat_init(ir); + + /* Check for polarizable models and flexible constraints */ + shellfc = init_shell_flexcon(fpLog_, + topGlobal_, + constr_ ? constr_->numFlexibleConstraints() : 0, + ir->nstcalcenergy, + haveDDAtomOrdering(*cr_), + runScheduleWork_->simulationWork.useGpuPme); + + if (haveDDAtomOrdering(*cr_)) + { + // Local state only becomes valid now. + dd_init_local_state(*cr_->dd, stateGlobal_, state_); + + /* Distribute the charge groups over the nodes from the main node */ + dd_partition_system(fpLog_, + mdLog_, + ir->init_step, + cr_, + TRUE, + stateGlobal_, + topGlobal_, + *ir, + mdModulesNotifiers_, + imdSession_, + pullWork_, + state_, + &f, + mdAtoms_, + top_, + fr_, + virtualSites_, + constr_, + nrnb_, + nullptr, + FALSE); + } + else + { + /* Copy the pointer to the global state */ + state_ = stateGlobal_; + + mdAlgorithmsSetupAtomData( + cr_, *ir, topGlobal_, top_, fr_, &f, mdAtoms_, constr_, virtualSites_, shellfc); + } + + auto* mdatoms = mdAtoms_->mdatoms(); + fr_->longRangeNonbondeds->updateAfterPartition(*mdatoms); + + // NOTE: The global state is no longer used at this point. + // But state_global is still used as temporary storage space for writing + // the global state to file and potentially for replica exchange. + // (Global topology should persist.) + + update_mdatoms(mdatoms, state_->lambda[FreeEnergyPerturbationCouplingType::Mass]); + + if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl != 0) + { + doFreeEnergyPerturbation = true; + } + + int64_t step = ir->init_step; + int64_t step_rel = 0; + + { + int cglo_flags = CGLO_GSTAT; + bool bSumEkinhOld = false; + t_vcm* vcm = nullptr; + compute_globals(gstat, + cr_, + ir, + fr_, + ekind_, + makeConstArrayRef(state_->x), + makeConstArrayRef(state_->v), + state_->box, + mdatoms, + nrnb_, + vcm, + nullptr, + enerd_, + force_vir, + shake_vir, + total_vir, + pres, + &nullSignaller, + state_->box, + &bSumEkinhOld, + cglo_flags, + step, + &observablesReducer); + // Clean up after pre-step use of compute_globals() + observablesReducer.markAsReadyToReduce(); + } + + if (MAIN(cr_)) + { + fprintf(stderr, + "starting md rerun '%s', reading coordinates from" + " input trajectory '%s'\n\n", + *(topGlobal_.name), + opt2fn("-rerun", nFile_, fnm_)); + if (mdrunOptions_.verbose) + { + fprintf(stderr, + "Calculated time to finish depends on nsteps from " + "run input file,\nwhich may not correspond to the time " + "needed to process input trajectory.\n\n"); + } + fprintf(fpLog_, "\n"); + } + + walltime_accounting_start_time(wallTimeAccounting_); + wallcycle_start(wallCycleCounters_, WallCycleCounter::Run); + print_start(fpLog_, cr_, wallTimeAccounting_, "mdrun"); + + /*********************************************************** + * + * Loop over MD steps + * + ************************************************************/ + + if (constr_) + { + GMX_LOG(mdLog_.info) + .asParagraph() + .appendText("Simulations has constraints. Rerun does not recalculate constraints."); + } + + rerun_fr.natoms = 0; + if (MAIN(cr_)) + { + isLastStep = !read_first_frame(oenv_, &status, opt2fn("-rerun", nFile_, fnm_), &rerun_fr, TRX_NEED_X); + if (rerun_fr.natoms != topGlobal_.natoms) + { + gmx_fatal(FARGS, + "Number of atoms in trajectory (%d) does not match the " + "run input file (%d)\n", + rerun_fr.natoms, + topGlobal_.natoms); + } + + if (ir->pbcType != PbcType::No) + { + if (!rerun_fr.bBox) + { + gmx_fatal(FARGS, + "Rerun trajectory frame step %" PRId64 + " time %f " + "does not contain a box, while pbc is used", + rerun_fr.step, + rerun_fr.time); + } + if (max_cutoff2(ir->pbcType, rerun_fr.box) < gmx::square(fr_->rlist)) + { + gmx_fatal(FARGS, + "Rerun trajectory frame step %" PRId64 + " time %f " + "has too small box dimensions", + rerun_fr.step, + rerun_fr.time); + } + } + } + + GMX_LOG(mdLog_.info) + .asParagraph() + .appendText( + "Rerun does not report kinetic energy, total energy, temperature, virial and " + "pressure."); + + if (PAR(cr_)) + { + rerun_parallel_comm(cr_, &rerun_fr, &isLastStep); + } + + if (ir->pbcType != PbcType::No) + { + /* Set the shift vectors. + * Necessary here when have a static box different from the tpr box. + */ + calc_shifts(rerun_fr.box, fr_->shift_vec); + } + + auto stopHandler = stopHandlerBuilder_->getStopHandlerMD( + compat::not_null(&signals[eglsSTOPCOND]), + false, + MAIN(cr_), + ir->nstlist, + mdrunOptions_.reproducible, + nstglobalcomm, + mdrunOptions_.maximumHoursToRun, + ir->nstlist == 0, + fpLog_, + step, + bNS, + wallTimeAccounting_); + + // we don't do counter resetting in rerun - finish will always be valid + walltime_accounting_set_valid_finish(wallTimeAccounting_); + + const DDBalanceRegionHandler ddBalanceRegionHandler(cr_); + + /* and stop now if we should */ + isLastStep = (isLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps)); + while (!isLastStep) + { + wallcycle_start(wallCycleCounters_, WallCycleCounter::Step); + + if (rerun_fr.bStep) + { + step = rerun_fr.step; + step_rel = step - ir->init_step; + } + if (rerun_fr.bTime) + { + t = rerun_fr.time; + } + else + { + t = step; + } + + if (ir->efep != FreeEnergyPerturbationType::No && MAIN(cr_)) + { + if (rerun_fr.bLambda) + { + ir->fepvals->init_lambda_without_states = rerun_fr.lambda; + } + else + { + if (rerun_fr.bFepState) + { + state_->fep_state = rerun_fr.fep_state; + } + } + + stateGlobal_->lambda = currentLambdas(step, *(ir->fepvals), state_->fep_state); + } + + if (MAIN(cr_)) + { + const bool constructVsites = + ((virtualSites_ != nullptr) && mdrunOptions_.rerunConstructVsites); + if (constructVsites && haveDDAtomOrdering(*cr_)) + { + gmx_fatal(FARGS, + "Vsite recalculation with -rerun is not implemented with domain " + "decomposition, " + "use a single rank"); + } + prepareRerunState(rerun_fr, stateGlobal_, constructVsites, virtualSites_); + } + + isLastStep = isLastStep || stopHandler->stoppingAfterCurrentStep(bNS); + + if (haveDDAtomOrdering(*cr_)) + { + /* Repartition the domain decomposition */ + const bool bMainState = true; + dd_partition_system(fpLog_, + mdLog_, + step, + cr_, + bMainState, + stateGlobal_, + topGlobal_, + *ir, + mdModulesNotifiers_, + imdSession_, + pullWork_, + state_, + &f, + mdAtoms_, + top_, + fr_, + virtualSites_, + constr_, + nrnb_, + wallCycleCounters_, + mdrunOptions_.verbose); + } + + if (MAIN(cr_)) + { + EnergyOutput::printHeader(fpLog_, step, t); /* can we improve the information printed here? */ + } + + if (ir->efep != FreeEnergyPerturbationType::No) + { + update_mdatoms(mdatoms, state_->lambda[FreeEnergyPerturbationCouplingType::Mass]); + } + + fr_->longRangeNonbondeds->updateAfterPartition(*mdatoms); + + force_flags = (GMX_FORCE_STATECHANGED | GMX_FORCE_DYNAMICBOX | GMX_FORCE_ALLFORCES + | GMX_FORCE_VIRIAL | // TODO: Get rid of this once #2649 and #3400 are solved + GMX_FORCE_ENERGY | (doFreeEnergyPerturbation ? GMX_FORCE_DHDL : 0)); + + const int shellfcFlags = force_flags | (mdrunOptions_.verbose ? GMX_FORCE_ENERGY : 0); + const int legacyForceFlags = ((shellfc) ? shellfcFlags : force_flags) | GMX_FORCE_NS; + + gmx_edsam* const ed = nullptr; + + if (bNS) + { + if (fr_->listedForcesGpu) + { + fr_->listedForcesGpu->updateHaveInteractions(top_->idef); + } + runScheduleWork_->domainWork = setupDomainLifetimeWorkload( + *ir, *fr_, pullWork_, ed, *mdatoms, runScheduleWork_->simulationWork); + } + + + runScheduleWork_->stepWork = setupStepWorkload(legacyForceFlags, + ir->mtsLevels, + step, + runScheduleWork_->domainWork, + runScheduleWork_->simulationWork); + + if (shellfc) + { + /* Now is the time to relax the shells */ + relax_shell_flexcon(fpLog_, + cr_, + ms_, + mdrunOptions_.verbose, + enforcedRotation_, + step, + ir, + mdModulesNotifiers_, + imdSession_, + pullWork_, + bNS, + top_, + constr_, + enerd_, + state_->numAtoms(), + state_->x.arrayRefWithPadding(), + state_->v.arrayRefWithPadding(), + state_->box, + state_->lambda, + &state_->hist, + &f.view(), + force_vir, + *mdatoms, + fr_->longRangeNonbondeds.get(), + nrnb_, + wallCycleCounters_, + shellfc, + fr_, + *runScheduleWork_, + t, + mu_tot, + virtualSites_, + ddBalanceRegionHandler); + } + else + { + /* The coordinates (x) are shifted (to get whole molecules) + * in do_force. + * This is parallellized as well, and does communication too. + * Check comments in sim_util.c + */ + Awh* awh = nullptr; + + try + { + do_force(fpLog_, + cr_, + ms_, + *ir, + mdModulesNotifiers_, + awh, + enforcedRotation_, + imdSession_, + pullWork_, + step, + nrnb_, + wallCycleCounters_, + top_, + state_->box, + state_->x.arrayRefWithPadding(), + {}, + &state_->hist, + &f.view(), + force_vir, + mdatoms, + enerd_, + state_->lambda, + fr_, + *runScheduleWork_, + virtualSites_, + mu_tot, + t, + ed, + fr_->longRangeNonbondeds.get(), + ddBalanceRegionHandler); + } + catch (const gmx::InternalError&) + { + GMX_LOG(mdLog_.warning) + .asParagraph() + .appendText( + "Continuing with next frame after catching invalid force in " + "previous frame"); + }; + } + + /* Now we have the energies and forces corresponding to the + * coordinates at time t. + */ + { + const bool isCheckpointingStep = false; + const bool doRerun = true; + do_md_trajectory_writing(fpLog_, + cr_, + nFile_, + fnm_, + step, + step_rel, + t, + ir, + state_, + stateGlobal_, + observablesHistory_, + topGlobal_, + fr_, + outf, + energyOutput, + ekind_, + f.view().force(), + isCheckpointingStep, + doRerun, + isLastStep, + mdrunOptions_.writeConfout, + EkindataState::NotUsed); + } + + stopHandler->setSignal(); + + { + const bool doInterSimSignal = false; + const bool doIntraSimSignal = true; + bool bSumEkinhOld = false; + t_vcm* vcm = nullptr; + SimulationSignaller signaller(&signals, cr_, ms_, doInterSimSignal, doIntraSimSignal); + + int cglo_flags = CGLO_GSTAT | CGLO_ENERGY; + compute_globals(gstat, + cr_, + ir, + fr_, + ekind_, + makeConstArrayRef(state_->x), + makeConstArrayRef(state_->v), + state_->box, + mdatoms, + nrnb_, + vcm, + wallCycleCounters_, + enerd_, + force_vir, + shake_vir, + total_vir, + pres, + &signaller, + state_->box, + &bSumEkinhOld, + cglo_flags, + step, + &observablesReducer); + // Clean up after pre-step use of compute_globals() + observablesReducer.markAsReadyToReduce(); + } + + /* Note: this is OK, but there are some numerical precision issues with using the convergence of + the virial that should probably be addressed eventually. state->veta has better properies, + but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could + generate the new shake_vir, but test the veta value for convergence. This will take some thought. */ + + /* Output stuff */ + if (MAIN(cr_)) + { + const bool bCalcEnerStep = true; + energyOutput.addDataAtEnergyStep(doFreeEnergyPerturbation, + bCalcEnerStep, + t, + mdatoms->tmass, + enerd_, + ir->fepvals.get(), + state_->box, + PTCouplingArrays({ state_->boxv, + state_->nosehoover_xi, + state_->nosehoover_vxi, + state_->nhpres_xi, + state_->nhpres_vxi }), + state_->fep_state, + total_vir, + pres, + ekind_, + mu_tot, + constr_); + + const bool do_ene = true; + const bool do_log = true; + Awh* awh = nullptr; + const bool do_dr = ir->nstdisreout != 0; + const bool do_or = ir->nstorireout != 0; + + EnergyOutput::printAnnealingTemperatures(do_log ? fpLog_ : nullptr, *groups, ir->opts, *ekind_); + energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), + do_ene, + do_dr, + do_or, + do_log ? fpLog_ : nullptr, + step, + t, + fr_->fcdata.get(), + awh); + + if (ir->bPull) + { + pull_print_output(pullWork_, step, t); + } + + if (do_per_step(step, ir->nstlog)) + { + if (fflush(fpLog_) != 0) + { + gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?"); + } + } + } + + /* Print the remaining wall clock time for the run */ + if (isMainSimMainRank(ms_, MAIN(cr_)) && (mdrunOptions_.verbose || gmx_got_usr_signal())) + { + if (shellfc) + { + fprintf(stderr, "\n"); + } + print_time(stderr, wallTimeAccounting_, step, ir, cr_); + } + + /* Ion/water position swapping. + * Not done in last step since trajectory writing happens before this call + * in the MD loop and exchanges would be lost anyway. */ + if ((ir->eSwapCoords != SwapType::No) && (step > 0) && !isLastStep + && do_per_step(step, ir->swap->nstswap)) + { + const bool doRerun = true; + do_swapcoords(cr_, + step, + t, + ir, + swap_, + wallCycleCounters_, + rerun_fr.x, + rerun_fr.box, + MAIN(cr_) && mdrunOptions_.verbose, + doRerun); + } + + 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); + } + + if (PAR(cr_)) + { + rerun_parallel_comm(cr_, &rerun_fr, &isLastStep); + } + + cycles = wallcycle_stop(wallCycleCounters_, WallCycleCounter::Step); + if (haveDDAtomOrdering(*cr_) && wallCycleCounters_) + { + dd_cycles_add(cr_->dd, cycles, ddCyclStep); + } + + if (!rerun_fr.bStep) + { + /* increase the MD step number */ + step++; + step_rel++; + } + observablesReducer.markAsReadyToReduce(); + } + /* End of main MD loop */ + + /* Closing TNG files can include compressing data. Therefore it is good to do that + * before stopping the time measurements. */ + mdoutf_tng_close(outf); + + /* Stop measuring walltime */ + walltime_accounting_end_time(wallTimeAccounting_); + + if (MAIN(cr_)) + { + close_trx(status); + } + + if (!thisRankHasDuty(cr_, DUTY_PME)) + { + /* Tell the PME only node to finish */ + gmx_pme_send_finish(cr_); + } + + done_mdoutf(outf); + + done_shellfc(fpLog_, shellfc, step_rel); + + walltime_accounting_set_nsteps_done(wallTimeAccounting_, step_rel); +} diff --git a/patches/gromacs-2025.0.diff/src/gromacs/mdrun/rerun.cpp.preplumed b/patches/gromacs-2025.0.diff/src/gromacs/mdrun/rerun.cpp.preplumed new file mode 100644 index 0000000000..cc2e4ae4c4 --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/mdrun/rerun.cpp.preplumed @@ -0,0 +1,942 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2018- The GROMACS Authors + * and the project initiators Erik Lindahl, Berk Hess and David van der Spoel. + * Consult the AUTHORS/COPYING files and https://www.gromacs.org for details. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * https://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at https://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out https://www.gromacs.org. + */ +/*! \internal \file + * + * \brief Implements the loop for simulation reruns + * + * \author Pascal Merz + * \ingroup module_mdrun + */ +#include "gmxpre.h" + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "gromacs/applied_forces/awh/awh.h" +#include "gromacs/commandline/filenm.h" +#include "gromacs/compat/pointers.h" +#include "gromacs/domdec/collect.h" +#include "gromacs/domdec/dlbtiming.h" +#include "gromacs/domdec/domdec.h" +#include "gromacs/domdec/domdec_network.h" +#include "gromacs/domdec/domdec_struct.h" +#include "gromacs/domdec/localtopologychecker.h" +#include "gromacs/domdec/mdsetup.h" +#include "gromacs/domdec/partition.h" +#include "gromacs/essentialdynamics/edsam.h" +#include "gromacs/ewald/pme_load_balancing.h" +#include "gromacs/ewald/pme_pp.h" +#include "gromacs/fileio/enxio.h" +#include "gromacs/fileio/trxio.h" +#include "gromacs/gmxlib/network.h" +#include "gromacs/gmxlib/nrnb.h" +#include "gromacs/gpu_utils/gpu_utils.h" +#include "gromacs/listed_forces/listed_forces.h" +#include "gromacs/listed_forces/listed_forces_gpu.h" +#include "gromacs/math/functions.h" +#include "gromacs/math/utilities.h" +#include "gromacs/math/vec.h" +#include "gromacs/math/vectypes.h" +#include "gromacs/mdlib/checkpointhandler.h" +#include "gromacs/mdlib/constr.h" +#include "gromacs/mdlib/ebin.h" +#include "gromacs/mdlib/enerdata_utils.h" +#include "gromacs/mdlib/energyoutput.h" +#include "gromacs/mdlib/expanded.h" +#include "gromacs/mdlib/force.h" +#include "gromacs/mdlib/force_flags.h" +#include "gromacs/mdlib/forcerec.h" +#include "gromacs/mdlib/freeenergyparameters.h" +#include "gromacs/mdlib/md_support.h" +#include "gromacs/mdlib/mdatoms.h" +#include "gromacs/mdlib/mdoutf.h" +#include "gromacs/mdlib/membed.h" +#include "gromacs/mdlib/resethandler.h" +#include "gromacs/mdlib/sighandler.h" +#include "gromacs/mdlib/simulationsignal.h" +#include "gromacs/mdlib/stat.h" +#include "gromacs/mdlib/stophandler.h" +#include "gromacs/mdlib/tgroup.h" +#include "gromacs/mdlib/trajectory_writing.h" +#include "gromacs/mdlib/update.h" +#include "gromacs/mdlib/vcm.h" +#include "gromacs/mdlib/vsite.h" +#include "gromacs/mdrunutility/handlerestart.h" +#include "gromacs/mdrunutility/multisim.h" +#include "gromacs/mdrunutility/printtime.h" +#include "gromacs/mdtypes/awh_history.h" +#include "gromacs/mdtypes/awh_params.h" +#include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/df_history.h" +#include "gromacs/mdtypes/energyhistory.h" +#include "gromacs/mdtypes/forcebuffers.h" +#include "gromacs/mdtypes/forcerec.h" +#include "gromacs/mdtypes/group.h" +#include "gromacs/mdtypes/inputrec.h" +#include "gromacs/mdtypes/interaction_const.h" +#include "gromacs/mdtypes/md_enums.h" +#include "gromacs/mdtypes/mdatom.h" +#include "gromacs/mdtypes/mdrunoptions.h" +#include "gromacs/mdtypes/multipletimestepping.h" +#include "gromacs/mdtypes/observableshistory.h" +#include "gromacs/mdtypes/observablesreducer.h" +#include "gromacs/mdtypes/simulation_workload.h" +#include "gromacs/mdtypes/state.h" +#include "gromacs/mimic/utilities.h" +#include "gromacs/pbcutil/pbc.h" +#include "gromacs/pulling/output.h" +#include "gromacs/pulling/pull.h" +#include "gromacs/swap/swapcoords.h" +#include "gromacs/taskassignment/include/gromacs/taskassignment/decidesimulationworkload.h" +#include "gromacs/timing/wallcycle.h" +#include "gromacs/timing/walltime_accounting.h" +#include "gromacs/topology/atoms.h" +#include "gromacs/topology/idef.h" +#include "gromacs/topology/ifunc.h" +#include "gromacs/topology/mtop_util.h" +#include "gromacs/topology/topology.h" +#include "gromacs/trajectory/trajectoryframe.h" +#include "gromacs/utility/arrayref.h" +#include "gromacs/utility/basedefinitions.h" +#include "gromacs/utility/cstringutil.h" +#include "gromacs/utility/enumerationhelpers.h" +#include "gromacs/utility/exceptions.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/gmxassert.h" +#include "gromacs/utility/logger.h" +#include "gromacs/utility/real.h" + +#include "legacysimulator.h" +#include "replicaexchange.h" +#include "shellfc.h" + +struct gmx_edsam; +struct gmx_mdoutf; +struct gmx_shellfc_t; + +using gmx::SimulationSignaller; +using gmx::VirtualSitesHandler; + +/*! \brief Copy the state from \p rerunFrame to \p globalState and, if requested, construct vsites + * + * \param[in] rerunFrame The trajectory frame to compute energy/forces for + * \param[in,out] globalState The global state container + * \param[in] constructVsites When true, vsite coordinates are constructed + * \param[in] vsite Vsite setup, can be nullptr when \p constructVsites = false + */ +static void prepareRerunState(const t_trxframe& rerunFrame, + t_state* globalState, + bool constructVsites, + const VirtualSitesHandler* vsite) +{ + auto x = makeArrayRef(globalState->x); + auto rerunX = arrayRefFromArray(reinterpret_cast(rerunFrame.x), globalState->numAtoms()); + std::copy(rerunX.begin(), rerunX.end(), x.begin()); + copy_mat(rerunFrame.box, globalState->box); + + if (constructVsites) + { + GMX_ASSERT(vsite, "Need valid vsite for constructing vsites"); + + vsite->construct(globalState->x, globalState->v, globalState->box, gmx::VSiteOperation::PositionsAndVelocities); + } +} + +void gmx::LegacySimulator::do_rerun() +{ + // TODO Historically, the EM and MD "integrators" used different + // names for the t_inputrec *parameter, but these must have the + // same name, now that it's a member of a struct. We use this ir + // alias to avoid a large ripple of nearly useless changes. + // t_inputrec is being replaced by IMdpOptionsProvider, so this + // will go away eventually. + const t_inputrec* ir = inputRec_; + double t; + bool isLastStep = false; + bool doFreeEnergyPerturbation = false; + unsigned int force_flags; + tensor force_vir, shake_vir, total_vir, pres; + t_trxstatus* status = nullptr; + rvec mu_tot; + t_trxframe rerun_fr; + ForceBuffers f; + gmx_global_stat_t gstat; + gmx_shellfc_t* shellfc; + + double cycles; + + GMX_RELEASE_ASSERT(cr_->dd == nullptr || !ddUsesUpdateGroups(*cr_->dd), + "Update groups are not supported with rerun"); + + SimulationSignals signals; + // Most global communication stages don't propagate mdrun + // signals, and will use this object to achieve that. + SimulationSignaller nullSignaller(nullptr, nullptr, nullptr, false, false); + + GMX_LOG(mdLog_.info) + .asParagraph() + .appendText( + "Note that it is planned that the command gmx mdrun -rerun will " + "be available in a different form in a future version of GROMACS, " + "e.g. gmx rerun -f."); + + if (ir->efep != FreeEnergyPerturbationType::No + && (mdAtoms_->mdatoms()->nMassPerturbed > 0 || (constr_ && constr_->havePerturbedConstraints()))) + { + gmx_fatal(FARGS, + "Perturbed masses or constraints are not supported by rerun. " + "Either make a .tpr without mass and constraint perturbation, " + "or use GROMACS 2018.4, 2018.5 or later 2018 version."); + } + if (ir->bExpanded) + { + gmx_fatal(FARGS, "Expanded ensemble not supported by rerun."); + } + if (ir->bSimTemp) + { + gmx_fatal(FARGS, "Simulated tempering not supported by rerun."); + } + if (ir->bDoAwh) + { + gmx_fatal(FARGS, "AWH not supported by rerun."); + } + if (replExParams_.exchangeInterval > 0) + { + gmx_fatal(FARGS, "Replica exchange not supported by rerun."); + } + if (opt2bSet("-ei", nFile_, fnm_) || observablesHistory_->edsamHistory != nullptr) + { + gmx_fatal(FARGS, "Essential dynamics not supported by rerun."); + } + if (ir->bIMD) + { + gmx_fatal(FARGS, "Interactive MD not supported by rerun."); + } + if (isMultiSim(ms_)) + { + gmx_fatal(FARGS, "Multiple simulations not supported by rerun."); + } + if (std::any_of(ir->opts.annealing, + ir->opts.annealing + ir->opts.ngtc, + [](SimulatedAnnealing i) { return i != SimulatedAnnealing::No; })) + { + gmx_fatal(FARGS, "Simulated annealing not supported by rerun."); + } + + /* Rerun can't work if an output file name is the same as the input file name. + * If this is the case, the user will get an error telling them what the issue is. + */ + if (strcmp(opt2fn("-rerun", nFile_, fnm_), opt2fn("-o", nFile_, fnm_)) == 0 + || strcmp(opt2fn("-rerun", nFile_, fnm_), opt2fn("-x", nFile_, fnm_)) == 0) + { + gmx_fatal(FARGS, + "When using mdrun -rerun, the name of the input trajectory file " + "%s cannot be identical to the name of an output file (whether " + "given explicitly with -o or -x, or by default)", + opt2fn("-rerun", nFile_, fnm_)); + } + + /* Settings for rerun */ + { + // TODO: Avoid changing inputrec (#3854) + auto* nonConstInputrec = const_cast(inputRec_); + nonConstInputrec->nstlist = 1; + nonConstInputrec->nstcalcenergy = 1; + nonConstInputrec->nstxout_compressed = 0; + } + int nstglobalcomm = 1; + const bool bNS = true; + + ObservablesReducer observablesReducer = observablesReducerBuilder_->build(); + + const SimulationGroups* groups = &topGlobal_.groups; + if (ir->eI == IntegrationAlgorithm::Mimic) + { + auto* nonConstGlobalTopology = const_cast(&topGlobal_); + nonConstGlobalTopology->intermolecularExclusionGroup = genQmmmIndices(topGlobal_); + } + int* fep_state = MAIN(cr_) ? &stateGlobal_->fep_state : nullptr; + gmx::ArrayRef lambda = MAIN(cr_) ? stateGlobal_->lambda : gmx::ArrayRef(); + initialize_lambdas( + fpLog_, ir->efep, ir->bSimTemp, *ir->fepvals, ir->simtempvals->temperatures, ekind_, MAIN(cr_), fep_state, lambda); + const bool simulationsShareState = false; + gmx_mdoutf* outf = init_mdoutf(fpLog_, + nFile_, + fnm_, + mdrunOptions_, + cr_, + outputProvider_, + mdModulesNotifiers_, + ir, + topGlobal_, + oenv_, + wallCycleCounters_, + StartingBehavior::NewSimulation, + simulationsShareState, + ms_); + gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), + topGlobal_, + *ir, + pullWork_, + mdoutf_get_fp_dhdl(outf), + true, + StartingBehavior::NewSimulation, + simulationsShareState, + mdModulesNotifiers_); + + gstat = global_stat_init(ir); + + /* Check for polarizable models and flexible constraints */ + shellfc = init_shell_flexcon(fpLog_, + topGlobal_, + constr_ ? constr_->numFlexibleConstraints() : 0, + ir->nstcalcenergy, + haveDDAtomOrdering(*cr_), + runScheduleWork_->simulationWork.useGpuPme); + + if (haveDDAtomOrdering(*cr_)) + { + // Local state only becomes valid now. + dd_init_local_state(*cr_->dd, stateGlobal_, state_); + + /* Distribute the charge groups over the nodes from the main node */ + dd_partition_system(fpLog_, + mdLog_, + ir->init_step, + cr_, + TRUE, + stateGlobal_, + topGlobal_, + *ir, + mdModulesNotifiers_, + imdSession_, + pullWork_, + state_, + &f, + mdAtoms_, + top_, + fr_, + virtualSites_, + constr_, + nrnb_, + nullptr, + FALSE); + } + else + { + /* Copy the pointer to the global state */ + state_ = stateGlobal_; + + mdAlgorithmsSetupAtomData( + cr_, *ir, topGlobal_, top_, fr_, &f, mdAtoms_, constr_, virtualSites_, shellfc); + } + + auto* mdatoms = mdAtoms_->mdatoms(); + fr_->longRangeNonbondeds->updateAfterPartition(*mdatoms); + + // NOTE: The global state is no longer used at this point. + // But state_global is still used as temporary storage space for writing + // the global state to file and potentially for replica exchange. + // (Global topology should persist.) + + update_mdatoms(mdatoms, state_->lambda[FreeEnergyPerturbationCouplingType::Mass]); + + if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl != 0) + { + doFreeEnergyPerturbation = true; + } + + int64_t step = ir->init_step; + int64_t step_rel = 0; + + { + int cglo_flags = CGLO_GSTAT; + bool bSumEkinhOld = false; + t_vcm* vcm = nullptr; + compute_globals(gstat, + cr_, + ir, + fr_, + ekind_, + makeConstArrayRef(state_->x), + makeConstArrayRef(state_->v), + state_->box, + mdatoms, + nrnb_, + vcm, + nullptr, + enerd_, + force_vir, + shake_vir, + total_vir, + pres, + &nullSignaller, + state_->box, + &bSumEkinhOld, + cglo_flags, + step, + &observablesReducer); + // Clean up after pre-step use of compute_globals() + observablesReducer.markAsReadyToReduce(); + } + + if (MAIN(cr_)) + { + fprintf(stderr, + "starting md rerun '%s', reading coordinates from" + " input trajectory '%s'\n\n", + *(topGlobal_.name), + opt2fn("-rerun", nFile_, fnm_)); + if (mdrunOptions_.verbose) + { + fprintf(stderr, + "Calculated time to finish depends on nsteps from " + "run input file,\nwhich may not correspond to the time " + "needed to process input trajectory.\n\n"); + } + fprintf(fpLog_, "\n"); + } + + walltime_accounting_start_time(wallTimeAccounting_); + wallcycle_start(wallCycleCounters_, WallCycleCounter::Run); + print_start(fpLog_, cr_, wallTimeAccounting_, "mdrun"); + + /*********************************************************** + * + * Loop over MD steps + * + ************************************************************/ + + if (constr_) + { + GMX_LOG(mdLog_.info) + .asParagraph() + .appendText("Simulations has constraints. Rerun does not recalculate constraints."); + } + + rerun_fr.natoms = 0; + if (MAIN(cr_)) + { + isLastStep = !read_first_frame(oenv_, &status, opt2fn("-rerun", nFile_, fnm_), &rerun_fr, TRX_NEED_X); + if (rerun_fr.natoms != topGlobal_.natoms) + { + gmx_fatal(FARGS, + "Number of atoms in trajectory (%d) does not match the " + "run input file (%d)\n", + rerun_fr.natoms, + topGlobal_.natoms); + } + + if (ir->pbcType != PbcType::No) + { + if (!rerun_fr.bBox) + { + gmx_fatal(FARGS, + "Rerun trajectory frame step %" PRId64 + " time %f " + "does not contain a box, while pbc is used", + rerun_fr.step, + rerun_fr.time); + } + if (max_cutoff2(ir->pbcType, rerun_fr.box) < gmx::square(fr_->rlist)) + { + gmx_fatal(FARGS, + "Rerun trajectory frame step %" PRId64 + " time %f " + "has too small box dimensions", + rerun_fr.step, + rerun_fr.time); + } + } + } + + GMX_LOG(mdLog_.info) + .asParagraph() + .appendText( + "Rerun does not report kinetic energy, total energy, temperature, virial and " + "pressure."); + + if (PAR(cr_)) + { + rerun_parallel_comm(cr_, &rerun_fr, &isLastStep); + } + + if (ir->pbcType != PbcType::No) + { + /* Set the shift vectors. + * Necessary here when have a static box different from the tpr box. + */ + calc_shifts(rerun_fr.box, fr_->shift_vec); + } + + auto stopHandler = stopHandlerBuilder_->getStopHandlerMD( + compat::not_null(&signals[eglsSTOPCOND]), + false, + MAIN(cr_), + ir->nstlist, + mdrunOptions_.reproducible, + nstglobalcomm, + mdrunOptions_.maximumHoursToRun, + ir->nstlist == 0, + fpLog_, + step, + bNS, + wallTimeAccounting_); + + // we don't do counter resetting in rerun - finish will always be valid + walltime_accounting_set_valid_finish(wallTimeAccounting_); + + const DDBalanceRegionHandler ddBalanceRegionHandler(cr_); + + /* and stop now if we should */ + isLastStep = (isLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps)); + while (!isLastStep) + { + wallcycle_start(wallCycleCounters_, WallCycleCounter::Step); + + if (rerun_fr.bStep) + { + step = rerun_fr.step; + step_rel = step - ir->init_step; + } + if (rerun_fr.bTime) + { + t = rerun_fr.time; + } + else + { + t = step; + } + + if (ir->efep != FreeEnergyPerturbationType::No && MAIN(cr_)) + { + if (rerun_fr.bLambda) + { + ir->fepvals->init_lambda_without_states = rerun_fr.lambda; + } + else + { + if (rerun_fr.bFepState) + { + state_->fep_state = rerun_fr.fep_state; + } + } + + stateGlobal_->lambda = currentLambdas(step, *(ir->fepvals), state_->fep_state); + } + + if (MAIN(cr_)) + { + const bool constructVsites = + ((virtualSites_ != nullptr) && mdrunOptions_.rerunConstructVsites); + if (constructVsites && haveDDAtomOrdering(*cr_)) + { + gmx_fatal(FARGS, + "Vsite recalculation with -rerun is not implemented with domain " + "decomposition, " + "use a single rank"); + } + prepareRerunState(rerun_fr, stateGlobal_, constructVsites, virtualSites_); + } + + isLastStep = isLastStep || stopHandler->stoppingAfterCurrentStep(bNS); + + if (haveDDAtomOrdering(*cr_)) + { + /* Repartition the domain decomposition */ + const bool bMainState = true; + dd_partition_system(fpLog_, + mdLog_, + step, + cr_, + bMainState, + stateGlobal_, + topGlobal_, + *ir, + mdModulesNotifiers_, + imdSession_, + pullWork_, + state_, + &f, + mdAtoms_, + top_, + fr_, + virtualSites_, + constr_, + nrnb_, + wallCycleCounters_, + mdrunOptions_.verbose); + } + + if (MAIN(cr_)) + { + EnergyOutput::printHeader(fpLog_, step, t); /* can we improve the information printed here? */ + } + + if (ir->efep != FreeEnergyPerturbationType::No) + { + update_mdatoms(mdatoms, state_->lambda[FreeEnergyPerturbationCouplingType::Mass]); + } + + fr_->longRangeNonbondeds->updateAfterPartition(*mdatoms); + + force_flags = (GMX_FORCE_STATECHANGED | GMX_FORCE_DYNAMICBOX | GMX_FORCE_ALLFORCES + | GMX_FORCE_VIRIAL | // TODO: Get rid of this once #2649 and #3400 are solved + GMX_FORCE_ENERGY | (doFreeEnergyPerturbation ? GMX_FORCE_DHDL : 0)); + + const int shellfcFlags = force_flags | (mdrunOptions_.verbose ? GMX_FORCE_ENERGY : 0); + const int legacyForceFlags = ((shellfc) ? shellfcFlags : force_flags) | GMX_FORCE_NS; + + gmx_edsam* const ed = nullptr; + + if (bNS) + { + if (fr_->listedForcesGpu) + { + fr_->listedForcesGpu->updateHaveInteractions(top_->idef); + } + runScheduleWork_->domainWork = setupDomainLifetimeWorkload( + *ir, *fr_, pullWork_, ed, *mdatoms, runScheduleWork_->simulationWork); + } + + + runScheduleWork_->stepWork = setupStepWorkload(legacyForceFlags, + ir->mtsLevels, + step, + runScheduleWork_->domainWork, + runScheduleWork_->simulationWork); + + if (shellfc) + { + /* Now is the time to relax the shells */ + relax_shell_flexcon(fpLog_, + cr_, + ms_, + mdrunOptions_.verbose, + enforcedRotation_, + step, + ir, + mdModulesNotifiers_, + imdSession_, + pullWork_, + bNS, + top_, + constr_, + enerd_, + state_->numAtoms(), + state_->x.arrayRefWithPadding(), + state_->v.arrayRefWithPadding(), + state_->box, + state_->lambda, + &state_->hist, + &f.view(), + force_vir, + *mdatoms, + fr_->longRangeNonbondeds.get(), + nrnb_, + wallCycleCounters_, + shellfc, + fr_, + *runScheduleWork_, + t, + mu_tot, + virtualSites_, + ddBalanceRegionHandler); + } + else + { + /* The coordinates (x) are shifted (to get whole molecules) + * in do_force. + * This is parallellized as well, and does communication too. + * Check comments in sim_util.c + */ + Awh* awh = nullptr; + + try + { + do_force(fpLog_, + cr_, + ms_, + *ir, + mdModulesNotifiers_, + awh, + enforcedRotation_, + imdSession_, + pullWork_, + step, + nrnb_, + wallCycleCounters_, + top_, + state_->box, + state_->x.arrayRefWithPadding(), + {}, + &state_->hist, + &f.view(), + force_vir, + mdatoms, + enerd_, + state_->lambda, + fr_, + *runScheduleWork_, + virtualSites_, + mu_tot, + t, + ed, + fr_->longRangeNonbondeds.get(), + ddBalanceRegionHandler); + } + catch (const gmx::InternalError&) + { + GMX_LOG(mdLog_.warning) + .asParagraph() + .appendText( + "Continuing with next frame after catching invalid force in " + "previous frame"); + }; + } + + /* Now we have the energies and forces corresponding to the + * coordinates at time t. + */ + { + const bool isCheckpointingStep = false; + const bool doRerun = true; + do_md_trajectory_writing(fpLog_, + cr_, + nFile_, + fnm_, + step, + step_rel, + t, + ir, + state_, + stateGlobal_, + observablesHistory_, + topGlobal_, + fr_, + outf, + energyOutput, + ekind_, + f.view().force(), + isCheckpointingStep, + doRerun, + isLastStep, + mdrunOptions_.writeConfout, + EkindataState::NotUsed); + } + + stopHandler->setSignal(); + + { + const bool doInterSimSignal = false; + const bool doIntraSimSignal = true; + bool bSumEkinhOld = false; + t_vcm* vcm = nullptr; + SimulationSignaller signaller(&signals, cr_, ms_, doInterSimSignal, doIntraSimSignal); + + int cglo_flags = CGLO_GSTAT | CGLO_ENERGY; + compute_globals(gstat, + cr_, + ir, + fr_, + ekind_, + makeConstArrayRef(state_->x), + makeConstArrayRef(state_->v), + state_->box, + mdatoms, + nrnb_, + vcm, + wallCycleCounters_, + enerd_, + force_vir, + shake_vir, + total_vir, + pres, + &signaller, + state_->box, + &bSumEkinhOld, + cglo_flags, + step, + &observablesReducer); + // Clean up after pre-step use of compute_globals() + observablesReducer.markAsReadyToReduce(); + } + + /* Note: this is OK, but there are some numerical precision issues with using the convergence of + the virial that should probably be addressed eventually. state->veta has better properies, + but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could + generate the new shake_vir, but test the veta value for convergence. This will take some thought. */ + + /* Output stuff */ + if (MAIN(cr_)) + { + const bool bCalcEnerStep = true; + energyOutput.addDataAtEnergyStep(doFreeEnergyPerturbation, + bCalcEnerStep, + t, + mdatoms->tmass, + enerd_, + ir->fepvals.get(), + state_->box, + PTCouplingArrays({ state_->boxv, + state_->nosehoover_xi, + state_->nosehoover_vxi, + state_->nhpres_xi, + state_->nhpres_vxi }), + state_->fep_state, + total_vir, + pres, + ekind_, + mu_tot, + constr_); + + const bool do_ene = true; + const bool do_log = true; + Awh* awh = nullptr; + const bool do_dr = ir->nstdisreout != 0; + const bool do_or = ir->nstorireout != 0; + + EnergyOutput::printAnnealingTemperatures(do_log ? fpLog_ : nullptr, *groups, ir->opts, *ekind_); + energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), + do_ene, + do_dr, + do_or, + do_log ? fpLog_ : nullptr, + step, + t, + fr_->fcdata.get(), + awh); + + if (ir->bPull) + { + pull_print_output(pullWork_, step, t); + } + + if (do_per_step(step, ir->nstlog)) + { + if (fflush(fpLog_) != 0) + { + gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?"); + } + } + } + + /* Print the remaining wall clock time for the run */ + if (isMainSimMainRank(ms_, MAIN(cr_)) && (mdrunOptions_.verbose || gmx_got_usr_signal())) + { + if (shellfc) + { + fprintf(stderr, "\n"); + } + print_time(stderr, wallTimeAccounting_, step, ir, cr_); + } + + /* Ion/water position swapping. + * Not done in last step since trajectory writing happens before this call + * in the MD loop and exchanges would be lost anyway. */ + if ((ir->eSwapCoords != SwapType::No) && (step > 0) && !isLastStep + && do_per_step(step, ir->swap->nstswap)) + { + const bool doRerun = true; + do_swapcoords(cr_, + step, + t, + ir, + swap_, + wallCycleCounters_, + rerun_fr.x, + rerun_fr.box, + MAIN(cr_) && mdrunOptions_.verbose, + doRerun); + } + + if (MAIN(cr_)) + { + /* read next frame from input trajectory */ + isLastStep = !read_next_frame(oenv_, status, &rerun_fr); + } + + if (PAR(cr_)) + { + rerun_parallel_comm(cr_, &rerun_fr, &isLastStep); + } + + cycles = wallcycle_stop(wallCycleCounters_, WallCycleCounter::Step); + if (haveDDAtomOrdering(*cr_) && wallCycleCounters_) + { + dd_cycles_add(cr_->dd, cycles, ddCyclStep); + } + + if (!rerun_fr.bStep) + { + /* increase the MD step number */ + step++; + step_rel++; + } + observablesReducer.markAsReadyToReduce(); + } + /* End of main MD loop */ + + /* Closing TNG files can include compressing data. Therefore it is good to do that + * before stopping the time measurements. */ + mdoutf_tng_close(outf); + + /* Stop measuring walltime */ + walltime_accounting_end_time(wallTimeAccounting_); + + if (MAIN(cr_)) + { + close_trx(status); + } + + if (!thisRankHasDuty(cr_, DUTY_PME)) + { + /* Tell the PME only node to finish */ + gmx_pme_send_finish(cr_); + } + + done_mdoutf(outf); + + done_shellfc(fpLog_, shellfc, step_rel); + + walltime_accounting_set_nsteps_done(wallTimeAccounting_, step_rel); +} diff --git a/patches/gromacs-2026.0.diff/src/gromacs/mdrun/rerun.cpp b/patches/gromacs-2026.0.diff/src/gromacs/mdrun/rerun.cpp new file mode 100644 index 0000000000..dcdec0d180 --- /dev/null +++ b/patches/gromacs-2026.0.diff/src/gromacs/mdrun/rerun.cpp @@ -0,0 +1,942 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2018- The GROMACS Authors + * and the project initiators Erik Lindahl, Berk Hess and David van der Spoel. + * Consult the AUTHORS/COPYING files and https://www.gromacs.org for details. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * https://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at https://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out https://www.gromacs.org. + */ +/*! \internal \file + * + * \brief Implements the loop for simulation reruns + * + * \author Pascal Merz + * \ingroup module_mdrun + */ +#include "gmxpre.h" + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "gromacs/applied_forces/awh/awh.h" +#include "gromacs/commandline/filenm.h" +#include "gromacs/compat/pointers.h" +#include "gromacs/domdec/collect.h" +#include "gromacs/domdec/dlbtiming.h" +#include "gromacs/domdec/domdec.h" +#include "gromacs/domdec/domdec_network.h" +#include "gromacs/domdec/domdec_struct.h" +#include "gromacs/domdec/localtopologychecker.h" +#include "gromacs/domdec/mdsetup.h" +#include "gromacs/domdec/partition.h" +#include "gromacs/essentialdynamics/edsam.h" +#include "gromacs/ewald/pme_load_balancing.h" +#include "gromacs/ewald/pme_pp.h" +#include "gromacs/fileio/enxio.h" +#include "gromacs/fileio/trxio.h" +#include "gromacs/gmxlib/network.h" +#include "gromacs/gmxlib/nrnb.h" +#include "gromacs/gpu_utils/gpu_utils.h" +#include "gromacs/listed_forces/listed_forces.h" +#include "gromacs/listed_forces/listed_forces_gpu.h" +#include "gromacs/math/functions.h" +#include "gromacs/math/utilities.h" +#include "gromacs/mdlib/checkpointhandler.h" +#include "gromacs/mdlib/constr.h" +#include "gromacs/mdlib/ebin.h" +#include "gromacs/mdlib/enerdata_utils.h" +#include "gromacs/mdlib/energyoutput.h" +#include "gromacs/mdlib/expanded.h" +#include "gromacs/mdlib/force.h" +#include "gromacs/mdlib/force_flags.h" +#include "gromacs/mdlib/forcerec.h" +#include "gromacs/mdlib/freeenergyparameters.h" +#include "gromacs/mdlib/md_support.h" +#include "gromacs/mdlib/mdatoms.h" +#include "gromacs/mdlib/mdoutf.h" +#include "gromacs/mdlib/membed.h" +#include "gromacs/mdlib/resethandler.h" +#include "gromacs/mdlib/sighandler.h" +#include "gromacs/mdlib/simulationsignal.h" +#include "gromacs/mdlib/stat.h" +#include "gromacs/mdlib/stophandler.h" +#include "gromacs/mdlib/tgroup.h" +#include "gromacs/mdlib/trajectory_writing.h" +#include "gromacs/mdlib/update.h" +#include "gromacs/mdlib/vcm.h" +#include "gromacs/mdlib/vsite.h" +#include "gromacs/mdrunutility/handlerestart.h" +#include "gromacs/mdrunutility/multisim.h" +#include "gromacs/mdrunutility/printtime.h" +#include "gromacs/mdtypes/awh_history.h" +#include "gromacs/mdtypes/awh_params.h" +#include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/df_history.h" +#include "gromacs/mdtypes/energyhistory.h" +#include "gromacs/mdtypes/forcebuffers.h" +#include "gromacs/mdtypes/forcerec.h" +#include "gromacs/mdtypes/group.h" +#include "gromacs/mdtypes/inputrec.h" +#include "gromacs/mdtypes/interaction_const.h" +#include "gromacs/mdtypes/md_enums.h" +#include "gromacs/mdtypes/mdatom.h" +#include "gromacs/mdtypes/mdrunoptions.h" +#include "gromacs/mdtypes/multipletimestepping.h" +#include "gromacs/mdtypes/observableshistory.h" +#include "gromacs/mdtypes/observablesreducer.h" +#include "gromacs/mdtypes/simulation_workload.h" +#include "gromacs/mdtypes/state.h" +#include "gromacs/mimic/utilities.h" +#include "gromacs/pbcutil/pbc.h" +#include "gromacs/pulling/output.h" +#include "gromacs/pulling/pull.h" +#include "gromacs/swap/swapcoords.h" +#include "gromacs/taskassignment/include/gromacs/taskassignment/decidesimulationworkload.h" +#include "gromacs/timing/wallcycle.h" +#include "gromacs/timing/walltime_accounting.h" +#include "gromacs/topology/atoms.h" +#include "gromacs/topology/idef.h" +#include "gromacs/topology/ifunc.h" +#include "gromacs/topology/mtop_util.h" +#include "gromacs/topology/topology.h" +#include "gromacs/trajectory/trajectoryframe.h" +#include "gromacs/utility/arrayref.h" +#include "gromacs/utility/basedefinitions.h" +#include "gromacs/utility/cstringutil.h" +#include "gromacs/utility/enumerationhelpers.h" +#include "gromacs/utility/exceptions.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/gmxassert.h" +#include "gromacs/utility/logger.h" +#include "gromacs/utility/real.h" +#include "gromacs/utility/vec.h" +#include "gromacs/utility/vectypes.h" + +#include "legacysimulator.h" +#include "replicaexchange.h" +#include "shellfc.h" + +struct gmx_edsam; +struct gmx_mdoutf; +struct gmx_shellfc_t; + +using gmx::SimulationSignaller; +using gmx::VirtualSitesHandler; + +/*! \brief Copy the state from \p rerunFrame to \p globalState and, if requested, construct vsites + * + * \param[in] rerunFrame The trajectory frame to compute energy/forces for + * \param[in,out] globalState The global state container + * \param[in] constructVsites When true, vsite coordinates are constructed + * \param[in] vsite Vsite setup, can be nullptr when \p constructVsites = false + */ +static void prepareRerunState(const t_trxframe& rerunFrame, + t_state* globalState, + bool constructVsites, + const VirtualSitesHandler* vsite) +{ + auto x = makeArrayRef(globalState->x); + auto rerunX = arrayRefFromArray(reinterpret_cast(rerunFrame.x), globalState->numAtoms()); + std::copy(rerunX.begin(), rerunX.end(), x.begin()); + copy_mat(rerunFrame.box, globalState->box); + + if (constructVsites) + { + GMX_ASSERT(vsite, "Need valid vsite for constructing vsites"); + + vsite->construct(globalState->x, globalState->v, globalState->box, gmx::VSiteOperation::PositionsAndVelocities); + } +} + +void gmx::LegacySimulator::do_rerun() +{ + const bool isMainRank = cr_->commMyGroup.isMainRank(); + + // TODO Historically, the EM and MD "integrators" used different + // names for the t_inputrec *parameter, but these must have the + // same name, now that it's a member of a struct. We use this ir + // alias to avoid a large ripple of nearly useless changes. + // t_inputrec is being replaced by IMdpOptionsProvider, so this + // will go away eventually. + const t_inputrec* ir = inputRec_; + double t; + bool isLastStep = false; + bool doFreeEnergyPerturbation = false; + unsigned int force_flags; + tensor force_vir, shake_vir, total_vir, pres; + t_trxstatus* status = nullptr; + rvec mu_tot; + t_trxframe rerun_fr; + ForceBuffers f; + gmx_global_stat_t gstat; + gmx_shellfc_t* shellfc; + + double cycles; + + GMX_RELEASE_ASSERT(cr_->dd == nullptr || !ddUsesUpdateGroups(*cr_->dd), + "Update groups are not supported with rerun"); + + SimulationSignals signals; + // Most global communication stages don't propagate mdrun + // signals, and will use this object to achieve that. + SimulationSignaller nullSignaller(nullptr, nullptr, nullptr, false, false); + + GMX_LOG(mdLog_.info) + .asParagraph() + .appendText( + "Note that it is planned that the command gmx mdrun -rerun will " + "be available in a different form in a future version of GROMACS, " + "e.g. gmx rerun -f."); + + if (ir->efep != FreeEnergyPerturbationType::No + && (mdAtoms_->mdatoms()->nMassPerturbed > 0 || (constr_ && constr_->havePerturbedConstraints()))) + { + gmx_fatal(FARGS, + "Perturbed masses or constraints are not supported by rerun. " + "Either make a .tpr without mass and constraint perturbation, " + "or use GROMACS 2018.4, 2018.5 or later 2018 version."); + } + if (ir->bExpanded) + { + gmx_fatal(FARGS, "Expanded ensemble not supported by rerun."); + } + if (ir->bSimTemp) + { + gmx_fatal(FARGS, "Simulated tempering not supported by rerun."); + } + if (ir->bDoAwh) + { + gmx_fatal(FARGS, "AWH not supported by rerun."); + } + if (replExParams_.exchangeInterval > 0) + { + gmx_fatal(FARGS, "Replica exchange not supported by rerun."); + } + if (opt2bSet("-ei", nFile_, fnm_) || observablesHistory_->edsamHistory != nullptr) + { + gmx_fatal(FARGS, "Essential dynamics not supported by rerun."); + } + if (ir->bIMD) + { + gmx_fatal(FARGS, "Interactive MD not supported by rerun."); + } + if (isMultiSim(ms_)) + { + gmx_fatal(FARGS, "Multiple simulations not supported by rerun."); + } + if (std::any_of(ir->opts.annealing, + ir->opts.annealing + ir->opts.ngtc, + [](SimulatedAnnealing i) { return i != SimulatedAnnealing::No; })) + { + gmx_fatal(FARGS, "Simulated annealing not supported by rerun."); + } + + /* Rerun can't work if an output file name is the same as the input file name. + * If this is the case, the user will get an error telling them what the issue is. + */ + if (std::strcmp(opt2fn("-rerun", nFile_, fnm_), opt2fn("-o", nFile_, fnm_)) == 0 + || std::strcmp(opt2fn("-rerun", nFile_, fnm_), opt2fn("-x", nFile_, fnm_)) == 0) + { + gmx_fatal(FARGS, + "When using mdrun -rerun, the name of the input trajectory file " + "%s cannot be identical to the name of an output file (whether " + "given explicitly with -o or -x, or by default)", + opt2fn("-rerun", nFile_, fnm_)); + } + + /* Settings for rerun */ + { + // TODO: Avoid changing inputrec (#3854) + auto* nonConstInputrec = const_cast(inputRec_); + nonConstInputrec->nstlist = 1; + nonConstInputrec->nstcalcenergy = 1; + nonConstInputrec->nstxout_compressed = 0; + } + int nstglobalcomm = 1; + const bool bNS = true; + + ObservablesReducer observablesReducer = observablesReducerBuilder_->build(); + + const SimulationGroups* groups = &topGlobal_.groups; + if (ir->eI == IntegrationAlgorithm::Mimic) + { + auto* nonConstGlobalTopology = const_cast(&topGlobal_); + nonConstGlobalTopology->intermolecularExclusionGroup = genQmmmIndices(topGlobal_); + } + int* fep_state = isMainRank ? &stateGlobal_->fep_state : nullptr; + gmx::ArrayRef lambda = isMainRank ? stateGlobal_->lambda : gmx::ArrayRef{}; + initialize_lambdas( + fpLog_, ir->efep, ir->bSimTemp, *ir->fepvals, ir->simtempvals->temperatures, ekind_, isMainRank, fep_state, lambda); + const bool simulationsShareState = false; + gmx_mdoutf* outf = init_mdoutf(fpLog_, + nFile_, + fnm_, + mdrunOptions_, + cr_, + outputProvider_, + mdModulesNotifiers_, + ir, + topGlobal_, + oenv_, + wallCycleCounters_, + StartingBehavior::NewSimulation, + simulationsShareState, + ms_); + gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), + topGlobal_, + *ir, + pullWork_, + mdoutf_get_fp_dhdl(outf), + true, + StartingBehavior::NewSimulation, + simulationsShareState, + mdModulesNotifiers_); + + gstat = global_stat_init(ir); + + /* Check for polarizable models and flexible constraints */ + shellfc = init_shell_flexcon(fpLog_, + topGlobal_, + constr_ ? constr_->numFlexibleConstraints() : 0, + ir->nstcalcenergy, + haveDDAtomOrdering(*cr_), + runScheduleWork_->simulationWork); + + if (haveDDAtomOrdering(*cr_)) + { + // Local state only becomes valid now. + dd_init_local_state(*cr_->dd, stateGlobal_, state_); + + /* Distribute the charge groups over the nodes from the main node */ + dd_partition_system(fpLog_, + mdLog_, + ir->init_step, + cr_->dd, + TRUE, + stateGlobal_, + topGlobal_, + *ir, + mdModulesNotifiers_, + imdSession_, + pullWork_, + state_, + &f, + mdAtoms_, + top_, + fr_, + virtualSites_, + constr_, + nrnb_, + nullptr, + FALSE); + } + else + { + /* Copy the pointer to the global state */ + state_ = stateGlobal_; + + mdAlgorithmsSetupAtomData( + cr_->dd, *ir, topGlobal_, top_, fr_, &f, mdAtoms_, constr_, virtualSites_, shellfc); + } + + auto* mdatoms = mdAtoms_->mdatoms(); + fr_->longRangeNonbondeds->updateAfterPartition(*mdatoms); + + // NOTE: The global state is no longer used at this point. + // But state_global is still used as temporary storage space for writing + // the global state to file and potentially for replica exchange. + // (Global topology should persist.) + + update_mdatoms(mdatoms, state_->lambda[FreeEnergyPerturbationCouplingType::Mass]); + + if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl != 0) + { + doFreeEnergyPerturbation = true; + } + + int64_t step = ir->init_step; + int64_t step_rel = 0; + + { + int cglo_flags = CGLO_GSTAT; + bool bSumEkinhOld = false; + t_vcm* vcm = nullptr; + compute_globals(gstat, + cr_->commMyGroup, + ir, + fr_, + ekind_, + makeConstArrayRef(state_->x), + makeConstArrayRef(state_->v), + state_->box, + mdatoms, + nrnb_, + vcm, + nullptr, + enerd_, + force_vir, + shake_vir, + total_vir, + pres, + &nullSignaller, + state_->box, + &bSumEkinhOld, + cglo_flags, + step, + &observablesReducer); + // Clean up after pre-step use of compute_globals() + observablesReducer.markAsReadyToReduce(); + } + + if (isMainRank) + { + fprintf(stderr, + "starting md rerun '%s', reading coordinates from" + " input trajectory '%s'\n\n", + *(topGlobal_.name), + opt2fn("-rerun", nFile_, fnm_)); + if (mdrunOptions_.verbose) + { + fprintf(stderr, + "Calculated time to finish depends on nsteps from " + "run input file,\nwhich may not correspond to the time " + "needed to process input trajectory.\n\n"); + } + fprintf(fpLog_, "\n"); + } + + walltime_accounting_start_time(wallTimeAccounting_); + wallcycle_start(wallCycleCounters_, WallCycleCounter::Run); + print_start(fpLog_, cr_, wallTimeAccounting_, "mdrun"); + + /*********************************************************** + * + * Loop over MD steps + * + ************************************************************/ + + if (constr_) + { + GMX_LOG(mdLog_.info) + .asParagraph() + .appendText("Simulations has constraints. Rerun does not recalculate constraints."); + } + + rerun_fr.natoms = 0; + if (isMainRank) + { + isLastStep = !read_first_frame(oenv_, &status, opt2fn("-rerun", nFile_, fnm_), &rerun_fr, TRX_NEED_X); + if (rerun_fr.natoms != topGlobal_.natoms) + { + gmx_fatal(FARGS, + "Number of atoms in trajectory (%d) does not match the " + "run input file (%d)\n", + rerun_fr.natoms, + topGlobal_.natoms); + } + + if (ir->pbcType != PbcType::No) + { + if (!rerun_fr.bBox) + { + gmx_fatal(FARGS, + "Rerun trajectory frame step %" PRId64 + " time %f " + "does not contain a box, while pbc is used", + rerun_fr.step, + rerun_fr.time); + } + if (max_cutoff2(ir->pbcType, rerun_fr.box) < gmx::square(fr_->rlist)) + { + gmx_fatal(FARGS, + "Rerun trajectory frame step %" PRId64 + " time %f " + "has too small box dimensions", + rerun_fr.step, + rerun_fr.time); + } + } + } + + GMX_LOG(mdLog_.info) + .asParagraph() + .appendText( + "Rerun does not report kinetic energy, total energy, temperature, virial and " + "pressure."); + + if (cr_->commMyGroup.isParallel()) + { + rerun_parallel_comm(cr_->commMyGroup, &rerun_fr, &isLastStep); + } + + if (ir->pbcType != PbcType::No) + { + /* Set the shift vectors. + * Necessary here when have a static box different from the tpr box. + */ + calc_shifts(rerun_fr.box, fr_->shift_vec); + } + + auto stopHandler = stopHandlerBuilder_->getStopHandlerMD( + compat::not_null(&signals[eglsSTOPCOND]), + false, + isMainRank, + 1, // rerun constructs the pairlist for each frame + mdrunOptions_.reproducible, + nstglobalcomm, + mdrunOptions_.maximumHoursToRun, + fpLog_, + step, + wallTimeAccounting_); + + // we don't do counter resetting in rerun - finish will always be valid + walltime_accounting_set_valid_finish(wallTimeAccounting_); + + const DDBalanceRegionHandler ddBalanceRegionHandler(cr_->dd); + + /* and stop now if we should */ + isLastStep = (isLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps)); + while (!isLastStep) + { + wallcycle_start(wallCycleCounters_, WallCycleCounter::Step); + + if (rerun_fr.bStep) + { + step = rerun_fr.step; + step_rel = step - ir->init_step; + } + if (rerun_fr.bTime) + { + t = rerun_fr.time; + } + else + { + t = step; + } + + if (ir->efep != FreeEnergyPerturbationType::No && isMainRank) + { + if (rerun_fr.bLambda) + { + ir->fepvals->init_lambda_without_states = rerun_fr.lambda; + } + else + { + if (rerun_fr.bFepState) + { + state_->fep_state = rerun_fr.fep_state; + } + } + + stateGlobal_->lambda = currentLambdas(step, *(ir->fepvals), state_->fep_state); + } + + if (isMainRank) + { + const bool constructVsites = + ((virtualSites_ != nullptr) && mdrunOptions_.rerunConstructVsites); + if (constructVsites && haveDDAtomOrdering(*cr_)) + { + gmx_fatal(FARGS, + "Vsite recalculation with -rerun is not implemented with domain " + "decomposition, " + "use a single rank"); + } + prepareRerunState(rerun_fr, stateGlobal_, constructVsites, virtualSites_); + } + + isLastStep = isLastStep || stopHandler->stoppingAfterCurrentStep(step); + + if (haveDDAtomOrdering(*cr_)) + { + /* Repartition the domain decomposition */ + const bool bMainState = true; + dd_partition_system(fpLog_, + mdLog_, + step, + cr_->dd, + bMainState, + stateGlobal_, + topGlobal_, + *ir, + mdModulesNotifiers_, + imdSession_, + pullWork_, + state_, + &f, + mdAtoms_, + top_, + fr_, + virtualSites_, + constr_, + nrnb_, + wallCycleCounters_, + mdrunOptions_.verbose); + } + + if (isMainRank) + { + EnergyOutput::printHeader(fpLog_, step, t); /* can we improve the information printed here? */ + } + + if (ir->efep != FreeEnergyPerturbationType::No) + { + update_mdatoms(mdatoms, state_->lambda[FreeEnergyPerturbationCouplingType::Mass]); + } + + fr_->longRangeNonbondeds->updateAfterPartition(*mdatoms); + + force_flags = (GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES + | GMX_FORCE_VIRIAL + // TODO: Get rid of this once #2649 and #3400 are solved + | GMX_FORCE_ENERGY | (doFreeEnergyPerturbation ? GMX_FORCE_DHDL : 0)); + + const int shellfcFlags = force_flags | (mdrunOptions_.verbose ? GMX_FORCE_ENERGY : 0); + const int legacyForceFlags = ((shellfc) ? shellfcFlags : force_flags) | GMX_FORCE_NS; + + gmx_edsam* const ed = nullptr; + + if (bNS) + { + if (fr_->listedForcesGpu) + { + fr_->listedForcesGpu->updateHaveInteractions(top_->idef); + } + runScheduleWork_->domainWork = setupDomainLifetimeWorkload( + *ir, *fr_, pullWork_, ed, *mdatoms, runScheduleWork_->simulationWork); + } + + + runScheduleWork_->stepWork = setupStepWorkload(legacyForceFlags, + ir->mtsLevels, + step, + runScheduleWork_->domainWork, + runScheduleWork_->simulationWork); + + if (shellfc) + { + /* Now is the time to relax the shells */ + relax_shell_flexcon(fpLog_, + cr_, + mdrunOptions_.verbose, + enforcedRotation_, + step, + ir, + mdModulesNotifiers_, + imdSession_, + pullWork_, + bNS, + top_, + constr_, + enerd_, + state_->numAtoms(), + state_->x.arrayRefWithPadding(), + state_->v.arrayRefWithPadding(), + state_->box, + state_->lambda, + &state_->hist, + &f.view(), + force_vir, + *mdatoms, + fr_->longRangeNonbondeds.get(), + nrnb_, + wallCycleCounters_, + shellfc, + fr_, + *runScheduleWork_, + t, + mu_tot, + virtualSites_, + ddBalanceRegionHandler); + } + else + { + /* The coordinates (x) are shifted (to get whole molecules) + * in do_force. + * This is parallellized as well, and does communication too. + * Check comments in sim_util.c + */ + Awh* awh = nullptr; + + try + { + do_force(fpLog_, + cr_, + *ir, + mdModulesNotifiers_, + awh, + enforcedRotation_, + imdSession_, + pullWork_, + step, + nrnb_, + wallCycleCounters_, + top_, + state_->box, + state_->x.arrayRefWithPadding(), + {}, + &state_->hist, + &f.view(), + force_vir, + mdatoms, + enerd_, + state_->lambda, + fr_, + *runScheduleWork_, + virtualSites_, + mu_tot, + t, + ed, + fr_->longRangeNonbondeds.get(), + ddBalanceRegionHandler); + } + catch (const gmx::InternalError&) + { + GMX_LOG(mdLog_.warning) + .asParagraph() + .appendText( + "Continuing with next frame after catching invalid force in " + "previous frame"); + }; + } + + /* Now we have the energies and forces corresponding to the + * coordinates at time t. + */ + { + const bool isCheckpointingStep = false; + const bool doRerun = true; + do_md_trajectory_writing(fpLog_, + cr_, + nFile_, + fnm_, + step, + step_rel, + t, + ir, + state_, + stateGlobal_, + observablesHistory_, + topGlobal_, + fr_, + outf, + energyOutput, + ekind_, + f.view().force(), + isCheckpointingStep, + doRerun, + isLastStep, + mdrunOptions_.writeConfout, + EkindataState::NotUsed); + } + + stopHandler->setSignal(); + + { + const bool doInterSimSignal = false; + const bool doIntraSimSignal = true; + bool bSumEkinhOld = false; + t_vcm* vcm = nullptr; + SimulationSignaller signaller(&signals, cr_, ms_, doInterSimSignal, doIntraSimSignal); + + int cglo_flags = CGLO_GSTAT | CGLO_ENERGY; + compute_globals(gstat, + cr_->commMyGroup, + ir, + fr_, + ekind_, + makeConstArrayRef(state_->x), + makeConstArrayRef(state_->v), + state_->box, + mdatoms, + nrnb_, + vcm, + wallCycleCounters_, + enerd_, + force_vir, + shake_vir, + total_vir, + pres, + &signaller, + state_->box, + &bSumEkinhOld, + cglo_flags, + step, + &observablesReducer); + // Clean up after pre-step use of compute_globals() + observablesReducer.markAsReadyToReduce(); + } + + /* Note: this is OK, but there are some numerical precision issues with using the convergence of + the virial that should probably be addressed eventually. state->veta has better properties, + but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could + generate the new shake_vir, but test the veta value for convergence. This will take some thought. */ + + /* Output stuff */ + if (isMainRank) + { + const bool bCalcEnerStep = true; + energyOutput.addDataAtEnergyStep(doFreeEnergyPerturbation, + bCalcEnerStep, + t, + mdatoms->tmass, + enerd_, + ir->fepvals.get(), + state_->box, + PTCouplingArrays({ state_->boxv, + state_->nosehoover_xi, + state_->nosehoover_vxi, + state_->nhpres_xi, + state_->nhpres_vxi }), + state_->fep_state, + total_vir, + pres, + ekind_, + mu_tot, + constr_); + + const bool do_ene = true; + const bool do_log = true; + Awh* awh = nullptr; + const bool do_dr = ir->nstdisreout != 0; + const bool do_or = ir->nstorireout != 0; + + EnergyOutput::printAnnealingTemperatures(do_log ? fpLog_ : nullptr, *groups, ir->opts, *ekind_); + energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), + do_ene, + do_dr, + do_or, + do_log ? fpLog_ : nullptr, + step, + t, + fr_->fcdata.get(), + awh); + + if (ir->bPull) + { + pull_print_output(pullWork_, step, t); + } + + if (do_per_step(step, ir->nstlog)) + { + if (std::fflush(fpLog_) != 0) + { + gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?"); + } + } + } + + /* Print the remaining wall clock time for the run */ + if (isMainSimMainRank(ms_, isMainRank) && (mdrunOptions_.verbose || gmx_got_usr_signal())) + { + if (shellfc) + { + fprintf(stderr, "\n"); + } + print_time(stderr, wallTimeAccounting_, step, ir, cr_->commMySim); + } + + /* Ion/water position swapping. + * Not done in last step since trajectory writing happens before this call + * in the MD loop and exchanges would be lost anyway. */ + if ((ir->eSwapCoords != SwapType::No) && (step > 0) && !isLastStep + && do_per_step(step, ir->swap->nstswap)) + { + const bool doRerun = true; + do_swapcoords(cr_->commMyGroup, + step, + t, + ir, + swap_, + wallCycleCounters_, + gmx::arrayRefFromArray(reinterpret_cast(rerun_fr.x), rerun_fr.natoms), + rerun_fr.box, + isMainRank && mdrunOptions_.verbose, + doRerun); + } + + 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); + } + + if (cr_->commMyGroup.isParallel()) + { + rerun_parallel_comm(cr_->commMyGroup, &rerun_fr, &isLastStep); + } + + cycles = wallcycle_stop(wallCycleCounters_, WallCycleCounter::Step); + if (haveDDAtomOrdering(*cr_) && wallCycleCounters_) + { + dd_cycles_add(cr_->dd, cycles, ddCyclStep); + } + + if (!rerun_fr.bStep) + { + /* increase the MD step number */ + step++; + step_rel++; + } + observablesReducer.markAsReadyToReduce(); + } + /* End of main MD loop */ + + /* Closing TNG files can include compressing data. Therefore it is good to do that + * before stopping the time measurements. */ + mdoutf_tng_close(outf); + + /* Stop measuring walltime */ + walltime_accounting_end_time(wallTimeAccounting_); + + if (isMainRank) + { + close_trx(status); + } + + if (!thisRankHasPmeDuty(cr_->dd)) + { + /* Tell the PME only node to finish */ + gmx_pme_send_finish(cr_->dd); + } + + done_mdoutf(outf); + + done_shellfc(fpLog_, shellfc, step_rel); + + walltime_accounting_set_nsteps_done(wallTimeAccounting_, step_rel); +} diff --git a/patches/gromacs-2026.0.diff/src/gromacs/mdrun/rerun.cpp.preplumed b/patches/gromacs-2026.0.diff/src/gromacs/mdrun/rerun.cpp.preplumed new file mode 100644 index 0000000000..e72dfcf464 --- /dev/null +++ b/patches/gromacs-2026.0.diff/src/gromacs/mdrun/rerun.cpp.preplumed @@ -0,0 +1,941 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2018- The GROMACS Authors + * and the project initiators Erik Lindahl, Berk Hess and David van der Spoel. + * Consult the AUTHORS/COPYING files and https://www.gromacs.org for details. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * https://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at https://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out https://www.gromacs.org. + */ +/*! \internal \file + * + * \brief Implements the loop for simulation reruns + * + * \author Pascal Merz + * \ingroup module_mdrun + */ +#include "gmxpre.h" + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "gromacs/applied_forces/awh/awh.h" +#include "gromacs/commandline/filenm.h" +#include "gromacs/compat/pointers.h" +#include "gromacs/domdec/collect.h" +#include "gromacs/domdec/dlbtiming.h" +#include "gromacs/domdec/domdec.h" +#include "gromacs/domdec/domdec_network.h" +#include "gromacs/domdec/domdec_struct.h" +#include "gromacs/domdec/localtopologychecker.h" +#include "gromacs/domdec/mdsetup.h" +#include "gromacs/domdec/partition.h" +#include "gromacs/essentialdynamics/edsam.h" +#include "gromacs/ewald/pme_load_balancing.h" +#include "gromacs/ewald/pme_pp.h" +#include "gromacs/fileio/enxio.h" +#include "gromacs/fileio/trxio.h" +#include "gromacs/gmxlib/network.h" +#include "gromacs/gmxlib/nrnb.h" +#include "gromacs/gpu_utils/gpu_utils.h" +#include "gromacs/listed_forces/listed_forces.h" +#include "gromacs/listed_forces/listed_forces_gpu.h" +#include "gromacs/math/functions.h" +#include "gromacs/math/utilities.h" +#include "gromacs/mdlib/checkpointhandler.h" +#include "gromacs/mdlib/constr.h" +#include "gromacs/mdlib/ebin.h" +#include "gromacs/mdlib/enerdata_utils.h" +#include "gromacs/mdlib/energyoutput.h" +#include "gromacs/mdlib/expanded.h" +#include "gromacs/mdlib/force.h" +#include "gromacs/mdlib/force_flags.h" +#include "gromacs/mdlib/forcerec.h" +#include "gromacs/mdlib/freeenergyparameters.h" +#include "gromacs/mdlib/md_support.h" +#include "gromacs/mdlib/mdatoms.h" +#include "gromacs/mdlib/mdoutf.h" +#include "gromacs/mdlib/membed.h" +#include "gromacs/mdlib/resethandler.h" +#include "gromacs/mdlib/sighandler.h" +#include "gromacs/mdlib/simulationsignal.h" +#include "gromacs/mdlib/stat.h" +#include "gromacs/mdlib/stophandler.h" +#include "gromacs/mdlib/tgroup.h" +#include "gromacs/mdlib/trajectory_writing.h" +#include "gromacs/mdlib/update.h" +#include "gromacs/mdlib/vcm.h" +#include "gromacs/mdlib/vsite.h" +#include "gromacs/mdrunutility/handlerestart.h" +#include "gromacs/mdrunutility/multisim.h" +#include "gromacs/mdrunutility/printtime.h" +#include "gromacs/mdtypes/awh_history.h" +#include "gromacs/mdtypes/awh_params.h" +#include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/df_history.h" +#include "gromacs/mdtypes/energyhistory.h" +#include "gromacs/mdtypes/forcebuffers.h" +#include "gromacs/mdtypes/forcerec.h" +#include "gromacs/mdtypes/group.h" +#include "gromacs/mdtypes/inputrec.h" +#include "gromacs/mdtypes/interaction_const.h" +#include "gromacs/mdtypes/md_enums.h" +#include "gromacs/mdtypes/mdatom.h" +#include "gromacs/mdtypes/mdrunoptions.h" +#include "gromacs/mdtypes/multipletimestepping.h" +#include "gromacs/mdtypes/observableshistory.h" +#include "gromacs/mdtypes/observablesreducer.h" +#include "gromacs/mdtypes/simulation_workload.h" +#include "gromacs/mdtypes/state.h" +#include "gromacs/mimic/utilities.h" +#include "gromacs/pbcutil/pbc.h" +#include "gromacs/pulling/output.h" +#include "gromacs/pulling/pull.h" +#include "gromacs/swap/swapcoords.h" +#include "gromacs/taskassignment/include/gromacs/taskassignment/decidesimulationworkload.h" +#include "gromacs/timing/wallcycle.h" +#include "gromacs/timing/walltime_accounting.h" +#include "gromacs/topology/atoms.h" +#include "gromacs/topology/idef.h" +#include "gromacs/topology/ifunc.h" +#include "gromacs/topology/mtop_util.h" +#include "gromacs/topology/topology.h" +#include "gromacs/trajectory/trajectoryframe.h" +#include "gromacs/utility/arrayref.h" +#include "gromacs/utility/basedefinitions.h" +#include "gromacs/utility/cstringutil.h" +#include "gromacs/utility/enumerationhelpers.h" +#include "gromacs/utility/exceptions.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/gmxassert.h" +#include "gromacs/utility/logger.h" +#include "gromacs/utility/real.h" +#include "gromacs/utility/vec.h" +#include "gromacs/utility/vectypes.h" + +#include "legacysimulator.h" +#include "replicaexchange.h" +#include "shellfc.h" + +struct gmx_edsam; +struct gmx_mdoutf; +struct gmx_shellfc_t; + +using gmx::SimulationSignaller; +using gmx::VirtualSitesHandler; + +/*! \brief Copy the state from \p rerunFrame to \p globalState and, if requested, construct vsites + * + * \param[in] rerunFrame The trajectory frame to compute energy/forces for + * \param[in,out] globalState The global state container + * \param[in] constructVsites When true, vsite coordinates are constructed + * \param[in] vsite Vsite setup, can be nullptr when \p constructVsites = false + */ +static void prepareRerunState(const t_trxframe& rerunFrame, + t_state* globalState, + bool constructVsites, + const VirtualSitesHandler* vsite) +{ + auto x = makeArrayRef(globalState->x); + auto rerunX = arrayRefFromArray(reinterpret_cast(rerunFrame.x), globalState->numAtoms()); + std::copy(rerunX.begin(), rerunX.end(), x.begin()); + copy_mat(rerunFrame.box, globalState->box); + + if (constructVsites) + { + GMX_ASSERT(vsite, "Need valid vsite for constructing vsites"); + + vsite->construct(globalState->x, globalState->v, globalState->box, gmx::VSiteOperation::PositionsAndVelocities); + } +} + +void gmx::LegacySimulator::do_rerun() +{ + const bool isMainRank = cr_->commMyGroup.isMainRank(); + + // TODO Historically, the EM and MD "integrators" used different + // names for the t_inputrec *parameter, but these must have the + // same name, now that it's a member of a struct. We use this ir + // alias to avoid a large ripple of nearly useless changes. + // t_inputrec is being replaced by IMdpOptionsProvider, so this + // will go away eventually. + const t_inputrec* ir = inputRec_; + double t; + bool isLastStep = false; + bool doFreeEnergyPerturbation = false; + unsigned int force_flags; + tensor force_vir, shake_vir, total_vir, pres; + t_trxstatus* status = nullptr; + rvec mu_tot; + t_trxframe rerun_fr; + ForceBuffers f; + gmx_global_stat_t gstat; + gmx_shellfc_t* shellfc; + + double cycles; + + GMX_RELEASE_ASSERT(cr_->dd == nullptr || !ddUsesUpdateGroups(*cr_->dd), + "Update groups are not supported with rerun"); + + SimulationSignals signals; + // Most global communication stages don't propagate mdrun + // signals, and will use this object to achieve that. + SimulationSignaller nullSignaller(nullptr, nullptr, nullptr, false, false); + + GMX_LOG(mdLog_.info) + .asParagraph() + .appendText( + "Note that it is planned that the command gmx mdrun -rerun will " + "be available in a different form in a future version of GROMACS, " + "e.g. gmx rerun -f."); + + if (ir->efep != FreeEnergyPerturbationType::No + && (mdAtoms_->mdatoms()->nMassPerturbed > 0 || (constr_ && constr_->havePerturbedConstraints()))) + { + gmx_fatal(FARGS, + "Perturbed masses or constraints are not supported by rerun. " + "Either make a .tpr without mass and constraint perturbation, " + "or use GROMACS 2018.4, 2018.5 or later 2018 version."); + } + if (ir->bExpanded) + { + gmx_fatal(FARGS, "Expanded ensemble not supported by rerun."); + } + if (ir->bSimTemp) + { + gmx_fatal(FARGS, "Simulated tempering not supported by rerun."); + } + if (ir->bDoAwh) + { + gmx_fatal(FARGS, "AWH not supported by rerun."); + } + if (replExParams_.exchangeInterval > 0) + { + gmx_fatal(FARGS, "Replica exchange not supported by rerun."); + } + if (opt2bSet("-ei", nFile_, fnm_) || observablesHistory_->edsamHistory != nullptr) + { + gmx_fatal(FARGS, "Essential dynamics not supported by rerun."); + } + if (ir->bIMD) + { + gmx_fatal(FARGS, "Interactive MD not supported by rerun."); + } + if (isMultiSim(ms_)) + { + gmx_fatal(FARGS, "Multiple simulations not supported by rerun."); + } + if (std::any_of(ir->opts.annealing, + ir->opts.annealing + ir->opts.ngtc, + [](SimulatedAnnealing i) { return i != SimulatedAnnealing::No; })) + { + gmx_fatal(FARGS, "Simulated annealing not supported by rerun."); + } + + /* Rerun can't work if an output file name is the same as the input file name. + * If this is the case, the user will get an error telling them what the issue is. + */ + if (std::strcmp(opt2fn("-rerun", nFile_, fnm_), opt2fn("-o", nFile_, fnm_)) == 0 + || std::strcmp(opt2fn("-rerun", nFile_, fnm_), opt2fn("-x", nFile_, fnm_)) == 0) + { + gmx_fatal(FARGS, + "When using mdrun -rerun, the name of the input trajectory file " + "%s cannot be identical to the name of an output file (whether " + "given explicitly with -o or -x, or by default)", + opt2fn("-rerun", nFile_, fnm_)); + } + + /* Settings for rerun */ + { + // TODO: Avoid changing inputrec (#3854) + auto* nonConstInputrec = const_cast(inputRec_); + nonConstInputrec->nstlist = 1; + nonConstInputrec->nstcalcenergy = 1; + nonConstInputrec->nstxout_compressed = 0; + } + int nstglobalcomm = 1; + const bool bNS = true; + + ObservablesReducer observablesReducer = observablesReducerBuilder_->build(); + + const SimulationGroups* groups = &topGlobal_.groups; + if (ir->eI == IntegrationAlgorithm::Mimic) + { + auto* nonConstGlobalTopology = const_cast(&topGlobal_); + nonConstGlobalTopology->intermolecularExclusionGroup = genQmmmIndices(topGlobal_); + } + int* fep_state = isMainRank ? &stateGlobal_->fep_state : nullptr; + gmx::ArrayRef lambda = isMainRank ? stateGlobal_->lambda : gmx::ArrayRef{}; + initialize_lambdas( + fpLog_, ir->efep, ir->bSimTemp, *ir->fepvals, ir->simtempvals->temperatures, ekind_, isMainRank, fep_state, lambda); + const bool simulationsShareState = false; + gmx_mdoutf* outf = init_mdoutf(fpLog_, + nFile_, + fnm_, + mdrunOptions_, + cr_, + outputProvider_, + mdModulesNotifiers_, + ir, + topGlobal_, + oenv_, + wallCycleCounters_, + StartingBehavior::NewSimulation, + simulationsShareState, + ms_); + gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), + topGlobal_, + *ir, + pullWork_, + mdoutf_get_fp_dhdl(outf), + true, + StartingBehavior::NewSimulation, + simulationsShareState, + mdModulesNotifiers_); + + gstat = global_stat_init(ir); + + /* Check for polarizable models and flexible constraints */ + shellfc = init_shell_flexcon(fpLog_, + topGlobal_, + constr_ ? constr_->numFlexibleConstraints() : 0, + ir->nstcalcenergy, + haveDDAtomOrdering(*cr_), + runScheduleWork_->simulationWork); + + if (haveDDAtomOrdering(*cr_)) + { + // Local state only becomes valid now. + dd_init_local_state(*cr_->dd, stateGlobal_, state_); + + /* Distribute the charge groups over the nodes from the main node */ + dd_partition_system(fpLog_, + mdLog_, + ir->init_step, + cr_->dd, + TRUE, + stateGlobal_, + topGlobal_, + *ir, + mdModulesNotifiers_, + imdSession_, + pullWork_, + state_, + &f, + mdAtoms_, + top_, + fr_, + virtualSites_, + constr_, + nrnb_, + nullptr, + FALSE); + } + else + { + /* Copy the pointer to the global state */ + state_ = stateGlobal_; + + mdAlgorithmsSetupAtomData( + cr_->dd, *ir, topGlobal_, top_, fr_, &f, mdAtoms_, constr_, virtualSites_, shellfc); + } + + auto* mdatoms = mdAtoms_->mdatoms(); + fr_->longRangeNonbondeds->updateAfterPartition(*mdatoms); + + // NOTE: The global state is no longer used at this point. + // But state_global is still used as temporary storage space for writing + // the global state to file and potentially for replica exchange. + // (Global topology should persist.) + + update_mdatoms(mdatoms, state_->lambda[FreeEnergyPerturbationCouplingType::Mass]); + + if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl != 0) + { + doFreeEnergyPerturbation = true; + } + + int64_t step = ir->init_step; + int64_t step_rel = 0; + + { + int cglo_flags = CGLO_GSTAT; + bool bSumEkinhOld = false; + t_vcm* vcm = nullptr; + compute_globals(gstat, + cr_->commMyGroup, + ir, + fr_, + ekind_, + makeConstArrayRef(state_->x), + makeConstArrayRef(state_->v), + state_->box, + mdatoms, + nrnb_, + vcm, + nullptr, + enerd_, + force_vir, + shake_vir, + total_vir, + pres, + &nullSignaller, + state_->box, + &bSumEkinhOld, + cglo_flags, + step, + &observablesReducer); + // Clean up after pre-step use of compute_globals() + observablesReducer.markAsReadyToReduce(); + } + + if (isMainRank) + { + fprintf(stderr, + "starting md rerun '%s', reading coordinates from" + " input trajectory '%s'\n\n", + *(topGlobal_.name), + opt2fn("-rerun", nFile_, fnm_)); + if (mdrunOptions_.verbose) + { + fprintf(stderr, + "Calculated time to finish depends on nsteps from " + "run input file,\nwhich may not correspond to the time " + "needed to process input trajectory.\n\n"); + } + fprintf(fpLog_, "\n"); + } + + walltime_accounting_start_time(wallTimeAccounting_); + wallcycle_start(wallCycleCounters_, WallCycleCounter::Run); + print_start(fpLog_, cr_, wallTimeAccounting_, "mdrun"); + + /*********************************************************** + * + * Loop over MD steps + * + ************************************************************/ + + if (constr_) + { + GMX_LOG(mdLog_.info) + .asParagraph() + .appendText("Simulations has constraints. Rerun does not recalculate constraints."); + } + + rerun_fr.natoms = 0; + if (isMainRank) + { + isLastStep = !read_first_frame(oenv_, &status, opt2fn("-rerun", nFile_, fnm_), &rerun_fr, TRX_NEED_X); + if (rerun_fr.natoms != topGlobal_.natoms) + { + gmx_fatal(FARGS, + "Number of atoms in trajectory (%d) does not match the " + "run input file (%d)\n", + rerun_fr.natoms, + topGlobal_.natoms); + } + + if (ir->pbcType != PbcType::No) + { + if (!rerun_fr.bBox) + { + gmx_fatal(FARGS, + "Rerun trajectory frame step %" PRId64 + " time %f " + "does not contain a box, while pbc is used", + rerun_fr.step, + rerun_fr.time); + } + if (max_cutoff2(ir->pbcType, rerun_fr.box) < gmx::square(fr_->rlist)) + { + gmx_fatal(FARGS, + "Rerun trajectory frame step %" PRId64 + " time %f " + "has too small box dimensions", + rerun_fr.step, + rerun_fr.time); + } + } + } + + GMX_LOG(mdLog_.info) + .asParagraph() + .appendText( + "Rerun does not report kinetic energy, total energy, temperature, virial and " + "pressure."); + + if (cr_->commMyGroup.isParallel()) + { + rerun_parallel_comm(cr_->commMyGroup, &rerun_fr, &isLastStep); + } + + if (ir->pbcType != PbcType::No) + { + /* Set the shift vectors. + * Necessary here when have a static box different from the tpr box. + */ + calc_shifts(rerun_fr.box, fr_->shift_vec); + } + + auto stopHandler = stopHandlerBuilder_->getStopHandlerMD( + compat::not_null(&signals[eglsSTOPCOND]), + false, + isMainRank, + 1, // rerun constructs the pairlist for each frame + mdrunOptions_.reproducible, + nstglobalcomm, + mdrunOptions_.maximumHoursToRun, + fpLog_, + step, + wallTimeAccounting_); + + // we don't do counter resetting in rerun - finish will always be valid + walltime_accounting_set_valid_finish(wallTimeAccounting_); + + const DDBalanceRegionHandler ddBalanceRegionHandler(cr_->dd); + + /* and stop now if we should */ + isLastStep = (isLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps)); + while (!isLastStep) + { + wallcycle_start(wallCycleCounters_, WallCycleCounter::Step); + + if (rerun_fr.bStep) + { + step = rerun_fr.step; + step_rel = step - ir->init_step; + } + if (rerun_fr.bTime) + { + t = rerun_fr.time; + } + else + { + t = step; + } + + if (ir->efep != FreeEnergyPerturbationType::No && isMainRank) + { + if (rerun_fr.bLambda) + { + ir->fepvals->init_lambda_without_states = rerun_fr.lambda; + } + else + { + if (rerun_fr.bFepState) + { + state_->fep_state = rerun_fr.fep_state; + } + } + + stateGlobal_->lambda = currentLambdas(step, *(ir->fepvals), state_->fep_state); + } + + if (isMainRank) + { + const bool constructVsites = + ((virtualSites_ != nullptr) && mdrunOptions_.rerunConstructVsites); + if (constructVsites && haveDDAtomOrdering(*cr_)) + { + gmx_fatal(FARGS, + "Vsite recalculation with -rerun is not implemented with domain " + "decomposition, " + "use a single rank"); + } + prepareRerunState(rerun_fr, stateGlobal_, constructVsites, virtualSites_); + } + + isLastStep = isLastStep || stopHandler->stoppingAfterCurrentStep(step); + + if (haveDDAtomOrdering(*cr_)) + { + /* Repartition the domain decomposition */ + const bool bMainState = true; + dd_partition_system(fpLog_, + mdLog_, + step, + cr_->dd, + bMainState, + stateGlobal_, + topGlobal_, + *ir, + mdModulesNotifiers_, + imdSession_, + pullWork_, + state_, + &f, + mdAtoms_, + top_, + fr_, + virtualSites_, + constr_, + nrnb_, + wallCycleCounters_, + mdrunOptions_.verbose); + } + + if (isMainRank) + { + EnergyOutput::printHeader(fpLog_, step, t); /* can we improve the information printed here? */ + } + + if (ir->efep != FreeEnergyPerturbationType::No) + { + update_mdatoms(mdatoms, state_->lambda[FreeEnergyPerturbationCouplingType::Mass]); + } + + fr_->longRangeNonbondeds->updateAfterPartition(*mdatoms); + + force_flags = (GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES + | GMX_FORCE_VIRIAL + // TODO: Get rid of this once #2649 and #3400 are solved + | GMX_FORCE_ENERGY | (doFreeEnergyPerturbation ? GMX_FORCE_DHDL : 0)); + + const int shellfcFlags = force_flags | (mdrunOptions_.verbose ? GMX_FORCE_ENERGY : 0); + const int legacyForceFlags = ((shellfc) ? shellfcFlags : force_flags) | GMX_FORCE_NS; + + gmx_edsam* const ed = nullptr; + + if (bNS) + { + if (fr_->listedForcesGpu) + { + fr_->listedForcesGpu->updateHaveInteractions(top_->idef); + } + runScheduleWork_->domainWork = setupDomainLifetimeWorkload( + *ir, *fr_, pullWork_, ed, *mdatoms, runScheduleWork_->simulationWork); + } + + + runScheduleWork_->stepWork = setupStepWorkload(legacyForceFlags, + ir->mtsLevels, + step, + runScheduleWork_->domainWork, + runScheduleWork_->simulationWork); + + if (shellfc) + { + /* Now is the time to relax the shells */ + relax_shell_flexcon(fpLog_, + cr_, + mdrunOptions_.verbose, + enforcedRotation_, + step, + ir, + mdModulesNotifiers_, + imdSession_, + pullWork_, + bNS, + top_, + constr_, + enerd_, + state_->numAtoms(), + state_->x.arrayRefWithPadding(), + state_->v.arrayRefWithPadding(), + state_->box, + state_->lambda, + &state_->hist, + &f.view(), + force_vir, + *mdatoms, + fr_->longRangeNonbondeds.get(), + nrnb_, + wallCycleCounters_, + shellfc, + fr_, + *runScheduleWork_, + t, + mu_tot, + virtualSites_, + ddBalanceRegionHandler); + } + else + { + /* The coordinates (x) are shifted (to get whole molecules) + * in do_force. + * This is parallellized as well, and does communication too. + * Check comments in sim_util.c + */ + Awh* awh = nullptr; + + try + { + do_force(fpLog_, + cr_, + *ir, + mdModulesNotifiers_, + awh, + enforcedRotation_, + imdSession_, + pullWork_, + step, + nrnb_, + wallCycleCounters_, + top_, + state_->box, + state_->x.arrayRefWithPadding(), + {}, + &state_->hist, + &f.view(), + force_vir, + mdatoms, + enerd_, + state_->lambda, + fr_, + *runScheduleWork_, + virtualSites_, + mu_tot, + t, + ed, + fr_->longRangeNonbondeds.get(), + ddBalanceRegionHandler); + } + catch (const gmx::InternalError&) + { + GMX_LOG(mdLog_.warning) + .asParagraph() + .appendText( + "Continuing with next frame after catching invalid force in " + "previous frame"); + }; + } + + /* Now we have the energies and forces corresponding to the + * coordinates at time t. + */ + { + const bool isCheckpointingStep = false; + const bool doRerun = true; + do_md_trajectory_writing(fpLog_, + cr_, + nFile_, + fnm_, + step, + step_rel, + t, + ir, + state_, + stateGlobal_, + observablesHistory_, + topGlobal_, + fr_, + outf, + energyOutput, + ekind_, + f.view().force(), + isCheckpointingStep, + doRerun, + isLastStep, + mdrunOptions_.writeConfout, + EkindataState::NotUsed); + } + + stopHandler->setSignal(); + + { + const bool doInterSimSignal = false; + const bool doIntraSimSignal = true; + bool bSumEkinhOld = false; + t_vcm* vcm = nullptr; + SimulationSignaller signaller(&signals, cr_, ms_, doInterSimSignal, doIntraSimSignal); + + int cglo_flags = CGLO_GSTAT | CGLO_ENERGY; + compute_globals(gstat, + cr_->commMyGroup, + ir, + fr_, + ekind_, + makeConstArrayRef(state_->x), + makeConstArrayRef(state_->v), + state_->box, + mdatoms, + nrnb_, + vcm, + wallCycleCounters_, + enerd_, + force_vir, + shake_vir, + total_vir, + pres, + &signaller, + state_->box, + &bSumEkinhOld, + cglo_flags, + step, + &observablesReducer); + // Clean up after pre-step use of compute_globals() + observablesReducer.markAsReadyToReduce(); + } + + /* Note: this is OK, but there are some numerical precision issues with using the convergence of + the virial that should probably be addressed eventually. state->veta has better properties, + but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could + generate the new shake_vir, but test the veta value for convergence. This will take some thought. */ + + /* Output stuff */ + if (isMainRank) + { + const bool bCalcEnerStep = true; + energyOutput.addDataAtEnergyStep(doFreeEnergyPerturbation, + bCalcEnerStep, + t, + mdatoms->tmass, + enerd_, + ir->fepvals.get(), + state_->box, + PTCouplingArrays({ state_->boxv, + state_->nosehoover_xi, + state_->nosehoover_vxi, + state_->nhpres_xi, + state_->nhpres_vxi }), + state_->fep_state, + total_vir, + pres, + ekind_, + mu_tot, + constr_); + + const bool do_ene = true; + const bool do_log = true; + Awh* awh = nullptr; + const bool do_dr = ir->nstdisreout != 0; + const bool do_or = ir->nstorireout != 0; + + EnergyOutput::printAnnealingTemperatures(do_log ? fpLog_ : nullptr, *groups, ir->opts, *ekind_); + energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), + do_ene, + do_dr, + do_or, + do_log ? fpLog_ : nullptr, + step, + t, + fr_->fcdata.get(), + awh); + + if (ir->bPull) + { + pull_print_output(pullWork_, step, t); + } + + if (do_per_step(step, ir->nstlog)) + { + if (std::fflush(fpLog_) != 0) + { + gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?"); + } + } + } + + /* Print the remaining wall clock time for the run */ + if (isMainSimMainRank(ms_, isMainRank) && (mdrunOptions_.verbose || gmx_got_usr_signal())) + { + if (shellfc) + { + fprintf(stderr, "\n"); + } + print_time(stderr, wallTimeAccounting_, step, ir, cr_->commMySim); + } + + /* Ion/water position swapping. + * Not done in last step since trajectory writing happens before this call + * in the MD loop and exchanges would be lost anyway. */ + if ((ir->eSwapCoords != SwapType::No) && (step > 0) && !isLastStep + && do_per_step(step, ir->swap->nstswap)) + { + const bool doRerun = true; + do_swapcoords(cr_->commMyGroup, + step, + t, + ir, + swap_, + wallCycleCounters_, + gmx::arrayRefFromArray(reinterpret_cast(rerun_fr.x), rerun_fr.natoms), + rerun_fr.box, + isMainRank && mdrunOptions_.verbose, + doRerun); + } + + if (isMainRank) + { + /* read next frame from input trajectory */ + isLastStep = !read_next_frame(oenv_, status, &rerun_fr); + } + + if (cr_->commMyGroup.isParallel()) + { + rerun_parallel_comm(cr_->commMyGroup, &rerun_fr, &isLastStep); + } + + cycles = wallcycle_stop(wallCycleCounters_, WallCycleCounter::Step); + if (haveDDAtomOrdering(*cr_) && wallCycleCounters_) + { + dd_cycles_add(cr_->dd, cycles, ddCyclStep); + } + + if (!rerun_fr.bStep) + { + /* increase the MD step number */ + step++; + step_rel++; + } + observablesReducer.markAsReadyToReduce(); + } + /* End of main MD loop */ + + /* Closing TNG files can include compressing data. Therefore it is good to do that + * before stopping the time measurements. */ + mdoutf_tng_close(outf); + + /* Stop measuring walltime */ + walltime_accounting_end_time(wallTimeAccounting_); + + if (isMainRank) + { + close_trx(status); + } + + if (!thisRankHasPmeDuty(cr_->dd)) + { + /* Tell the PME only node to finish */ + gmx_pme_send_finish(cr_->dd); + } + + done_mdoutf(outf); + + done_shellfc(fpLog_, shellfc, step_rel); + + walltime_accounting_set_nsteps_done(wallTimeAccounting_, step_rel); +} From 74b08d949869e7cc9143fbbbe3e0446c39778dad Mon Sep 17 00:00:00 2001 From: Michele Invernizzi Date: Wed, 14 Oct 2026 18:04:35 +0200 Subject: [PATCH 09/10] GREX: complete the exchange of positions when saving them GREX::savePositions() shares the positions of all atoms (shareAll), but then waits only for the inputs that are active on this step (waitData). When no action needs atoms on an exchange step, the domain decomposition started its exchange but never completed it. With asynchronous sharing (the default on fewer than 10 ranks), the sends were never received: the next exchange waited for them forever with large messages, and with small ones the following receives could get the positions of an earlier step. PLUMED 2.9 always waited. Found by a GROMACS regression test: -hrex with two TIP4P water replicas on 2 ranks each hung at the step after the first exchange attempt with the PLUMED 2.10 kernel, and ran with the 2.9 one. Co-Authored-By: Claude Opus 5.5 --- CHANGES/v2.10.md | 3 +++ src/core/GREX.cpp | 10 ++++++++++ 2 files changed, 13 insertions(+) diff --git a/CHANGES/v2.10.md b/CHANGES/v2.10.md index 04d0a51440..0565c516eb 100644 --- a/CHANGES/v2.10.md +++ b/CHANGES/v2.10.md @@ -116,6 +116,9 @@ Before switching to version 2.10, users are invited to carefully read the follow - 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 diff --git a/src/core/GREX.cpp b/src/core/GREX.cpp index 007cdc6c7b..1fe322dfd2 100644 --- a/src/core/GREX.cpp +++ b/src/core/GREX.cpp @@ -21,6 +21,8 @@ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ #include "GREX.h" #include "PlumedMain.h" +#include "ActionSet.h" +#include "DomainDecomposition.h" #include "tools/Tools.h" #include "tools/Communicator.h" #include @@ -188,6 +190,14 @@ void GREX::savePositions() { plumedMain.resetActive(true); plumedMain.shareAll(); plumedMain.waitData(); + // shareAll() also starts the exchange of the domain decomposition when no action needs atoms on this + // step, but waitData() only completes the exchanges of the active inputs. Complete it here, or the + // next exchange waits for messages that were never received (with asynchronous sharing). + for(const auto & dd : plumedMain.getActionSet().select()) { + if(!dd->isActive() && dd->hasBeenSet()) { + dd->wait(); + } + } std::ostringstream o; plumedMain.writeBinary(o); buffer=o.str(); From 608e30de631d0d8372ba919206ff1209f69d12b0 Mon Sep 17 00:00:00 2001 From: Michele Invernizzi Date: Thu, 15 Oct 2026 09:22:18 +0200 Subject: [PATCH 10/10] Docs: what the GROMACS interfaces support A table in the list of MD codes compares the PLUMED module of GROMACS 2025 and later with the patches for GROMACS 2025/2026 and for the older versions (ENERGY, multiple walkers, checkpoints, stop requests, shells, replica exchange, lambda dynamics), as tested by the regression tests of the PLUMED testcenter. The descriptions of the 2025/2026 patches also list stop requests and shells. Co-Authored-By: Claude Opus 5.5 --- patches/gromacs-2025.0.config | 5 ++++- patches/gromacs-2026.0.config | 7 +++++-- user-doc/Introduction.md | 18 +++++++++++++++++- 3 files changed, 26 insertions(+), 4 deletions(-) diff --git a/patches/gromacs-2025.0.config b/patches/gromacs-2025.0.config index d1115d0b4f..b5df204a4e 100644 --- a/patches/gromacs-2025.0.config +++ b/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. diff --git a/patches/gromacs-2026.0.config b/patches/gromacs-2026.0.config index 56c3cd7147..1a0baff93e 100644 --- a/patches/gromacs-2026.0.config +++ b/patches/gromacs-2026.0.config @@ -9,7 +9,7 @@ function plumed_patch_info(){ cat << EOF A basic PLUMED interface is already present in GROMACS -This patch previes extra features to be implemented in the official PLUMED integration in the next GROMACS version, including: +This patch provides extra features to be implemented in the official PLUMED integration in the next GROMACS version, including: - multiple walkers (GREX) - OpenMP thread count passed to PLUMED @@ -17,10 +17,13 @@ This patch previes extra features to be implemented in the official PLUMED integ 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. diff --git a/user-doc/Introduction.md b/user-doc/Introduction.md index c82c4f7ec4..3edb74f2e2 100644 --- a/user-doc/Introduction.md +++ b/user-doc/Introduction.md @@ -46,7 +46,7 @@ As far as we know, the following MD codes can be used with PLUMED out of the box - [Metalwalls](https://gitlab.com/ampere2/metalwalls) - [ASE](https://wiki.fysik.dtu.dk/ase/) - [GPUMD](https://gpumd.org/) -- [GROMACS](https://www.gromacs.org/) as of version 2025, with limited support (no replica exchange, no \ref ENERGY collective variable, no lambda dynamics). +- [GROMACS](https://www.gromacs.org/) as of version 2025, with limited support (see the table below). Please refer to the documentation of the MD code to know how to use it with the latest PLUMED release. If you maintain another MD code that is PLUMED-ready let us know and we will add it to this list. @@ -55,6 +55,22 @@ Additionally, we provide patching procedures for the following codes: @CODESL@ +The patches for GROMACS 2025 and 2026 add features to the PLUMED interface that GROMACS includes since version 2025, +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 | +| PLUMED stopping the simulation (e.g. \ref COMMITTOR) | no | yes | yes | +| polarizable models with shells | no | yes | no | +| replica exchange with a PLUMED bias (`-replex`, `-hrex`) | no | no | yes | +| lambda dynamics (expanded ensemble) | no | no | yes | + +The interface with GROMACS is tested on the [PLUMED testcenter](https://plumed-testcenter.github.io/browse.html). + Alternatively, one can use PLUMED as a \ref tools for post processing the results from molecular dynamics or enhanced sampling calculations. Notice that PLUMED can be used as an analysis tool