diff --git a/CHANGES/v2.10.md b/CHANGES/v2.10.md index 78d6a4d29f..0565c516eb 100644 --- a/CHANGES/v2.10.md +++ b/CHANGES/v2.10.md @@ -113,3 +113,26 @@ Before switching to version 2.10, users are invited to carefully read the follow - Fixed the GROMACS 2025 patch so that `-multidir` simulations without `-plumed` do not have to run in lockstep, which made them deadlock when they had different numbers of steps. As in the GROMACS 2026 patch, the multi-simulation record is now used only when PLUMED is active. + - Fixed the virial with a bias on ENERGY (e.g. OPES multithermal or the well-tempered ensemble) and + domain decomposition: only the first rank rescaled its share of the virial of the MD code, so the + pressure was wrong, by orders of magnitude in some cases. As in PLUMED 2.9, it is now rescaled on all ranks. + - Replica exchange with more than one rank per replica (e.g. GROMACS -replex or -hrex) hung after an exchange + attempt on which no action needed the atoms: the exchange of positions started for the exchange was never + completed. With small systems, the following exchanges could instead receive the positions of an earlier step. + - GROMACS 2025 and 2026 patches: PLUMED can now stop the simulation, e.g. with COMMITTOR. + The request reaches all ranks and the simulations that share state, and mdrun stops at the next + step at which GROMACS communicates signals (every nstglobalcomm steps, and at least every 200 steps + between simulations that share state; with the modular simulator, at the next neighbor-search step + after that), writing a checkpoint there. + - GROMACS 2025 and 2026 patches: when the PLUMED input needs the potential energy, the forces + that PLUMED applied on virtual sites were lost. They are now spread on the constructing atoms, + on the steps when PLUMED uses virtual sites. + - GROMACS 2025 and 2026 patches: with multiple time stepping, the potential energy passed to PLUMED + on the fast steps lacked the slow terms (e.g. the PME mesh, by 650-1250 kJ/mol in a water test). + mdrun now stops with an error if PLUMED needs the energy on such a step. + - GROMACS 2025 and 2026 patches: errors raised on purpose by the PLUMED module, e.g. for replica + exchange with a bias, are no longer reported as internal errors. + - GROMACS 2025 and 2026 patches: a request of PLUMED to stop (e.g. \ref COMMITTOR) was ignored in reruns. + - GROMACS 2025 and 2026 patches: with shell (polarizable) particles, PLUMED also ran on each shell relaxation + step, as step 1: e.g. \ref PRINT printed extra lines. It now only calculates its forces again on those steps, + without updating. 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-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 { 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 34b3ba7c22..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 @@ -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: @@ -147,6 +151,11 @@ try : plumed_(std::make_unique()),replex_(options.replex_) plumed_->cmd("setTimestep", &options.simulationTimeStep_); plumed_->cmd("init", nullptr); } +catch (const GromacsException&) +{ + // Keep the type and message of the errors raised on purpose, e.g. NotImplementedError + throw; +} catch (const std::exception& ex) { GMX_THROW(InternalError( @@ -161,6 +170,11 @@ try checkpointRequested_ = 1; } } +catch (const GromacsException&) +{ + // Keep the type and message of the errors raised on purpose, e.g. NotImplementedError + throw; +} catch (const std::exception& ex) { GMX_THROW(InternalError( @@ -174,7 +188,12 @@ try // setup: these instructions in the original patch are BEFORE do_force() // now this is called within do_force(), but this does not impact the results const t_commrec* cr = &(forceProviderInput.cr_); - long int lstep = forceProviderInput.step_; + long int lstep = prepareStepIsCalled_ ? currentStep_ : forceProviderInput.step_; + // With shell particles, relax_shell_flexcon() calls do_force() again on the same step (passing + // step 1): PLUMED then calculates its forces again, without updating (depositing, printing...). + // Only in MD, where prepareStep() gives the step: e.g. reruns can have frames with the same step + repeatedCall_ = prepareStepIsCalled_ && (calculatedStep_ == lstep); + calculatedStep_ = lstep; plumed_->cmd("setStepLong", &lstep); if (haveDDAtomOrdering(*cr)) { @@ -201,7 +220,10 @@ try } preparedStep_.reset(); - plumedWantsToStop_ = 0; + if (!repeatedCall_) + { + plumedWantsToStop_ = 0; + } plumed_->cmd("setStopFlag", &plumedWantsToStop_); if (!prepareStepIsCalled_) @@ -229,19 +251,42 @@ try clear_mat(plumedVir_); plumed_->cmd("setVirial", &plumedVir_[0][0]); - plumed_->cmd("performCalc", nullptr); + plumed_->cmd(repeatedCall_ ? "performCalcNoUpdate" : "performCalc", nullptr); checkReplicaExchangeBias(); msmul(plumedVir_, 0.5, plumedVir_); forceProviderOutput->forceWithVirial_.addVirialContribution(plumedVir_); } +catch (const GromacsException&) +{ + // Keep the type and message of the errors raised on purpose, e.g. NotImplementedError + throw; +} catch (const std::exception& ex) { GMX_THROW(InternalError( std::string("An error occurred while PLUMED was calculating the forces\n:") + ex.what())); } +bool PlumedForceProvider::usesVirtualSites() +{ + if (isVirtualSite_.empty()) + { + return false; + } + /* PLUMED applies forces only on the atoms it uses. Their list is the same on all ranks, + so that all take the same decision. It is empty when PLUMED only biases the energy. */ + int numAtoms = 0; + plumed_->cmd("createFullList", &numAtoms); + const int* atoms = nullptr; + plumed_->cmd("getFullList", &atoms); + const bool found = + std::any_of(atoms, atoms + numAtoms, [this](int atom) { return isVirtualSite_[atom]; }); + plumed_->cmd("clearFullList", nullptr); + return found; +} + bool PlumedForceProvider::prepareStep(int64_t step, bool isCheckpointingStep) try { @@ -252,6 +297,7 @@ try long int lstep = step; plumed_->cmd("setStepLong", &lstep); prepareStepIsCalled_ = true; + currentStep_ = step; if (plumedAPIversion_ > 3) { // Set on every step: PLUMED resets it only on the steps when some action is active @@ -264,12 +310,22 @@ try preparedStep_ = step; return plumedNeedsEnergy_ != 0; } +catch (const GromacsException&) +{ + // Keep the type and message of the errors raised on purpose, e.g. NotImplementedError + throw; +} catch (const std::exception& ex) { GMX_THROW(InternalError( std::string("An error occurred while PLUMED was preparing the step\n:") + ex.what())); } +bool PlumedForceProvider::requestsStop() const +{ + return plumedWantsToStop_ != 0; +} + void PlumedForceProvider::checkReplicaExchangeBias() { if (replex_) @@ -284,7 +340,7 @@ void PlumedForceProvider::checkReplicaExchangeBias() } } -void PlumedForceProvider::applyAfterPotentialEnergy(bool energyWasComputed, +bool PlumedForceProvider::applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef force, tensor virial) @@ -292,7 +348,7 @@ try { if (!plumedNeedsEnergy_) { - return; + return false; } if (!energyWasComputed || potentialEnergy == nullptr) { @@ -312,11 +368,17 @@ try // Match the legacy GROMACS patch: pass 2*virial into PLUMED, then replace. msmul(virial, 2.0, plumedVir_); plumed_->cmd("setVirial", &plumedVir_[0][0]); - plumed_->cmd("performCalc", nullptr); + plumed_->cmd(repeatedCall_ ? "performCalcNoUpdate" : "performCalc", nullptr); checkReplicaExchangeBias(); msmul(plumedVir_, 0.5, virial); + return usesVirtualSites(); +} +catch (const GromacsException&) +{ + // Keep the type and message of the errors raised on purpose, e.g. NotImplementedError + throw; } catch (const std::exception& ex) { 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..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 @@ -43,6 +43,7 @@ #include #include +#include #include "gromacs/mdtypes/iforceprovider.h" @@ -81,17 +82,22 @@ class PlumedForceProvider final : public IForceProvider void calculateForces(const ForceProviderInput& forceProviderInput, ForceProviderOutput* forceProviderOutput) override; bool prepareStep(int64_t step, bool isCheckpointingStep) override; - void applyAfterPotentialEnergy(bool energyWasComputed, + bool requestsStop() const override; + bool applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef 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 } }; @@ -101,6 +107,12 @@ class PlumedForceProvider final : public IForceProvider bool prepareStepIsCalled_ = false; //! Set on the main rank by writeCheckpointData(), passed to PLUMED by calculateForces(). int checkpointRequested_ = 0; + //! Step for which prepareStep() was last called + int64_t currentStep_ = 0; + //! Last step that PLUMED calculated + std::optional 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-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/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-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-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..c53d8c4879 100644 --- a/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.cpp +++ b/patches/gromacs-2025.0.diff/src/gromacs/mdtypes/iforceprovider.cpp @@ -109,13 +109,29 @@ bool ForceProviders::prepareStep(int64_t step, bool isCheckpointingStep) const return needed; } -void ForceProviders::applyAfterPotentialEnergy(bool energyWasComputed, +bool ForceProviders::requestsStop() const +{ + for (const auto& provider : impl_->providers_) + { + if (provider.first->requestsStop()) + { + return true; + } + } + return false; +} + +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 56916fd4f3..e8adf8c1cb 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. @@ -204,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: @@ -252,8 +263,12 @@ class ForceProviders //! Prepares all providers for \p step, and returns whether any needs the potential energy. bool prepareStep(int64_t step, bool isCheckpointingStep) const; - //! Forwards applyAfterPotentialEnergy() to all providers. - void applyAfterPotentialEnergy(bool energyWasComputed, + //! Whether any provider asks to stop the simulation. + bool requestsStop() const; + + //! 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.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/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 8251a0394d..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 @@ -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: @@ -162,6 +165,11 @@ try : plumed_(std::make_unique()), replex_(options.replex_) plumed_->cmd("setTimestep", &options.simulationTimeStep_); plumed_->cmd("init", nullptr); } +catch (const GromacsException&) +{ + // Keep the type and message of the errors raised on purpose, e.g. NotImplementedError + throw; +} catch (const std::exception& ex) { GMX_THROW(InternalError( @@ -181,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( @@ -193,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()) { @@ -220,7 +238,10 @@ try } preparedStep_.reset(); - plumedWantsToStop_ = 0; + if (!repeatedCall_) + { + plumedWantsToStop_ = 0; + } plumed_->cmd("setStopFlag", &plumedWantsToStop_); if (!prepareStepIsCalled_) @@ -255,19 +276,29 @@ try // in the original patch do_force() was called HERE // Do the work - plumed_->cmd("performCalc", nullptr); + plumed_->cmd(repeatedCall_ ? "performCalcNoUpdate" : "performCalc", nullptr); checkReplicaExchangeBias(); msmul(plumedVir_, 0.5, plumedVir_); forceProviderOutput->forceWithVirial_.addVirialContribution(plumedVir_); } +catch (const GromacsException&) +{ + // Keep the type and message of the errors raised on purpose, e.g. NotImplementedError + throw; +} catch (const std::exception& ex) { GMX_THROW(InternalError( std::string("An error occurred while PLUMED was calculating the forces\n:") + ex.what())); } +bool PlumedForceProvider::requestsStop() const +{ + return plumedWantsToStop_ != 0; +} + void PlumedForceProvider::checkReplicaExchangeBias() { if (replex_) @@ -282,6 +313,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 { @@ -292,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 @@ -304,13 +354,18 @@ try preparedStep_ = step; return plumedNeedsEnergy_ != 0; } +catch (const GromacsException&) +{ + // Keep the type and message of the errors raised on purpose, e.g. NotImplementedError + throw; +} catch (const std::exception& ex) { GMX_THROW(InternalError( std::string("An error occurred while PLUMED was preparing the step\n:") + ex.what())); } -void PlumedForceProvider::applyAfterPotentialEnergy(bool energyWasComputed, +bool PlumedForceProvider::applyAfterPotentialEnergy(bool energyWasComputed, real* potentialEnergy, ArrayRef force, tensor virial) @@ -318,7 +373,7 @@ try { if (!plumedNeedsEnergy_) { - return; + return false; } if (!energyWasComputed || potentialEnergy == nullptr) { @@ -338,11 +393,17 @@ 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(); 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) { 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..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 @@ -43,6 +43,7 @@ #include #include +#include #include "gromacs/mdtypes/iforceprovider.h" @@ -91,13 +92,16 @@ 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 * \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; @@ -105,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; @@ -120,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 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/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 */ 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); +} 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..c53d8c4879 100644 --- a/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.cpp +++ b/patches/gromacs-2026.0.diff/src/gromacs/mdtypes/iforceprovider.cpp @@ -109,13 +109,29 @@ bool ForceProviders::prepareStep(int64_t step, bool isCheckpointingStep) const return needed; } -void ForceProviders::applyAfterPotentialEnergy(bool energyWasComputed, +bool ForceProviders::requestsStop() const +{ + for (const auto& provider : impl_->providers_) + { + if (provider.first->requestsStop()) + { + return true; + } + } + return false; +} + +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 3084fb526b..64a7fdb20f 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. @@ -208,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: @@ -256,8 +267,12 @@ class ForceProviders //! Prepares all providers for \p step, and returns whether any needs the potential energy. bool prepareStep(int64_t step, bool isCheckpointingStep) const; - //! Forwards applyAfterPotentialEnergy() to all providers. - void applyAfterPotentialEnergy(bool energyWasComputed, + //! Whether any provider asks to stop the simulation. + bool requestsStop() const; + + //! 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/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; } } 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(); 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