From 215dd74563c634ce41e32ef2007caf840f90bfc7 Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Mon, 2 Dec 2024 14:45:23 +0100 Subject: [PATCH 1/7] setting up the patch for gmx 2025 --- patches/gromacs-2025.0.config | 34 +++ .../applied_forces/plumed/plumedMDModule.cpp | 158 +++++++++++++ .../plumed/plumedMDModule.cpp.preplumed | 154 +++++++++++++ .../applied_forces/plumed/plumedOptions.cpp | 104 +++++++++ .../plumed/plumedOptions.cpp.preplumed | 99 ++++++++ .../applied_forces/plumed/plumedOptions.h | 117 ++++++++++ .../plumed/plumedOptions.h.preplumed | 111 +++++++++ .../plumed/plumedforceprovider.cpp | 215 ++++++++++++++++++ .../plumed/plumedforceprovider.cpp.preplumed | 206 +++++++++++++++++ 9 files changed, 1198 insertions(+) create mode 100644 patches/gromacs-2025.0.config create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedMDModule.cpp create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedMDModule.cpp.preplumed create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.cpp create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.cpp.preplumed create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h.preplumed create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp.preplumed diff --git a/patches/gromacs-2025.0.config b/patches/gromacs-2025.0.config new file mode 100644 index 0000000000..b998dcb06d --- /dev/null +++ b/patches/gromacs-2025.0.config @@ -0,0 +1,34 @@ + + +function plumed_preliminary_test(){ +# check if the README contains the word GROMACS and if gromacs has been already configured + grep -q GROMACS README 1>/dev/null 2>/dev/null +} + +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 + +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 +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. + +To enable PLUMED in a gromacs simulation one should use +mdrun with an extra -plumed flag. The flag can be used to +specify the name of the PLUMED input file, e.g.: + +gmx mdrun -plumed plumed.dat + +For more information on gromacs you should visit http://www.gromacs.org + +EOF +} + +plumed_before_patch(){ + plumed_patch_info +} 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 new file mode 100644 index 0000000000..57ded02a22 --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedMDModule.cpp @@ -0,0 +1,158 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2024- 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 Plumed MDModule class + * + * \author Daniele Rapetti + * \ingroup module_applied_forces + */ +#include "gmxpre.h" + +#include "plumedMDModule.h" + +#include +#include + +#include "gromacs/domdec/localatomsetmanager.h" +#include "gromacs/fileio/checkpoint.h" +#include "gromacs/mdrunutility/mdmodulesnotifiers.h" +#include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/imdmodule.h" +#include "gromacs/utility/keyvaluetreebuilder.h" + +#include "plumedOptions.h" +#include "plumedforceprovider.h" + +namespace gmx +{ + +namespace +{ + +/*! \internal + * \brief Plumed module + * + * Class that implements the plumed MDModule + */ +class PlumedMDModule final : public IMDModule +{ +public: + //! \brief Construct the plumed module. + explicit PlumedMDModule() = default; + // Now callbacks for several kinds of MdModuleNotification are created + // and subscribed, and will be dispatched correctly at run time + // based on the type of the parameter required by the lambda. + + /*! \brief Requests to be notified during pre-processing. + * + * Plumed does not act during the preprocessing phase of a simulation, so the input are ignored + */ + void subscribeToPreProcessingNotifications(MDModulesNotifiers* /*notifier*/) override {} + + /*! \brief Subscribe to MDModules notifications for information needed just before the simulation. + */ + void subscribeToSimulationSetupNotifications(MDModulesNotifiers* notifier) override + { + // TODO: add a check for threadmpi (see #5104, https://gitlab.com/gromacs/gromacs/-/merge_requests/4367#note_2102475958, the manual and the force provider for the details) + + // Access the plumed filename this is used to activate the plumed module + notifier->simulationSetupNotifier_.subscribe( + [this](const PlumedInputFilename& plumedFilename) + { this->options_.setPlumedFile(plumedFilename.plumedFilename_); }); + // Access the temperature if it is constant during the simulation + notifier->simulationSetupNotifier_.subscribe( + [this](const EnsembleTemperature& ensembleT) + { this->options_.setEnsembleTemperature(ensembleT); }); + // Access of the topology + notifier->simulationSetupNotifier_.subscribe([this](const gmx_mtop_t& mtop) + { this->options_.setTopology(mtop); }); + // Retrieve the Communication Record during simulations setup + notifier->simulationSetupNotifier_.subscribe([this](const t_commrec& cr) + { this->options_.setComm(cr); }); + // setting the simulation time step + notifier->simulationSetupNotifier_.subscribe( + [this](const SimulationTimeStep& simulationTimeStep) + { this->options_.setSimulationTimeStep(simulationTimeStep.delta_t); }); + // Retrieve the starting behavior + notifier->simulationSetupNotifier_.subscribe( + [this](const StartingBehavior& startingBehavior) + { this->options_.setStartingBehavior(startingBehavior); }); + // Retrieve the Multisim options + notifier->simulationSetupNotifier_.subscribe( + [this](const gmx_multisim_t* ms) + { this->options_.setMultisim(ms); }); + // writing checkpoint data + notifier->checkpointingNotifier_.subscribe( + [this](MDModulesWriteCheckpointData /*checkpointData*/) + { + if (options_.active()) + { + plumedForceProvider_->writeCheckpointData(); + } + }); + } + + //! From IMDModule + IMdpOptionProvider* mdpOptionProvider() override { return nullptr; } + //! From IMDModule + IMDOutputProvider* outputProvider() override + { // Plumed provide its own output + return nullptr; + } + //! From IMDModule - Adds this module to the force providers if active + void initForceProviders(ForceProviders* forceProviders) override + { + if (options_.active()) + { + plumedForceProvider_ = std::make_unique(options_.options()); + forceProviders->addForceProvider(plumedForceProvider_.get(), "Plumed"); + } + } + + +private: + //! Parameters that become available at simulation setup time. + PlumedOptionProvider options_{}; + //! Object that evaluates the forces + std::unique_ptr plumedForceProvider_{}; +}; + +} // namespace + +std::unique_ptr PlumedModuleInfo::create() +{ + return std::make_unique(); +} +} // namespace gmx diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedMDModule.cpp.preplumed b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedMDModule.cpp.preplumed new file mode 100644 index 0000000000..d5adc14adf --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedMDModule.cpp.preplumed @@ -0,0 +1,154 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2024- 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 Plumed MDModule class + * + * \author Daniele Rapetti + * \ingroup module_applied_forces + */ +#include "gmxpre.h" + +#include "plumedMDModule.h" + +#include +#include + +#include "gromacs/domdec/localatomsetmanager.h" +#include "gromacs/fileio/checkpoint.h" +#include "gromacs/mdrunutility/mdmodulesnotifiers.h" +#include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/imdmodule.h" +#include "gromacs/utility/keyvaluetreebuilder.h" + +#include "plumedOptions.h" +#include "plumedforceprovider.h" + +namespace gmx +{ + +namespace +{ + +/*! \internal + * \brief Plumed module + * + * Class that implements the plumed MDModule + */ +class PlumedMDModule final : public IMDModule +{ +public: + //! \brief Construct the plumed module. + explicit PlumedMDModule() = default; + // Now callbacks for several kinds of MdModuleNotification are created + // and subscribed, and will be dispatched correctly at run time + // based on the type of the parameter required by the lambda. + + /*! \brief Requests to be notified during pre-processing. + * + * Plumed does not act during the preprocessing phase of a simulation, so the input are ignored + */ + void subscribeToPreProcessingNotifications(MDModulesNotifiers* /*notifier*/) override {} + + /*! \brief Subscribe to MDModules notifications for information needed just before the simulation. + */ + void subscribeToSimulationSetupNotifications(MDModulesNotifiers* notifier) override + { + // TODO: add a check for threadmpi (see #5104, https://gitlab.com/gromacs/gromacs/-/merge_requests/4367#note_2102475958, the manual and the force provider for the details) + + // Access the plumed filename this is used to activate the plumed module + notifier->simulationSetupNotifier_.subscribe( + [this](const PlumedInputFilename& plumedFilename) + { this->options_.setPlumedFile(plumedFilename.plumedFilename_); }); + // Access the temperature if it is constant during the simulation + notifier->simulationSetupNotifier_.subscribe( + [this](const EnsembleTemperature& ensembleT) + { this->options_.setEnsembleTemperature(ensembleT); }); + // Access of the topology + notifier->simulationSetupNotifier_.subscribe([this](const gmx_mtop_t& mtop) + { this->options_.setTopology(mtop); }); + // Retrieve the Communication Record during simulations setup + notifier->simulationSetupNotifier_.subscribe([this](const t_commrec& cr) + { this->options_.setComm(cr); }); + // setting the simulation time step + notifier->simulationSetupNotifier_.subscribe( + [this](const SimulationTimeStep& simulationTimeStep) + { this->options_.setSimulationTimeStep(simulationTimeStep.delta_t); }); + // Retrieve the starting behavior + notifier->simulationSetupNotifier_.subscribe( + [this](const StartingBehavior& startingBehavior) + { this->options_.setStartingBehavior(startingBehavior); }); + // writing checkpoint data + notifier->checkpointingNotifier_.subscribe( + [this](MDModulesWriteCheckpointData /*checkpointData*/) + { + if (options_.active()) + { + plumedForceProvider_->writeCheckpointData(); + } + }); + } + + //! From IMDModule + IMdpOptionProvider* mdpOptionProvider() override { return nullptr; } + //! From IMDModule + IMDOutputProvider* outputProvider() override + { // Plumed provide its own output + return nullptr; + } + //! From IMDModule - Adds this module to the force providers if active + void initForceProviders(ForceProviders* forceProviders) override + { + if (options_.active()) + { + plumedForceProvider_ = std::make_unique(options_.options()); + forceProviders->addForceProvider(plumedForceProvider_.get(), "Plumed"); + } + } + + +private: + //! Parameters that become available at simulation setup time. + PlumedOptionProvider options_{}; + //! Object that evaluates the forces + std::unique_ptr plumedForceProvider_{}; +}; + +} // namespace + +std::unique_ptr PlumedModuleInfo::create() +{ + return std::make_unique(); +} +} // 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 new file mode 100644 index 0000000000..af644d9076 --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.cpp @@ -0,0 +1,104 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2024- 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 + * Defines options for Plumed. This class handles parameters set during + * pre-processing time. + * + * \author Daniele Rapetti + * \ingroup module_applied_forces + */ +#include "plumedOptions.h" + +#include "gromacs/math/vec.h" +#include "gromacs/mdrunutility/handlerestart.h" +#include "gromacs/mdrunutility/mdmodulesnotifiers.h" +#include "gromacs/topology/topology.h" + +namespace gmx +{ + +void PlumedOptionProvider::setTopology(const gmx_mtop_t& mtop) +{ + opts_.natoms_ = mtop.natoms; +} + +void PlumedOptionProvider::setEnsembleTemperature(const EnsembleTemperature& temp) +{ + if (temp.constantEnsembleTemperature_) + { + opts_.ensembleTemperature_ = temp.constantEnsembleTemperature_; + } +} + +void PlumedOptionProvider::setPlumedFile(const std::optional& fname) +{ + if (fname.has_value()) + { + opts_.active_ = true; + opts_.plumedFile_ = fname.value(); + } +} +const PlumedOptions& PlumedOptionProvider::options() const +{ + return opts_; +} + +void PlumedOptionProvider::setSimulationTimeStep(double timeStep) +{ + opts_.simulationTimeStep_ = timeStep; +} + +void PlumedOptionProvider::setStartingBehavior(const StartingBehavior& behavior) +{ + opts_.startingBehavior_ = behavior; +} + +void PlumedOptionProvider::setComm(const t_commrec& cr) +{ + opts_.cr_ = &cr; +} + +void PlumedOptionProvider::setMultisim(const gmx_multisim_t* ms) +{ + opts_.ms_ = ms; +} + +bool PlumedOptionProvider::active() const +{ + return opts_.active_; +} + + +} // namespace gmx diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.cpp.preplumed b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.cpp.preplumed new file mode 100644 index 0000000000..f3de891cc0 --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.cpp.preplumed @@ -0,0 +1,99 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2024- 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 + * Defines options for Plumed. This class handles parameters set during + * pre-processing time. + * + * \author Daniele Rapetti + * \ingroup module_applied_forces + */ +#include "plumedOptions.h" + +#include "gromacs/math/vec.h" +#include "gromacs/mdrunutility/handlerestart.h" +#include "gromacs/mdrunutility/mdmodulesnotifiers.h" +#include "gromacs/topology/topology.h" + +namespace gmx +{ + +void PlumedOptionProvider::setTopology(const gmx_mtop_t& mtop) +{ + opts_.natoms_ = mtop.natoms; +} + +void PlumedOptionProvider::setEnsembleTemperature(const EnsembleTemperature& temp) +{ + if (temp.constantEnsembleTemperature_) + { + opts_.ensembleTemperature_ = temp.constantEnsembleTemperature_; + } +} + +void PlumedOptionProvider::setPlumedFile(const std::optional& fname) +{ + if (fname.has_value()) + { + opts_.active_ = true; + opts_.plumedFile_ = fname.value(); + } +} +const PlumedOptions& PlumedOptionProvider::options() const +{ + return opts_; +} + +void PlumedOptionProvider::setSimulationTimeStep(double timeStep) +{ + opts_.simulationTimeStep_ = timeStep; +} + +void PlumedOptionProvider::setStartingBehavior(const StartingBehavior& behavior) +{ + opts_.startingBehavior_ = behavior; +} + +void PlumedOptionProvider::setComm(const t_commrec& cr) +{ + opts_.cr_ = &cr; +} + +bool PlumedOptionProvider::active() const +{ + return opts_.active_; +} + + +} // namespace gmx 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 new file mode 100644 index 0000000000..199e4e44a4 --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h @@ -0,0 +1,117 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2024- 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 + * Declares options for PLUMED. This class handles parameters set during + * pre-processing time. + * + * \author Daniele Rapetti + * \ingroup module_applied_forces + */ +#ifndef GMX_APPLIED_FORCES_PLUMEDOPTIONPROVIDER_H +#define GMX_APPLIED_FORCES_PLUMEDOPTIONPROVIDER_H + +#include +#include + +#include "gromacs/utility/real.h" + +struct gmx_mtop_t; +struct t_commrec; +struct gmx_multisim_t; + +namespace gmx +{ +enum class StartingBehavior; +struct EnsembleTemperature; + +struct PlumedOptions +{ + std::string plumedFile_; + int natoms_; + const t_commrec* cr_; + const gmx_multisim_t* ms_; + real simulationTimeStep_; + std::optional ensembleTemperature_{}; + StartingBehavior startingBehavior_{}; + bool active_{ false }; +}; + +class PlumedOptionProvider +{ +public: + /*! @brief Sets the needed informations from the topology object + * + * As now oly hte number of atoms is fetched + * @param mtop topology object + */ + void setTopology(const gmx_mtop_t& mtop); + /*! @brief Sets the (eventual) ensemble temperature + * @param temp the object with the optional temperature + */ + void setEnsembleTemperature(const EnsembleTemperature& temp); + /*! @brief Sets the name of the PLUMED file to read + * + * When called, with a non empty string, it activates the PLUMED module + * this simulation. + * @param fname the (optional) name of the file + */ + void setPlumedFile(const std::optional& fname); + /*! @brief Sets the timestep + * @param timeStep the timestep value + */ + void setSimulationTimeStep(double timeStep); + /*! @brief Sets the starting beahviour of the simulation + * @param startingBehavior the starting behaviopur object + */ + void setStartingBehavior(const StartingBehavior& startingBehavior); + /*! @brief Sets the address to the communication record object + * @param cr the Communication Record object + */ + void setComm(const t_commrec& cr); + /*! @brief Sets the address to the multisimulation object + * @param ms the address to the multisimulation object + */ + void setMultisim(const gmx_multisim_t* ms); + //! @brief returns the active status of the module + bool active() const; + + //! @brief returns a reference to the internal PlumedOptions element + const PlumedOptions& options() const; + +private: + PlumedOptions opts_; +}; +} // namespace gmx +#endif // GMX_APPLIED_FORCES_PLUMEDOPTIONPROVIDER_H diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h.preplumed b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h.preplumed new file mode 100644 index 0000000000..8451762f2b --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedOptions.h.preplumed @@ -0,0 +1,111 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2024- 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 + * Declares options for PLUMED. This class handles parameters set during + * pre-processing time. + * + * \author Daniele Rapetti + * \ingroup module_applied_forces + */ +#ifndef GMX_APPLIED_FORCES_PLUMEDOPTIONPROVIDER_H +#define GMX_APPLIED_FORCES_PLUMEDOPTIONPROVIDER_H + +#include +#include + +#include "gromacs/utility/real.h" + +struct gmx_mtop_t; +struct t_commrec; + +namespace gmx +{ +enum class StartingBehavior; +struct EnsembleTemperature; + +struct PlumedOptions +{ + std::string plumedFile_; + int natoms_; + const t_commrec* cr_; + real simulationTimeStep_; + std::optional ensembleTemperature_{}; + StartingBehavior startingBehavior_{}; + bool active_{ false }; +}; + +class PlumedOptionProvider +{ +public: + /*! @brief Sets the needed informations from the topology object + * + * As now oly hte number of atoms is fetched + * @param mtop topology object + */ + void setTopology(const gmx_mtop_t& mtop); + /*! @brief Sets the (eventual) ensemble temperature + * @param temp the object with the optional temperature + */ + void setEnsembleTemperature(const EnsembleTemperature& temp); + /*! @brief Sets the name of the PLUMED file to read + * + * When called, with a non empty string, it activates the PLUMED module + * this simulation. + * @param fname the (optional) name of the file + */ + void setPlumedFile(const std::optional& fname); + /*! @brief Sets the timestep + * @param timeStep the timestep value + */ + void setSimulationTimeStep(double timeStep); + /*! @brief Sets the starting beahviour of the simulation + * @param startingBehavior the starting behaviopur object + */ + void setStartingBehavior(const StartingBehavior& startingBehavior); + /*! @brief Sets the address to the communication record object + * @param cr the Communication Record object + */ + void setComm(const t_commrec& cr); + //! @brief returns the active status of the module + bool active() const; + + //! @brief returns a reference to the internal PlumedOptions element + const PlumedOptions& options() const; + +private: + PlumedOptions opts_; +}; +} // namespace gmx +#endif // GMX_APPLIED_FORCES_PLUMEDOPTIONPROVIDER_H 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 new file mode 100644 index 0000000000..fc9a00d1ca --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp @@ -0,0 +1,215 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2024- 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 + * Implementation of the Plumed force provider class + * + * \author Daniele Rapetti + * \ingroup module_applied_forces + */ + +#include "plumedforceprovider.h" + +#include "gromacs/domdec/domdec.h" +#include "gromacs/domdec/domdec_struct.h" +#include "gromacs/math/units.h" +#include "gromacs/mdrunutility/handlerestart.h" +#include "gromacs/mdrunutility/multisim.h" +#include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/enerdata.h" +#include "gromacs/mdtypes/forceoutput.h" +#include "gromacs/utility/exceptions.h" + +#include "plumedOptions.h" + +#define __PLUMED_WRAPPER_FORTRAN 0 // NOLINT(bugprone-reserved-identifier) + +#define __PLUMED_WRAPPER_LINK_RUNTIME 1 // NOLINT(bugprone-reserved-identifier) +#define __PLUMED_WRAPPER_EXTERN 0 // NOLINT(bugprone-reserved-identifier) + +#define __PLUMED_WRAPPER_CXX 1 // NOLINT(bugprone-reserved-identifier) +#define __PLUMED_WRAPPER_LIBCXX11 1 // NOLINT(bugprone-reserved-identifier) +#define __PLUMED_WRAPPER_LIBCXX17 1 // NOLINT(bugprone-reserved-identifier) +#define __PLUMED_WRAPPER_IMPLEMENTATION 1 // NOLINT(bugprone-reserved-identifier) +#define __PLUMED_HAS_DLOPEN // NOLINT(bugprone-reserved-identifier) + +#include "external/plumed/Plumed.h" + +namespace gmx +{ + +PlumedForceProvider::~PlumedForceProvider() = default; +PlumedForceProvider::PlumedForceProvider(const PlumedOptions& options) +try : plumed_(std::make_unique()) +{ + // I prefer to pass a struct with data because it stops the coupling + // at the implementation and not at the function signature: + // less code to edit when adding new options :) +#if GMX_THREAD_MPI + if (options.cr_->nnodes > 1) + { + GMX_THROW(InvalidInputError( + "plumed MPI interface is not compatible with THREAD_MPI when uses more than one " + "rank")); + } +#endif + int real_precision = sizeof(real); + plumed_->cmd("setRealPrecision", &real_precision); + real energyUnits = 1.0; + plumed_->cmd("setMDEnergyUnits", &energyUnits); + real lengthUnits = 1.0; + plumed_->cmd("setMDLengthUnits", &lengthUnits); + real timeUnits = 1.0; + plumed_->cmd("setMDTimeUnits", &timeUnits); + + plumed_->cmd("setPlumedDat", options.plumedFile_.c_str()); + + plumed_->cmd("getApiVersion", &plumedAPIversion_); + if (plumedAPIversion_ > 1) + { + /* setting kbT is only implemented with api>1) */ + if (options.ensembleTemperature_.has_value()) + { + real kbT = options.ensembleTemperature_.value() * gmx::c_boltz; + plumed_->cmd("setKbT", &kbT); + } + } + + if (plumedAPIversion_ > 2) + { + if ((options.startingBehavior_ != StartingBehavior::NewSimulation)) + { + int res = 1; + plumed_->cmd("setRestart", &res); + } + } + + if (isMultiSim(options.ms_)) + { + if (MAIN(options.cr_)) + plumed_->cmd("GREX setMPIIntercomm", &options.ms_->mainRanksComm_); + plumed_->cmd("GREX setMPIIntracomm", &options.cr_->mpi_comm_mygroup); + plumed_->cmd("GREX init", nullptr); + } + + if (PAR(options.cr_)) + { + if (havePPDomainDecomposition(options.cr_)) + { + plumed_->cmd("setMPIComm", &options.cr_->dd->mpi_comm_all); + } + } + plumed_->cmd("setNatoms", options.natoms_); + plumed_->cmd("setMDEngine", "gromacs"); + plumed_->cmd("setLogFile", "PLUMED.OUT"); + + plumed_->cmd("setTimestep", &options.simulationTimeStep_); + plumed_->cmd("init", nullptr); + + if (haveDDAtomOrdering(*options.cr_)) + { + int nat_home = dd_numHomeAtoms(*options.cr_->dd); + plumed_->cmd("setAtomsNlocal", &nat_home); + plumed_->cmd("setAtomsGatindex", options.cr_->dd->globalAtomIndices.data()); + } +} +catch (const std::exception& ex) +{ + GMX_THROW(InternalError( + std::string("An error occurred while initializing the PLUMED force provider:\n") + ex.what())); +} + +void PlumedForceProvider::writeCheckpointData() +try +{ + if (plumedAPIversion_ > 3) + { + int checkp = 1; + plumed_->cmd("doCheckPoint", &checkp); + } +} +catch (const std::exception& ex) +{ + GMX_THROW(InternalError( + std::string("An error occurred while PLUMED was writing the checkpoint data\n:") + ex.what())); +} + +void PlumedForceProvider::calculateForces(const ForceProviderInput& forceProviderInput, + ForceProviderOutput* forceProviderOutput) +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_; + plumed_->cmd("setStepLong", &lstep); + if (haveDDAtomOrdering(*cr)) + { + int nat_home = dd_numHomeAtoms(*cr->dd); + plumed_->cmd("setAtomsNlocal", &nat_home); + plumed_->cmd("setAtomsGatindex", cr->dd->globalAtomIndices.data()); + } + + plumed_->cmd("setPositions", &(forceProviderInput.x_.data()->as_vec()[0])); + plumed_->cmd("setMasses", forceProviderInput.massT_.data()); + plumed_->cmd("setCharges", forceProviderInput.chargeA_.data()); + plumed_->cmd("setBox", &forceProviderInput.box_[0][0]); + + plumed_->cmd("prepareCalc", nullptr); + + int plumedWantsToStop = 0; + plumed_->cmd("setStopFlag", &plumedWantsToStop); + + real* fOut = &(forceProviderOutput->forceWithVirial_.force_.data()->as_vec()[0]); + plumed_->cmd("setForces", fOut); + + matrix plumed_vir; + clear_mat(plumed_vir); + plumed_->cmd("setVirial", &plumed_vir[0][0]); + + // end setup: these instructions in the original patch are BEFORE do_force() + // in the original patch do_force() was called HERE + + // Do the work + plumed_->cmd("performCalc", nullptr); + + msmul(plumed_vir, 0.5, plumed_vir); + forceProviderOutput->forceWithVirial_.addVirialContribution(plumed_vir); +} +catch (const std::exception& ex) +{ + GMX_THROW(InternalError( + std::string("An error occurred while PLUMED was calculating the forces\n:") + ex.what())); +} +} // namespace gmx diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp.preplumed b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp.preplumed new file mode 100644 index 0000000000..6be72ac96c --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/plumedforceprovider.cpp.preplumed @@ -0,0 +1,206 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2024- 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 + * Implementation of the Plumed force provider class + * + * \author Daniele Rapetti + * \ingroup module_applied_forces + */ + +#include "plumedforceprovider.h" + +#include "gromacs/domdec/domdec.h" +#include "gromacs/domdec/domdec_struct.h" +#include "gromacs/math/units.h" +#include "gromacs/mdrunutility/handlerestart.h" +#include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/enerdata.h" +#include "gromacs/mdtypes/forceoutput.h" +#include "gromacs/utility/exceptions.h" + +#include "plumedOptions.h" + +#define __PLUMED_WRAPPER_FORTRAN 0 // NOLINT(bugprone-reserved-identifier) + +#define __PLUMED_WRAPPER_LINK_RUNTIME 1 // NOLINT(bugprone-reserved-identifier) +#define __PLUMED_WRAPPER_EXTERN 0 // NOLINT(bugprone-reserved-identifier) + +#define __PLUMED_WRAPPER_CXX 1 // NOLINT(bugprone-reserved-identifier) +#define __PLUMED_WRAPPER_LIBCXX11 1 // NOLINT(bugprone-reserved-identifier) +#define __PLUMED_WRAPPER_LIBCXX17 1 // NOLINT(bugprone-reserved-identifier) +#define __PLUMED_WRAPPER_IMPLEMENTATION 1 // NOLINT(bugprone-reserved-identifier) +#define __PLUMED_HAS_DLOPEN // NOLINT(bugprone-reserved-identifier) + +#include "external/plumed/Plumed.h" + +namespace gmx +{ + +PlumedForceProvider::~PlumedForceProvider() = default; +PlumedForceProvider::PlumedForceProvider(const PlumedOptions& options) +try : plumed_(std::make_unique()) +{ + // I prefer to pass a struct with data because it stops the coupling + // at the implementation and not at the function signature: + // less code to edit when adding new options :) +#if GMX_THREAD_MPI + if (options.cr_->nnodes > 1) + { + GMX_THROW(InvalidInputError( + "plumed MPI interface is not compatible with THREAD_MPI when uses more than one " + "rank")); + } +#endif + int real_precision = sizeof(real); + plumed_->cmd("setRealPrecision", &real_precision); + real energyUnits = 1.0; + plumed_->cmd("setMDEnergyUnits", &energyUnits); + real lengthUnits = 1.0; + plumed_->cmd("setMDLengthUnits", &lengthUnits); + real timeUnits = 1.0; + plumed_->cmd("setMDTimeUnits", &timeUnits); + + plumed_->cmd("setPlumedDat", options.plumedFile_.c_str()); + + plumed_->cmd("getApiVersion", &plumedAPIversion_); + if (plumedAPIversion_ > 1) + { + /* setting kbT is only implemented with api>1) */ + if (options.ensembleTemperature_.has_value()) + { + real kbT = options.ensembleTemperature_.value() * gmx::c_boltz; + plumed_->cmd("setKbT", &kbT); + } + } + + if (plumedAPIversion_ > 2) + { + if ((options.startingBehavior_ != StartingBehavior::NewSimulation)) + { + int res = 1; + plumed_->cmd("setRestart", &res); + } + } + + if (PAR(options.cr_)) + { + if (haveDDAtomOrdering(*options.cr_)) + { + plumed_->cmd("setMPIComm", &options.cr_->dd->mpi_comm_all); + } + } + plumed_->cmd("setNatoms", options.natoms_); + plumed_->cmd("setMDEngine", "gromacs"); + plumed_->cmd("setLogFile", "PLUMED.OUT"); + + plumed_->cmd("setTimestep", &options.simulationTimeStep_); + plumed_->cmd("init", nullptr); + + if (haveDDAtomOrdering(*options.cr_)) + { + int nat_home = dd_numHomeAtoms(*options.cr_->dd); + plumed_->cmd("setAtomsNlocal", &nat_home); + plumed_->cmd("setAtomsGatindex", options.cr_->dd->globalAtomIndices.data()); + } +} +catch (const std::exception& ex) +{ + GMX_THROW(InternalError( + std::string("An error occurred while initializing the PLUMED force provider:\n") + ex.what())); +} + +void PlumedForceProvider::writeCheckpointData() +try +{ + if (plumedAPIversion_ > 3) + { + int checkp = 1; + plumed_->cmd("doCheckPoint", &checkp); + } +} +catch (const std::exception& ex) +{ + GMX_THROW(InternalError( + std::string("An error occurred while PLUMED was writing the checkpoint data\n:") + ex.what())); +} + +void PlumedForceProvider::calculateForces(const ForceProviderInput& forceProviderInput, + ForceProviderOutput* forceProviderOutput) +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_; + plumed_->cmd("setStepLong", &lstep); + if (haveDDAtomOrdering(*cr)) + { + int nat_home = dd_numHomeAtoms(*cr->dd); + plumed_->cmd("setAtomsNlocal", &nat_home); + plumed_->cmd("setAtomsGatindex", cr->dd->globalAtomIndices.data()); + } + + plumed_->cmd("setPositions", &(forceProviderInput.x_.data()->as_vec()[0])); + plumed_->cmd("setMasses", forceProviderInput.massT_.data()); + plumed_->cmd("setCharges", forceProviderInput.chargeA_.data()); + plumed_->cmd("setBox", &forceProviderInput.box_[0][0]); + + plumed_->cmd("prepareCalc", nullptr); + + int plumedWantsToStop = 0; + plumed_->cmd("setStopFlag", &plumedWantsToStop); + + real* fOut = &(forceProviderOutput->forceWithVirial_.force_.data()->as_vec()[0]); + plumed_->cmd("setForces", fOut); + + matrix plumed_vir; + clear_mat(plumed_vir); + plumed_->cmd("setVirial", &plumed_vir[0][0]); + + // end setup: these instructions in the original patch are BEFORE do_force() + // in the original patch do_force() was called HERE + + // Do the work + plumed_->cmd("performCalc", nullptr); + + msmul(plumed_vir, 0.5, plumed_vir); + forceProviderOutput->forceWithVirial_.addVirialContribution(plumed_vir); +} +catch (const std::exception& ex) +{ + GMX_THROW(InternalError( + std::string("An error occurred while PLUMED was calculating the forces\n:") + ex.what())); +} +} // namespace gmx From 1c38aa22467fab682f76a8f6adad3f3402fb3801 Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Mon, 2 Dec 2024 14:45:57 +0100 Subject: [PATCH 2/7] adding the possibility to ignore the addition of the extra files --- patches/gromacs-2025.0.config | 2 ++ patches/patch.sh | 24 ++++++++++++++---------- 2 files changed, 16 insertions(+), 10 deletions(-) diff --git a/patches/gromacs-2025.0.config b/patches/gromacs-2025.0.config index b998dcb06d..9d3bce3bcb 100644 --- a/patches/gromacs-2025.0.config +++ b/patches/gromacs-2025.0.config @@ -32,3 +32,5 @@ EOF plumed_before_patch(){ plumed_patch_info } + +PLUMED_ONLY_PATCH=true diff --git a/patches/patch.sh b/patches/patch.sh index 4a454077d0..df77f44442 100755 --- a/patches/patch.sh +++ b/patches/patch.sh @@ -284,16 +284,20 @@ case "$action" in test -n "$quiet" || echo "Executing plumed_before_patch function" plumed_before_patch fi - if test -n "$include" ; then - test -n "$quiet" || echo "Including Plumed.h, Plumed.inc, and Plumed.cmake ($mode mode)" - echo "#include \"$PLUMED_INCLUDEDIR/$PLUMED_PROGRAM_NAME/wrapper/Plumed.h\"" > Plumed.h - echo "include $PLUMED_ROOT/src/lib/Plumed.inc.$mode" > Plumed.inc - echo "include($PLUMED_ROOT/src/lib/Plumed.cmake.$mode)" > Plumed.cmake - else - test -n "$quiet" || echo "Linking Plumed.h, Plumed.inc, and Plumed.cmake ($mode mode)" - ln -fs "$PLUMED_INCLUDEDIR/$PLUMED_PROGRAM_NAME/wrapper/Plumed.h" Plumed.h - ln -fs "$PLUMED_ROOT/src/lib/Plumed.inc.$mode" Plumed.inc - ln -fs "$PLUMED_ROOT/src/lib/Plumed.cmake.$mode" Plumed.cmake + + if [[ -z $PLUMED_ONLY_PATCH ]]; then + #PLUMED_ONLY_PATCH is special for programs that have an embedded plumed integration + if test -n "$include" ; then + test -n "$quiet" || echo "Including Plumed.h, Plumed.inc, and Plumed.cmake ($mode mode)" + echo "#include \"$PLUMED_INCLUDEDIR/$PLUMED_PROGRAM_NAME/wrapper/Plumed.h\"" > Plumed.h + echo "include $PLUMED_ROOT/src/lib/Plumed.inc.$mode" > Plumed.inc + echo "include($PLUMED_ROOT/src/lib/Plumed.cmake.$mode)" > Plumed.cmake + else + test -n "$quiet" || echo "Linking Plumed.h, Plumed.inc, and Plumed.cmake ($mode mode)" + ln -fs "$PLUMED_INCLUDEDIR/$PLUMED_PROGRAM_NAME/wrapper/Plumed.h" Plumed.h + ln -fs "$PLUMED_ROOT/src/lib/Plumed.inc.$mode" Plumed.inc + ln -fs "$PLUMED_ROOT/src/lib/Plumed.cmake.$mode" Plumed.cmake + fi fi if [ -d "$diff" ]; then From 57abcd11e93d7057f0419605e70c428ab62fcee9 Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Mon, 27 Jan 2025 10:38:33 +0100 Subject: [PATCH 3/7] typo --- patches/patch.sh | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/patches/patch.sh b/patches/patch.sh index df77f44442..0a7e4bfa50 100755 --- a/patches/patch.sh +++ b/patches/patch.sh @@ -33,7 +33,7 @@ Options: -d FILE, --diff FILE set the path to diff file (default: ROOT/patches/ENGINE.diff) (*) -q, --quiet - do not write loggin information; useful with -i to print just + do not write logging information; useful with -i to print just the patching information -I, --include use include files rather than symbolic links @@ -225,6 +225,7 @@ test -n "$quiet" || echo "diff file: $diff" if [ -f "$config" ] then test -n "$quiet" || echo "sourcing config file: $config" + # shellcheck disable=SC1090 source "$config" fi @@ -257,7 +258,7 @@ case "$action" in exit 1 fi fi - if [ -e Plumed.h -o -e Plumed.inc -o -e Plumed.cmake ] + if [ -e Plumed.h ] || [ -e Plumed.inc ] || [ -e Plumed.cmake ] then if ( type -t plumed_before_revert 1>/dev/null || type -t plumed_after_revert 1>/dev/null) && ( type -t plumed_before_patch 1>/dev/null || type -t plumed_after_patch 1>/dev/null) then From d0727272cacd6891ff36db232c297d379846527c Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Wed, 29 Jan 2025 16:59:03 +0100 Subject: [PATCH 4/7] used patch --save-originals to have a base --- .../src/external/plumed/PlumedInclude.h | 14 + .../external/plumed/PlumedInclude.h.preplumed | 0 .../src/external/plumed/PlumedKernel.cpp | 3 + .../plumed/PlumedKernel.cpp.preplumed | 0 .../applied_forces/plumed/CMakeLists.txt | 55 + .../plumed/CMakeLists.txt.preplumed | 52 + .../applied_forces/plumed/PlumedOutside.cpp | 32 + .../plumed/PlumedOutside.cpp.preplumed | 0 .../applied_forces/plumed/PlumedOutside.h | 39 + .../plumed/PlumedOutside.h.preplumed | 0 .../applied_forces/plumed/plumedMDModule.cpp | 13 +- .../applied_forces/plumed/plumedOptions.cpp | 8 +- .../applied_forces/plumed/plumedOptions.h | 7 +- .../plumed/plumedforceprovider.cpp | 39 +- .../src/gromacs/mdrun/replicaexchange.cpp | 1483 +++++++++++++++++ .../mdrun/replicaexchange.cpp.preplumed | 1424 ++++++++++++++++ patches/patch.sh | 13 +- 17 files changed, 3151 insertions(+), 31 deletions(-) create mode 100644 patches/gromacs-2025.0.diff/src/external/plumed/PlumedInclude.h create mode 100644 patches/gromacs-2025.0.diff/src/external/plumed/PlumedInclude.h.preplumed create mode 100644 patches/gromacs-2025.0.diff/src/external/plumed/PlumedKernel.cpp create mode 100644 patches/gromacs-2025.0.diff/src/external/plumed/PlumedKernel.cpp.preplumed create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/CMakeLists.txt create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/CMakeLists.txt.preplumed create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.cpp create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.cpp.preplumed create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.h create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.h.preplumed create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/mdrun/replicaexchange.cpp create mode 100644 patches/gromacs-2025.0.diff/src/gromacs/mdrun/replicaexchange.cpp.preplumed diff --git a/patches/gromacs-2025.0.diff/src/external/plumed/PlumedInclude.h b/patches/gromacs-2025.0.diff/src/external/plumed/PlumedInclude.h new file mode 100644 index 0000000000..caa5f66dab --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/external/plumed/PlumedInclude.h @@ -0,0 +1,14 @@ +#define __PLUMED_WRAPPER_FORTRAN 0 // NOLINT(bugprone-reserved-identifier) + +#define __PLUMED_WRAPPER_LINK_RUNTIME 1 // NOLINT(bugprone-reserved-identifier) +#define __PLUMED_WRAPPER_EXTERN 1 // NOLINT(bugprone-reserved-identifier) + +#define __PLUMED_WRAPPER_CXX 1 // NOLINT(bugprone-reserved-identifier) +#define __PLUMED_WRAPPER_LIBCXX11 1 // NOLINT(bugprone-reserved-identifier) +#define __PLUMED_WRAPPER_LIBCXX17 1 // NOLINT(bugprone-reserved-identifier) +#ifndef __PLUMED_WRAPPER_IMPLEMENTATION // NOLINT(bugprone-reserved-identifier) +# define __PLUMED_WRAPPER_IMPLEMENTATION 0 // NOLINT(bugprone-reserved-identifier) +#endif +#define __PLUMED_HAS_DLOPEN // NOLINT(bugprone-reserved-identifier) + +#include "external/plumed/Plumed.h" \ No newline at end of file diff --git a/patches/gromacs-2025.0.diff/src/external/plumed/PlumedInclude.h.preplumed b/patches/gromacs-2025.0.diff/src/external/plumed/PlumedInclude.h.preplumed new file mode 100644 index 0000000000..e69de29bb2 diff --git a/patches/gromacs-2025.0.diff/src/external/plumed/PlumedKernel.cpp b/patches/gromacs-2025.0.diff/src/external/plumed/PlumedKernel.cpp new file mode 100644 index 0000000000..e15b5fa82a --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/external/plumed/PlumedKernel.cpp @@ -0,0 +1,3 @@ +#define __PLUMED_WRAPPER_IMPLEMENTATION 1 // NOLINT(bugprone-reserved-identifier) + +#include "PlumedInclude.h" \ No newline at end of file diff --git a/patches/gromacs-2025.0.diff/src/external/plumed/PlumedKernel.cpp.preplumed b/patches/gromacs-2025.0.diff/src/external/plumed/PlumedKernel.cpp.preplumed new file mode 100644 index 0000000000..e69de29bb2 diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/CMakeLists.txt b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/CMakeLists.txt new file mode 100644 index 0000000000..9150ee89ca --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/CMakeLists.txt @@ -0,0 +1,55 @@ +# +# This file is part of the GROMACS molecular simulation package. +# +# Copyright 2024- 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. + +# If Plumed is not activated then compile just the stub implementation plumedMDModule_stub.cpp +# In other case compile regular plumedMDModule.cpp and tests +if(GMX_PLUMED_ACTIVE) + set(PLUMED_SOURCES + ../../../external/plumed/PlumedKernel.cpp + PlumedOutside.cpp + plumedforceprovider.cpp + plumedOptions.cpp + plumedMDModule.cpp + ) + gmx_add_libgromacs_sources( + ${PLUMED_SOURCES} + ) + + if (BUILD_TESTING) + add_subdirectory(tests) + endif() +else() +gmx_add_libgromacs_sources( + plumedMDModule_stub.cpp + ) +endif() diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/CMakeLists.txt.preplumed b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/CMakeLists.txt.preplumed new file mode 100644 index 0000000000..20e4ed371f --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/CMakeLists.txt.preplumed @@ -0,0 +1,52 @@ +# +# This file is part of the GROMACS molecular simulation package. +# +# Copyright 2024- 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. + +# If Plumed is not activated then compile just the stub implementation plumedMDModule_stub.cpp +# In other case compile regular plumedMDModule.cpp and tests +if(GMX_PLUMED_ACTIVE) + set(PLUMED_SOURCES plumedforceprovider.cpp + plumedOptions.cpp + plumedMDModule.cpp + ) + gmx_add_libgromacs_sources( + ${PLUMED_SOURCES} + ) + + if (BUILD_TESTING) + add_subdirectory(tests) + endif() +else() +gmx_add_libgromacs_sources( + plumedMDModule_stub.cpp + ) +endif() diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.cpp b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.cpp new file mode 100644 index 0000000000..301f5efb00 --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.cpp @@ -0,0 +1,32 @@ +#include "PlumedOutside.h" + +/*PLUMED*/ +namespace PLMD +{ +PlumedOutside::PlumedOutside() {} +PlumedOutside::~PlumedOutside() {} +void PlumedOutside::setPlumed(Plumed* plumed) +{ + plumed_ = plumed; + active_ = true; +} + +bool PlumedOutside::plumedHREX() +{ + return false; + // not yet integrated + // return plumedHREX_; +} + +bool PlumedOutside::active() +{ + return active_; +} + +PlumedOutside& plumedOutside() +{ + static PlumedOutside outside; + return outside; +} /* */ +} // namespace PLMD +/*ENDPLUMED*/ \ No newline at end of file diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.cpp.preplumed b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.cpp.preplumed new file mode 100644 index 0000000000..e69de29bb2 diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.h b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.h new file mode 100644 index 0000000000..572a16429a --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.h @@ -0,0 +1,39 @@ +#ifndef PLMD_OUTSIDE_H +#define PLMD_OUTSIDE_H +#include "external/plumed/PlumedInclude.h" + +#include "gromacs/utility/exceptions.h" +namespace PLMD +{ +class Plumed; +class PlumedOutside +{ + Plumed* plumed_{ nullptr }; + bool plumedHREX_{ false }; + bool active_{ false }; + +public: + PlumedOutside(); + ~PlumedOutside(); + void setPlumed(Plumed* plumed); + template + void cmd(Key&& key, T&& val) + try + { + plumed_->cmd(std::forward(key), std::forward(val)); + } + catch (const std::exception& ex) + { + GMX_THROW(gmx::InternalError( + std::string("An error occurred while running a command from the PLUMED patch:\n") + + ex.what())); + } + + bool plumedHREX(); + bool active(); +}; + +PlumedOutside& plumedOutside(); +} // namespace PLMD + +#endif // PLMD_OUTSIDE_H \ No newline at end of file diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.h.preplumed b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.h.preplumed new file mode 100644 index 0000000000..e69de29bb2 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 57ded02a22..c56c259a6a 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,9 +55,12 @@ #include "plumedOptions.h" #include "plumedforceprovider.h" +class gmx_multisim_t; namespace gmx { + + namespace { @@ -91,6 +94,12 @@ class PlumedMDModule final : public IMDModule notifier->simulationSetupNotifier_.subscribe( [this](const PlumedInputFilename& plumedFilename) { this->options_.setPlumedFile(plumedFilename.plumedFilename_); }); + + // Retrieve the Multisim options + notifier->simulationSetupNotifier_.subscribe( + [this](const gmx_multisim_t* ms) { this->options_.setMultisim(ms); }); + + // Access the temperature if it is constant during the simulation notifier->simulationSetupNotifier_.subscribe( [this](const EnsembleTemperature& ensembleT) @@ -109,10 +118,6 @@ class PlumedMDModule final : public IMDModule notifier->simulationSetupNotifier_.subscribe( [this](const StartingBehavior& startingBehavior) { this->options_.setStartingBehavior(startingBehavior); }); - // Retrieve the Multisim options - notifier->simulationSetupNotifier_.subscribe( - [this](const gmx_multisim_t* ms) - { this->options_.setMultisim(ms); }); // writing checkpoint data notifier->checkpointingNotifier_.subscribe( [this](MDModulesWriteCheckpointData /*checkpointData*/) 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 af644d9076..0bbe6bd3c9 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 @@ -90,14 +90,14 @@ void PlumedOptionProvider::setComm(const t_commrec& cr) opts_.cr_ = &cr; } -void PlumedOptionProvider::setMultisim(const gmx_multisim_t* ms) +bool PlumedOptionProvider::active() const { - opts_.ms_ = ms; + return opts_.active_; } -bool PlumedOptionProvider::active() const +void PlumedOptionProvider::setMultisim(const gmx_multisim_t* ms) { - return opts_.active_; + opts_.ms_ = ms; } 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 199e4e44a4..097aaa1bd5 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 @@ -48,8 +48,8 @@ #include "gromacs/utility/real.h" struct gmx_mtop_t; -struct t_commrec; struct gmx_multisim_t; +struct t_commrec; namespace gmx { @@ -61,7 +61,7 @@ struct PlumedOptions std::string plumedFile_; int natoms_; const t_commrec* cr_; - const gmx_multisim_t* ms_; + const gmx_multisim_t* ms_; real simulationTimeStep_; std::optional ensembleTemperature_{}; StartingBehavior startingBehavior_{}; @@ -99,11 +99,12 @@ class PlumedOptionProvider /*! @brief Sets the address to the communication record object * @param cr the Communication Record object */ - void setComm(const t_commrec& cr); + void setComm(const t_commrec& cr); /*! @brief Sets the address to the multisimulation object * @param ms the address to the multisimulation object */ void setMultisim(const gmx_multisim_t* ms); + //! @brief returns the active status of the module bool active() const; 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 fc9a00d1ca..6b5a0fbae8 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,31 +41,21 @@ #include "plumedforceprovider.h" +#include "external/plumed/PlumedInclude.h" + #include "gromacs/domdec/domdec.h" #include "gromacs/domdec/domdec_struct.h" #include "gromacs/math/units.h" +#include "gromacs/mdlib/gmx_omp_nthreads.h" #include "gromacs/mdrunutility/handlerestart.h" -#include "gromacs/mdrunutility/multisim.h" #include "gromacs/mdtypes/commrec.h" #include "gromacs/mdtypes/enerdata.h" #include "gromacs/mdtypes/forceoutput.h" +#include "gromacs/mdrunutility/multisim.h" #include "gromacs/utility/exceptions.h" +#include "PlumedOutside.h" #include "plumedOptions.h" - -#define __PLUMED_WRAPPER_FORTRAN 0 // NOLINT(bugprone-reserved-identifier) - -#define __PLUMED_WRAPPER_LINK_RUNTIME 1 // NOLINT(bugprone-reserved-identifier) -#define __PLUMED_WRAPPER_EXTERN 0 // NOLINT(bugprone-reserved-identifier) - -#define __PLUMED_WRAPPER_CXX 1 // NOLINT(bugprone-reserved-identifier) -#define __PLUMED_WRAPPER_LIBCXX11 1 // NOLINT(bugprone-reserved-identifier) -#define __PLUMED_WRAPPER_LIBCXX17 1 // NOLINT(bugprone-reserved-identifier) -#define __PLUMED_WRAPPER_IMPLEMENTATION 1 // NOLINT(bugprone-reserved-identifier) -#define __PLUMED_HAS_DLOPEN // NOLINT(bugprone-reserved-identifier) - -#include "external/plumed/Plumed.h" - namespace gmx { @@ -73,6 +63,7 @@ PlumedForceProvider::~PlumedForceProvider() = default; PlumedForceProvider::PlumedForceProvider(const PlumedOptions& options) try : plumed_(std::make_unique()) { + PLMD::plumedOutside().setPlumed(plumed_.get()); // I prefer to pass a struct with data because it stops the coupling // at the implementation and not at the function signature: // less code to edit when adding new options :) @@ -115,6 +106,21 @@ try : plumed_(std::make_unique()) } } + if (plumedAPIversion_ > 5) + { + int nth = gmx_omp_nthreads_get(ModuleMultiThread::Default); + plumed_->cmd("setNumOMPthreads", &nth); + } + + /* + if(plumedAPIversion_>9) { + //todo::get this info + plumed_->cmd("setGpuDeviceId", &deviceId); + } + */ + + + if (isMultiSim(options.ms_)) { if (MAIN(options.cr_)) @@ -123,9 +129,10 @@ try : plumed_(std::make_unique()) plumed_->cmd("GREX init", nullptr); } + if (PAR(options.cr_)) { - if (havePPDomainDecomposition(options.cr_)) + if (haveDDAtomOrdering(*options.cr_)) { plumed_->cmd("setMPIComm", &options.cr_->dd->mpi_comm_all); } diff --git a/patches/gromacs-2025.0.diff/src/gromacs/mdrun/replicaexchange.cpp b/patches/gromacs-2025.0.diff/src/gromacs/mdrun/replicaexchange.cpp new file mode 100644 index 0000000000..22cc844a00 --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/mdrun/replicaexchange.cpp @@ -0,0 +1,1483 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 1991- 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 replica exchange routines. + * + * \author David van der Spoel + * \author Mark Abraham + * \ingroup module_mdrun + */ +#include "gmxpre.h" + +#include "replicaexchange.h" + +#include "config.h" + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "gromacs/domdec/collect.h" +#include "gromacs/gmxlib/network.h" +#include "gromacs/math/functions.h" +#include "gromacs/math/units.h" +#include "gromacs/math/vec.h" +#include "gromacs/math/vectypes.h" +#include "gromacs/mdrunutility/multisim.h" +#include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/enerdata.h" +#include "gromacs/mdtypes/inputrec.h" +#include "gromacs/mdtypes/md_enums.h" +#include "gromacs/mdtypes/state.h" +#include "gromacs/random/seed.h" +#include "gromacs/random/threefry.h" +#include "gromacs/random/uniformintdistribution.h" +#include "gromacs/random/uniformrealdistribution.h" +#include "gromacs/topology/ifunc.h" +#include "gromacs/utility/arrayref.h" +#include "gromacs/utility/enumerationhelpers.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/gmxassert.h" +#include "gromacs/utility/pleasecite.h" +#include "gromacs/utility/smalloc.h" +/*PLUMED*/ +#include "gromacs/applied_forces/plumed/PlumedOutside.h" +/*ENDPLUMED*/ +//! Helps cut off probability values. +constexpr int c_probabilityCutoff = 100; + +/* we don't bother evaluating if events are more rare than exp(-100) = 3.7x10^-44 */ + +//! Rank in the multisimulation +#define MSRANK(ms, nodeid) (nodeid) + +//! Enum for replica exchange flavours +enum class ReplicaExchangeType : int +{ + Temperature, + Lambda, + EndSingle, + TemperatureLambda, + Count +}; +/*! \brief Strings describing replica exchange flavours. + * + * end_single_marker merely notes the end of single variable replica + * exchange. All types higher than it are multiple replica exchange + * methods. + * + * Eventually, should add 'pressure', 'temperature and pressure', + * 'lambda_and_pressure', 'temperature_lambda_pressure'?; Let's wait + * until we feel better about the pressure control methods giving + * exact ensembles. Right now, we assume constant pressure */ +static const char* enumValueToString(ReplicaExchangeType enumValue) +{ + constexpr gmx::EnumerationArray replicateExchangeTypeNames = { + "temperature", "lambda", "end_single_marker", "temperature and lambda" + }; + return replicateExchangeTypeNames[enumValue]; +} + +//! Working data for replica exchange. +struct gmx_repl_ex +{ + //! Replica ID + int repl; + //! Total number of replica + int nrepl; + //! Temperature + real temp; + //! Replica exchange type from ReplicaExchangeType enum + ReplicaExchangeType type; + //! Quantity, e.g. temperature or lambda; first index is ere, second index is replica ID + gmx::EnumerationArray q; + //! Use constant pressure and temperature + gmx_bool bNPT; + //! Replica pressures + real* pres; + //! Replica indices + int* ind; + //! Used for keeping track of all the replica swaps + int* allswaps; + //! Replica exchange interval (number of steps) + int nst; + //! Number of exchanges per interval + int nex; + //! Random seed + int seed; + //! Number of even and odd replica change attempts + int nattempt[2]; + //! Sum of probabilities + real* prob_sum; + //! Number of moves between replicas i and j + int** nmoves; + //! i-th element of the array is the number of exchanges between replica i-1 and i + int* nexchange; + + /*! \brief Helper arrays for replica exchange; allocated here + * so they don't have to be allocated each time */ + //! \{ + int* destinations; + int** cyclic; + int** order; + int* tmpswap; + gmx_bool* incycle; + gmx_bool* bEx; + //! \} + + //! Helper arrays to hold the quantities that are exchanged. + //! \{ + real* prob; + real* Epot; + real* beta; + real* Vol; + real** de; + //! \} +}; + +// TODO We should add Doxygen here some time. +//! \cond + +static gmx_bool repl_quantity(const gmx_multisim_t* ms, struct gmx_repl_ex* re, ReplicaExchangeType ere, real q) +{ + real* qall; + /*PLUMED*/ + constexpr gmx_bool bDiff = TRUE; + /*ENDPLUMED*/ + int s; + + snew(qall, ms->numSimulations_); + qall[re->repl] = q; + gmx_sum_sim(ms->numSimulations_, qall, ms); + + /*PLUMED*/ /*ENDPLUMED*/ + + if (bDiff) + { + /* Set the replica exchange type and quantities */ + re->type = ere; + + snew(re->q[ere], re->nrepl); + for (s = 0; s < ms->numSimulations_; s++) + { + re->q[ere][s] = qall[s]; + } + } + sfree(qall); + return bDiff; +} + +gmx_repl_ex_t init_replica_exchange(FILE* fplog, + const gmx_multisim_t* ms, + int numAtomsInSystem, + const t_inputrec* ir, + const ReplicaExchangeParameters& replExParams) +{ + real pres; + int i, j; + struct gmx_repl_ex* re; + gmx_bool bTemp; + gmx_bool bLambda = FALSE; + + fprintf(fplog, "\nInitializing Replica Exchange\n"); + + if (!isMultiSim(ms) || ms->numSimulations_ == 1) + { + gmx_fatal(FARGS, + "Nothing to exchange with only one replica, maybe you forgot to set the " + "-multidir option of mdrun?"); + } + if (replExParams.numExchanges < 0) + { + gmx_fatal(FARGS, "Replica exchange number of exchanges needs to be positive"); + } + + if (!EI_DYNAMICS(ir->eI)) + { + gmx_fatal(FARGS, "Replica exchange is only supported by dynamical simulations"); + /* Note that PAR(cr) is defined by cr->nnodes > 1, which is + * distinct from isMultiSim(ms). A multi-simulation only runs + * with real MPI parallelism, but this does not imply PAR(cr) + * is true! + * + * Since we are using a dynamical integrator, the only + * decomposition is DD, so PAR(cr) and haveDDAtomOrdering(*cr) are + * synonymous. The only way for cr->nnodes > 1 to be true is + * if we are using DD. */ + } + + if (!haveConstantEnsembleTemperature(*ir)) + { + gmx_fatal(FARGS, + "Replica exchange is only supported for systems that have a constant ensemble " + "temperature"); + } + + snew(re, 1); + + re->repl = ms->simulationIndex_; + re->nrepl = ms->numSimulations_; + + fprintf(fplog, "Repl There are %d replicas:\n", re->nrepl); + + /* We only check that the number of atoms in the systms match. + * This, of course, do not guarantee that the systems are the same, + * but it does guarantee that we can perform replica exchange. + */ + check_multi_int(fplog, ms, numAtomsInSystem, "the number of atoms", FALSE); + check_multi_int(fplog, ms, static_cast(ir->eI), "the integrator", FALSE); + check_multi_int64(fplog, ms, ir->init_step + ir->nsteps, "init_step+nsteps", FALSE); + const int nst = replExParams.exchangeInterval; + check_multi_int64(fplog, + ms, + gmx::divideRoundUp(ir->init_step, nst), + "first exchange step: init_step/-replex", + FALSE); + check_multi_int(fplog, ms, static_cast(ir->etc), "the temperature coupling", FALSE); + check_multi_int(fplog, ms, ir->opts.ngtc, "the number of temperature coupling groups", FALSE); + check_multi_int( + fplog, ms, static_cast(ir->pressureCouplingOptions.epc), "the pressure coupling", FALSE); + check_multi_int(fplog, ms, static_cast(ir->efep), "free energy", FALSE); + check_multi_int(fplog, ms, ir->fepvals->n_lambda, "number of lambda states", FALSE); + + re->temp = constantEnsembleTemperature(*ir); + for (i = 1; (i < ir->opts.ngtc); i++) + { + if (ir->opts.ref_t[i] != re->temp) + { + fprintf(fplog, + "\nWARNING: The temperatures of the different temperature coupling groups are " + "not identical\n\n"); + fprintf(stderr, + "\nWARNING: The temperatures of the different temperature coupling groups are " + "not identical\n\n"); + } + } + + re->type = ReplicaExchangeType::Count; + bTemp = repl_quantity(ms, re, ReplicaExchangeType::Temperature, re->temp); + if (ir->efep != FreeEnergyPerturbationType::No) + { + bLambda = repl_quantity( + ms, re, ReplicaExchangeType::Lambda, static_cast(ir->fepvals->init_fep_state)); + } + if (re->type == ReplicaExchangeType::Count) /* nothing was assigned */ + { + gmx_fatal(FARGS, + "The properties of the %d systems are all the same, there is nothing to exchange", + re->nrepl); + } + if (bLambda && bTemp) + { + re->type = ReplicaExchangeType::TemperatureLambda; + } + + if (bTemp) + { + please_cite(fplog, "Sugita1999a"); + if (ir->pressureCouplingOptions.epc != PressureCoupling::No) + { + re->bNPT = TRUE; + fprintf(fplog, "Repl Using Constant Pressure REMD.\n"); + please_cite(fplog, "Okabe2001a"); + } + if (ir->etc == TemperatureCoupling::Berendsen) + { + gmx_fatal(FARGS, + "REMD with the %s thermostat does not produce correct potential energy " + "distributions, consider using the %s thermostat instead", + enumValueToString(ir->etc), + enumValueToString(TemperatureCoupling::VRescale)); + } + } + if (bLambda) + { + if (ir->fepvals->delta_lambda != 0) /* check this? */ + { + gmx_fatal(FARGS, "delta_lambda is not zero"); + } + } + if (re->bNPT) + { + snew(re->pres, re->nrepl); + if (ir->pressureCouplingOptions.epct == PressureCouplingType::SurfaceTension) + { + pres = ir->pressureCouplingOptions.ref_p[ZZ][ZZ]; + } + else + { + pres = 0; + j = 0; + for (i = 0; i < DIM; i++) + { + if (ir->pressureCouplingOptions.compress[i][i] != 0) + { + pres += ir->pressureCouplingOptions.ref_p[i][i]; + j++; + } + } + pres /= j; + } + re->pres[re->repl] = pres; + gmx_sum_sim(re->nrepl, re->pres, ms); + } + + /* Make an index for increasing replica order */ + /* only makes sense if one or the other is varying, not both! + if both are varying, we trust the order the person gave. */ + snew(re->ind, re->nrepl); + for (i = 0; i < re->nrepl; i++) + { + re->ind[i] = i; + } + + /*PLUMED*/ /*ENDPLUMED*/ + + /* keep track of all the swaps, starting with the initial placement. */ + snew(re->allswaps, re->nrepl); + for (i = 0; i < re->nrepl; i++) + { + re->allswaps[i] = re->ind[i]; + } + + switch (re->type) + { + case ReplicaExchangeType::Temperature: + fprintf(fplog, "\nReplica exchange in temperature\n"); + for (i = 0; i < re->nrepl; i++) + { + fprintf(fplog, " %5.1f", re->q[re->type][re->ind[i]]); + } + fprintf(fplog, "\n"); + break; + case ReplicaExchangeType::Lambda: + fprintf(fplog, "\nReplica exchange in lambda\n"); + for (i = 0; i < re->nrepl; i++) + { + fprintf(fplog, " %3d", static_cast(re->q[re->type][re->ind[i]])); + } + fprintf(fplog, "\n"); + break; + case ReplicaExchangeType::TemperatureLambda: + fprintf(fplog, "\nReplica exchange in temperature and lambda state\n"); + for (i = 0; i < re->nrepl; i++) + { + fprintf(fplog, " %5.1f", re->q[ReplicaExchangeType::Temperature][re->ind[i]]); + } + fprintf(fplog, "\n"); + for (i = 0; i < re->nrepl; i++) + { + fprintf(fplog, " %5d", static_cast(re->q[ReplicaExchangeType::Lambda][re->ind[i]])); + } + fprintf(fplog, "\n"); + break; + default: gmx_incons("Unknown replica exchange quantity"); + } + if (re->bNPT) + { + fprintf(fplog, "\nRepl p"); + for (i = 0; i < re->nrepl; i++) + { + fprintf(fplog, " %5.2f", re->pres[re->ind[i]]); + } + + for (i = 0; i < re->nrepl; i++) + { + if ((i > 0) && (re->pres[re->ind[i]] < re->pres[re->ind[i - 1]])) + { + fprintf(fplog, + "\nWARNING: The reference pressures decrease with increasing " + "temperatures\n\n"); + fprintf(stderr, + "\nWARNING: The reference pressures decrease with increasing " + "temperatures\n\n"); + } + } + } + re->nst = nst; + if (replExParams.randomSeed == -1) + { + if (isMainSim(ms)) + { + re->seed = static_cast(gmx::makeRandomSeed()); + } + else + { + re->seed = 0; + } + gmx_sumi_sim(1, &(re->seed), ms); + } + else + { + re->seed = replExParams.randomSeed; + } + fprintf(fplog, "\nReplica exchange interval: %d\n", re->nst); + fprintf(fplog, "\nReplica random seed: %d\n", re->seed); + + re->nattempt[0] = 0; + re->nattempt[1] = 0; + + snew(re->prob_sum, re->nrepl); + snew(re->nexchange, re->nrepl); + snew(re->nmoves, re->nrepl); + for (i = 0; i < re->nrepl; i++) + { + snew(re->nmoves[i], re->nrepl); + } + fprintf(fplog, "Replica exchange information below: ex and x = exchange, pr = probability\n"); + + /* generate space for the helper functions so we don't have to snew each time */ + + snew(re->destinations, re->nrepl); + snew(re->incycle, re->nrepl); + snew(re->tmpswap, re->nrepl); + snew(re->cyclic, re->nrepl); + snew(re->order, re->nrepl); + for (i = 0; i < re->nrepl; i++) + { + snew(re->cyclic[i], re->nrepl + 1); + snew(re->order[i], re->nrepl); + } + /* allocate space for the functions storing the data for the replicas */ + /* not all of these arrays needed in all cases, but they don't take + up much space, since the max size is nrepl**2 */ + snew(re->prob, re->nrepl); + snew(re->bEx, re->nrepl); + snew(re->beta, re->nrepl); + snew(re->Vol, re->nrepl); + snew(re->Epot, re->nrepl); + snew(re->de, re->nrepl); + for (i = 0; i < re->nrepl; i++) + { + snew(re->de[i], re->nrepl); + } + re->nex = replExParams.numExchanges; + return re; +} + +static void exchange_reals(const gmx_multisim_t gmx_unused* ms, int gmx_unused b, real* v, int n) +{ + real* buf; + int i; + + if (v) + { + snew(buf, n); +#if GMX_MPI + /* + MPI_Sendrecv(v, n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0, + buf,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0, + ms->mainRanksComm_,MPI_STATUS_IGNORE); + */ + { + MPI_Request mpi_req; + + MPI_Isend(v, n * sizeof(real), MPI_BYTE, MSRANK(ms, b), 0, ms->mainRanksComm_, &mpi_req); + MPI_Recv(buf, n * sizeof(real), MPI_BYTE, MSRANK(ms, b), 0, ms->mainRanksComm_, MPI_STATUS_IGNORE); + MPI_Wait(&mpi_req, MPI_STATUS_IGNORE); + } +#endif + for (i = 0; i < n; i++) + { + v[i] = buf[i]; + } + sfree(buf); + } +} + + +static void exchange_doubles(const gmx_multisim_t gmx_unused* ms, int gmx_unused b, double* v, int n) +{ + double* buf; + int i; + + if (v) + { + snew(buf, n); +#if GMX_MPI + /* + MPI_Sendrecv(v, n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0, + buf,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0, + ms->mainRanksComm_,MPI_STATUS_IGNORE); + */ + { + MPI_Request mpi_req; + + MPI_Isend(v, n * sizeof(double), MPI_BYTE, MSRANK(ms, b), 0, ms->mainRanksComm_, &mpi_req); + MPI_Recv(buf, n * sizeof(double), MPI_BYTE, MSRANK(ms, b), 0, ms->mainRanksComm_, MPI_STATUS_IGNORE); + MPI_Wait(&mpi_req, MPI_STATUS_IGNORE); + } +#endif + for (i = 0; i < n; i++) + { + v[i] = buf[i]; + } + sfree(buf); + } +} + +static void exchange_rvecs(const gmx_multisim_t gmx_unused* ms, int gmx_unused b, rvec* v, int n) +{ + rvec* buf; + int i; + + if (v) + { + snew(buf, n); +#if GMX_MPI + /* + MPI_Sendrecv(v[0], n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0, + buf[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0, + ms->mainRanksComm_,MPI_STATUS_IGNORE); + */ + { + MPI_Request mpi_req; + + MPI_Isend(v[0], n * sizeof(rvec), MPI_BYTE, MSRANK(ms, b), 0, ms->mainRanksComm_, &mpi_req); + MPI_Recv(buf[0], n * sizeof(rvec), MPI_BYTE, MSRANK(ms, b), 0, ms->mainRanksComm_, MPI_STATUS_IGNORE); + MPI_Wait(&mpi_req, MPI_STATUS_IGNORE); + } +#endif + for (i = 0; i < n; i++) + { + copy_rvec(buf[i], v[i]); + } + sfree(buf); + } +} + +static void exchange_state(const gmx_multisim_t* ms, int b, t_state* state) +{ + /* When t_state changes, this code should be updated. */ + int ngtc, nnhpres; + ngtc = state->ngtc * state->nhchainlength; + nnhpres = state->nnhpres * state->nhchainlength; + exchange_rvecs(ms, b, state->box, DIM); + exchange_rvecs(ms, b, state->box_rel, DIM); + exchange_rvecs(ms, b, state->boxv, DIM); + exchange_reals(ms, b, &(state->veta), 1); + exchange_reals(ms, b, &(state->vol0), 1); + exchange_rvecs(ms, b, state->svir_prev, DIM); + exchange_rvecs(ms, b, state->fvir_prev, DIM); + exchange_rvecs(ms, b, state->pres_prev, DIM); + exchange_doubles(ms, b, state->nosehoover_xi.data(), ngtc); + exchange_doubles(ms, b, state->nosehoover_vxi.data(), ngtc); + exchange_doubles(ms, b, state->nhpres_xi.data(), nnhpres); + exchange_doubles(ms, b, state->nhpres_vxi.data(), nnhpres); + exchange_doubles(ms, b, state->therm_integral.data(), state->ngtc); + exchange_doubles(ms, b, &state->baros_integral, 1); + exchange_rvecs(ms, b, state->x.rvec_array(), state->numAtoms()); + exchange_rvecs(ms, b, state->v.rvec_array(), state->numAtoms()); +} + +static void copy_state_serial(const t_state* src, t_state* dest) +{ + if (dest != src) + { + /* Currently the local state is always a pointer to the global + * in serial, so we should never end up here. + * TODO: Implement a (trivial) t_state copy once converted to C++. + */ + GMX_RELEASE_ASSERT(false, "State copying is currently not implemented in replica exchange"); + } +} + +static void scale_velocities(gmx::ArrayRef velocities, real fac) +{ + for (auto& v : velocities) + { + v *= fac; + } +} + +static void print_transition_matrix(FILE* fplog, int n, int** nmoves, const int* nattempt) +{ + int i, j, ntot; + float Tprint; + + ntot = nattempt[0] + nattempt[1]; + fprintf(fplog, "\n"); + fprintf(fplog, "Repl"); + for (i = 0; i < n; i++) + { + fprintf(fplog, " "); /* put the title closer to the center */ + } + fprintf(fplog, "Empirical Transition Matrix\n"); + + fprintf(fplog, "Repl"); + for (i = 0; i < n; i++) + { + fprintf(fplog, "%8d", (i + 1)); + } + fprintf(fplog, "\n"); + + for (i = 0; i < n; i++) + { + fprintf(fplog, "Repl"); + for (j = 0; j < n; j++) + { + Tprint = 0.0; + if (nmoves[i][j] > 0) + { + Tprint = nmoves[i][j] / (2.0 * ntot); + } + fprintf(fplog, "%8.4f", Tprint); + } + fprintf(fplog, "%3d\n", i); + } +} + +static void print_ind(FILE* fplog, const char* leg, int n, int* ind, const gmx_bool* bEx) +{ + int i; + + fprintf(fplog, "Repl %2s %2d", leg, ind[0]); + for (i = 1; i < n; i++) + { + fprintf(fplog, " %c %2d", (bEx != nullptr && bEx[i]) ? 'x' : ' ', ind[i]); + } + fprintf(fplog, "\n"); +} + +static void print_allswitchind(FILE* fplog, int n, int* pind, int* allswaps, int* tmpswap) +{ + int i; + + for (i = 0; i < n; i++) + { + tmpswap[i] = allswaps[i]; + } + for (i = 0; i < n; i++) + { + allswaps[i] = tmpswap[pind[i]]; + } + + fprintf(fplog, "\nAccepted Exchanges: "); + for (i = 0; i < n; i++) + { + fprintf(fplog, "%d ", pind[i]); + } + fprintf(fplog, "\n"); + + /* the "Order After Exchange" is the state label corresponding to the configuration that + started in state listed in order, i.e. + + 3 0 1 2 + + means that the: + configuration starting in simulation 3 is now in simulation 0, + configuration starting in simulation 0 is now in simulation 1, + configuration starting in simulation 1 is now in simulation 2, + configuration starting in simulation 2 is now in simulation 3 + */ + fprintf(fplog, "Order After Exchange: "); + for (i = 0; i < n; i++) + { + fprintf(fplog, "%d ", allswaps[i]); + } + fprintf(fplog, "\n\n"); +} + +static void print_prob(FILE* fplog, const char* leg, int n, real* prob) +{ + int i; + char buf[8]; + + fprintf(fplog, "Repl %2s ", leg); + for (i = 1; i < n; i++) + { + if (prob[i] >= 0) + { + sprintf(buf, "%4.2f", prob[i]); + fprintf(fplog, " %3s", buf[0] == '1' ? "1.0" : buf + 1); + } + else + { + fprintf(fplog, " "); + } + } + fprintf(fplog, "\n"); +} + +static void print_count(FILE* fplog, const char* leg, int n, int* count) +{ + int i; + + fprintf(fplog, "Repl %2s ", leg); + for (i = 1; i < n; i++) + { + fprintf(fplog, " %4d", count[i]); + } + fprintf(fplog, "\n"); +} + +static real calc_delta(FILE* fplog, gmx_bool bPrint, struct gmx_repl_ex* re, int a, int b, int ap, int bp) +{ + + real ediff, dpV, delta = 0; + real* Epot = re->Epot; + real* Vol = re->Vol; + real** de = re->de; + real* beta = re->beta; + + /* Two cases; we are permuted and not. In all cases, setting ap = a and bp = b will reduce + to the non permuted case */ + + switch (re->type) + { + case ReplicaExchangeType::Temperature: + /* + * Okabe et. al. Chem. Phys. Lett. 335 (2001) 435-439 + */ + ediff = Epot[b] - Epot[a]; + delta = -(beta[bp] - beta[ap]) * ediff; + break; + case ReplicaExchangeType::Lambda: + /* two cases: when we are permuted, and not. */ + /* non-permuted: + ediff = E_new - E_old + = [H_b(x_a) + H_a(x_b)] - [H_b(x_b) + H_a(x_a)] + = [H_b(x_a) - H_a(x_a)] + [H_a(x_b) - H_b(x_b)] + = de[b][a] + de[a][b] */ + + /* permuted: + ediff = E_new - E_old + = [H_bp(x_a) + H_ap(x_b)] - [H_bp(x_b) + H_ap(x_a)] + = [H_bp(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_bp(x_b)] + = [H_bp(x_a) - H_a(x_a) + H_a(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_b(x_b) + H_b(x_b) - H_bp(x_b)] + = [H_bp(x_a) - H_a(x_a)] - [H_ap(x_a) - H_a(x_a)] + [H_ap(x_b) - H_b(x_b)] - H_bp(x_b) - H_b(x_b)] + = (de[bp][a] - de[ap][a]) + (de[ap][b] - de[bp][b]) */ + /* but, in the current code implementation, we flip configurations, not indices . . . + So let's examine that. + = [H_b(x_ap) - H_a(x_a)] - [H_a(x_ap) - H_a(x_a)] + [H_a(x_bp) - H_b(x_b)] - H_b(x_bp) - H_b(x_b)] + = [H_b(x_ap) - H_a(x_ap)] + [H_a(x_bp) - H_b(x_pb)] + = (de[b][ap] - de[a][ap]) + (de[a][bp] - de[b][bp] + So, if we exchange b<=> bp and a<=> ap, we return to the same result. + So the simple solution is to flip the + position of perturbed and original indices in the tests. + */ + + ediff = (de[bp][a] - de[ap][a]) + (de[ap][b] - de[bp][b]); + delta = ediff * beta[a]; /* assume all same temperature in this case */ + break; + case ReplicaExchangeType::TemperatureLambda: + /* not permuted: */ + /* delta = reduced E_new - reduced E_old + = [beta_b H_b(x_a) + beta_a H_a(x_b)] - [beta_b H_b(x_b) + beta_a H_a(x_a)] + = [beta_b H_b(x_a) - beta_a H_a(x_a)] + [beta_a H_a(x_b) - beta_b H_b(x_b)] + = [beta_b dH_b(x_a) + beta_b H_a(x_a) - beta_a H_a(x_a)] + + [beta_a dH_a(x_b) + beta_a H_b(x_b) - beta_b H_b(x_b)] + = [beta_b dH_b(x_a) + [beta_a dH_a(x_b) + + beta_b (H_a(x_a) - H_b(x_b)]) - beta_a (H_a(x_a) - H_b(x_b)) + = beta_b dH_b(x_a) + beta_a dH_a(x_b) - (beta_b - beta_a)(H_b(x_b) - H_a(x_a) */ + /* delta = beta[b]*de[b][a] + beta[a]*de[a][b] - (beta[b] - beta[a])*(Epot[b] - Epot[a]; */ + /* permuted (big breath!) */ + /* delta = reduced E_new - reduced E_old + = [beta_bp H_bp(x_a) + beta_ap H_ap(x_b)] - [beta_bp H_bp(x_b) + beta_ap H_ap(x_a)] + = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)] + = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)] + - beta_pb H_a(x_a) + beta_ap H_a(x_a) + beta_pb H_a(x_a) - beta_ap H_a(x_a) + - beta_ap H_b(x_b) + beta_bp H_b(x_b) + beta_ap H_b(x_b) - beta_bp H_b(x_b) + = [(beta_bp H_bp(x_a) - beta_bp H_a(x_a)) - (beta_ap H_ap(x_a) - beta_ap H_a(x_a))] + + [(beta_ap H_ap(x_b) - beta_ap H_b(x_b)) - (beta_bp H_bp(x_b) - beta_bp H_b(x_b))] + + beta_pb H_a(x_a) - beta_ap H_a(x_a) + beta_ap H_b(x_b) - beta_bp H_b(x_b) + = [beta_bp (H_bp(x_a) - H_a(x_a)) - beta_ap (H_ap(x_a) - H_a(x_a))] + + [beta_ap (H_ap(x_b) - H_b(x_b)) - beta_bp (H_bp(x_b) - H_b(x_b))] + + beta_pb (H_a(x_a) - H_b(x_b)) - beta_ap (H_a(x_a) - H_b(x_b)) + = ([beta_bp de[bp][a] - beta_ap de[ap][a]) + beta_ap de[ap][b] - beta_bp de[bp][b]) + + (beta_pb-beta_ap)(H_a(x_a) - H_b(x_b)) */ + delta = beta[bp] * (de[bp][a] - de[bp][b]) + beta[ap] * (de[ap][b] - de[ap][a]) + - (beta[bp] - beta[ap]) * (Epot[b] - Epot[a]); + break; + default: gmx_incons("Unknown replica exchange quantity"); + } + if (bPrint) + { + fprintf(fplog, "Repl %d <-> %d dE_term = %10.3e (kT)\n", a, b, delta); + } + /*PLUMED*/ + if (PLMD::plumedOutside().plumedHREX()) + { + /* this is necessary because with plumed HREX the energy contribution is + already taken into account */ + delta = 0.0; + } + /*ENDPLUMED*/ + if (re->bNPT) + { + /* revist the calculation for 5.0. Might be some improvements. */ + dpV = (beta[ap] * re->pres[ap] - beta[bp] * re->pres[bp]) * (Vol[b] - Vol[a]) / gmx::c_presfac; + if (bPrint) + { + fprintf(fplog, " dpV = %10.3e d = %10.3e\n", dpV, delta + dpV); + } + delta += dpV; + } + return delta; +} + +static void test_for_replica_exchange(FILE* fplog, + const gmx_multisim_t* ms, + struct gmx_repl_ex* re, + const gmx_enerdata_t* enerd, + real vol, + int64_t step, + real time) +{ + int m, i, j, a, b, ap, bp, i0, i1, tmp; + real delta = 0; + gmx_bool bPrint, bMultiEx; + gmx_bool* bEx = re->bEx; + real* prob = re->prob; + int* pind = re->destinations; /* permuted index */ + gmx_bool bEpot = FALSE; + gmx_bool bDLambda = FALSE; + gmx_bool bVol = FALSE; + gmx::ThreeFry2x64<64> rng(re->seed, gmx::RandomDomain::ReplicaExchange); + gmx::UniformRealDistribution uniformRealDist; + gmx::UniformIntDistribution uniformNreplDist(0, re->nrepl - 1); + + bMultiEx = (re->nex > 1); /* multiple exchanges at each state */ + fprintf(fplog, "Replica exchange at step %" PRId64 " time %.5f\n", step, time); + + if (re->bNPT) + { + for (i = 0; i < re->nrepl; i++) + { + re->Vol[i] = 0; + } + bVol = TRUE; + re->Vol[re->repl] = vol; + } + if ((re->type == ReplicaExchangeType::Temperature || re->type == ReplicaExchangeType::TemperatureLambda)) + { + for (i = 0; i < re->nrepl; i++) + { + re->Epot[i] = 0; + } + bEpot = TRUE; + re->Epot[re->repl] = enerd->term[F_EPOT]; + /* temperatures of different states*/ + for (i = 0; i < re->nrepl; i++) + { + re->beta[i] = 1.0 / (re->q[ReplicaExchangeType::Temperature][i] * gmx::c_boltz); + } + } + else + { + for (i = 0; i < re->nrepl; i++) + { + re->beta[i] = 1.0 / (re->temp * gmx::c_boltz); /* we have a single temperature */ + } + } + if (re->type == ReplicaExchangeType::Lambda || re->type == ReplicaExchangeType::TemperatureLambda) + { + bDLambda = TRUE; + /* lambda differences. */ + /* de[i][j] is the energy of the jth simulation in the ith Hamiltonian + minus the energy of the jth simulation in the jth Hamiltonian */ + for (i = 0; i < re->nrepl; i++) + { + for (j = 0; j < re->nrepl; j++) + { + re->de[i][j] = 0; + } + } + for (i = 0; i < re->nrepl; i++) + { + re->de[i][re->repl] = + enerd->foreignLambdaTerms.deltaH(re->q[ReplicaExchangeType::Lambda][i]); + } + } + + /* now actually do the communication */ + if (bVol) + { + gmx_sum_sim(re->nrepl, re->Vol, ms); + } + if (bEpot) + { + gmx_sum_sim(re->nrepl, re->Epot, ms); + } + if (bDLambda) + { + for (i = 0; i < re->nrepl; i++) + { + gmx_sum_sim(re->nrepl, re->de[i], ms); + } + } + + /* make a duplicate set of indices for shuffling */ + for (i = 0; i < re->nrepl; i++) + { + pind[i] = re->ind[i]; + } + + rng.restart(step, 0); + + /* PLUMED */ + int plumed_test_exchange_pattern = 0; + if (plumed_test_exchange_pattern && PLMD::PlumedOutside().plumedHREX()) + gmx_fatal(FARGS, "hrex not compatible with ad hoc exchange patterns"); + /* END PLUMED */ + + if (bMultiEx) + { + /* multiple random switch exchange */ + int nself = 0; + + + for (i = 0; i < re->nex + nself; i++) + { + // For now this is superfluous, but just in case we ever add more + // calls in different branches it is safer to always reset the distribution. + uniformNreplDist.reset(); + + /* randomly select a pair */ + /* in theory, could reduce this by identifying only which switches had a nonneglibible + probability of occurring (log p > -100) and only operate on those switches */ + /* find out which state it is from, and what label that state currently has. Likely + more work that useful. */ + i0 = uniformNreplDist(rng); + i1 = uniformNreplDist(rng); + if (i0 == i1) + { + nself++; + continue; /* self-exchange, back up and do it again */ + } + + a = re->ind[i0]; /* what are the indices of these states? */ + b = re->ind[i1]; + ap = pind[i0]; + bp = pind[i1]; + + bPrint = FALSE; /* too noisy */ + /* calculate the energy difference */ + /* if the code changes to flip the STATES, rather than the configurations, + use the commented version of the code */ + /* delta = calc_delta(fplog,bPrint,re,a,b,ap,bp); */ + delta = calc_delta(fplog, bPrint, re, ap, bp, a, b); + + /* we actually only use the first space in the prob and bEx array, + since there are actually many switches between pairs. */ + + if (delta <= 0) + { + /* accepted */ + prob[0] = 1; + bEx[0] = TRUE; + } + else + { + if (delta > c_probabilityCutoff) + { + prob[0] = 0; + } + else + { + prob[0] = std::exp(-delta); + } + // roll a number to determine if accepted. For now it is superfluous to + // reset, but just in case we ever add more calls in different branches + // it is safer to always reset the distribution. + uniformRealDist.reset(); + bEx[0] = uniformRealDist(rng) < prob[0]; + } + re->prob_sum[0] += prob[0]; + + if (bEx[0]) + { + /* swap the states */ + tmp = pind[i0]; + pind[i0] = pind[i1]; + pind[i1] = tmp; + } + } + re->nattempt[0]++; /* keep track of total permutation trials here */ + print_allswitchind(fplog, re->nrepl, pind, re->allswaps, re->tmpswap); + } + else + { + /* standard nearest neighbor replica exchange */ + + /* PLUMED */ + if (PLMD::PlumedOutside().active()) + { + int partner = re->repl; + PLMD::PlumedOutside().cmd("getExchangesFlag", &plumed_test_exchange_pattern); + if (plumed_test_exchange_pattern > 0) + { + int* list; + snew(list, re->nrepl); + PLMD::PlumedOutside().cmd("setNumberOfReplicas", &(re->nrepl)); + PLMD::PlumedOutside().cmd("getExchangesList", list); + for (i = 0; i < re->nrepl; i++) + re->ind[i] = list[i]; + sfree(list); + } + + for (i = 1; i < re->nrepl; i++) + { + if (i % 2 != m) + continue; + a = re->ind[i - 1]; + b = re->ind[i]; + if (re->repl == a) + partner = b; + if (re->repl == b) + partner = a; + } + PLMD::PlumedOutside().cmd("GREX setPartner", &partner); + PLMD::PlumedOutside().cmd("GREX calculate", nullptr); + PLMD::PlumedOutside().cmd("GREX shareAllDeltaBias", nullptr); + } + /* END PLUMED */ + + m = (step / re->nst) % 2; + for (i = 1; i < re->nrepl; i++) + { + a = re->ind[i - 1]; + b = re->ind[i]; + + bPrint = (re->repl == a || re->repl == b); + if (i % 2 == m) + { + delta = calc_delta(fplog, bPrint, re, a, b, a, b); + /* PLUMED */ + if (PLMD::PlumedOutside().active()) + { + real adb, bdb, dplumed; + char buf[300]; + sprintf(buf, "GREX getDeltaBias %d", a); + PLMD::PlumedOutside().cmd(buf, &adb); + sprintf(buf, "GREX getDeltaBias %d", b); + PLMD::PlumedOutside().cmd(buf, &bdb); + dplumed = adb * re->beta[a] + bdb * re->beta[b]; + delta += dplumed; + if (bPrint) + fprintf(fplog, "dplumed = %10.3e dE_Term = %10.3e (kT)\n", dplumed, delta); + } + /* END PLUMED */ + if (delta <= 0) + { + /* accepted */ + prob[i] = 1; + bEx[i] = TRUE; + } + else + { + if (delta > c_probabilityCutoff) + { + prob[i] = 0; + } + else + { + prob[i] = std::exp(-delta); + } + // roll a number to determine if accepted. For now it is superfluous to + // reset, but just in case we ever add more calls in different branches + // it is safer to always reset the distribution. + uniformRealDist.reset(); + bEx[i] = uniformRealDist(rng) < prob[i]; + } + re->prob_sum[i] += prob[i]; + + if (bEx[i]) + { + /* PLUMED */ + if (!plumed_test_exchange_pattern) + { + /* standard neighbour swapping */ + /* swap these two */ + tmp = pind[i - 1]; + pind[i - 1] = pind[i]; + pind[i] = tmp; + re->nexchange[i]++; /* statistics for back compatibility */ + } + else + { + /* alternative swapping patterns */ + tmp = pind[a]; + pind[a] = pind[b]; + pind[b] = tmp; + re->nexchange[i]++; /* statistics for back compatibility */ + } + /* END PLUMED */ + } + } + else + { + prob[i] = -1; + bEx[i] = FALSE; + } + } + /* print some statistics */ + print_ind(fplog, "ex", re->nrepl, re->ind, bEx); + print_prob(fplog, "pr", re->nrepl, prob); + fprintf(fplog, "\n"); + re->nattempt[m]++; + } + + /* PLUMED */ + if (plumed_test_exchange_pattern > 0) + { + for (i = 0; i < re->nrepl; i++) + { + re->ind[i] = i; + } + } + /* END PLUMED */ + + /* record which moves were made and accepted */ + for (i = 0; i < re->nrepl; i++) + { + re->nmoves[re->ind[i]][pind[i]] += 1; + re->nmoves[pind[i]][re->ind[i]] += 1; + } + fflush(fplog); /* make sure we can see what the last exchange was */ +} + +static void cyclic_decomposition(const int* destinations, int** cyclic, gmx_bool* incycle, const int nrepl, int* nswap) +{ + + int i, j, c, p; + int maxlen = 1; + for (i = 0; i < nrepl; i++) + { + incycle[i] = FALSE; + } + for (i = 0; i < nrepl; i++) /* one cycle for each replica */ + { + if (incycle[i]) + { + cyclic[i][0] = -1; + continue; + } + cyclic[i][0] = i; + incycle[i] = TRUE; + c = 1; + p = i; + for (j = 0; j < nrepl; j++) /* potentially all cycles are part, but we will break first */ + { + p = destinations[p]; /* start permuting */ + if (p == i) + { + cyclic[i][c] = -1; + if (c > maxlen) + { + maxlen = c; + } + break; /* we've reached the original element, the cycle is complete, and we marked the end. */ + } + else + { + cyclic[i][c] = p; /* each permutation gives a new member of the cycle */ + incycle[p] = TRUE; + c++; + } + } + } + *nswap = maxlen - 1; + + if (debug) + { + for (i = 0; i < nrepl; i++) + { + fprintf(debug, "Cycle %d:", i); + for (j = 0; j < nrepl; j++) + { + if (cyclic[i][j] < 0) + { + break; + } + fprintf(debug, "%2d", cyclic[i][j]); + } + fprintf(debug, "\n"); + } + fflush(debug); + } +} + +static void compute_exchange_order(int** cyclic, int** order, const int nrepl, const int maxswap) +{ + int i, j; + + for (j = 0; j < maxswap; j++) + { + for (i = 0; i < nrepl; i++) + { + if (cyclic[i][j + 1] >= 0) + { + order[cyclic[i][j + 1]][j] = cyclic[i][j]; + order[cyclic[i][j]][j] = cyclic[i][j + 1]; + } + } + for (i = 0; i < nrepl; i++) + { + if (order[i][j] < 0) + { + order[i][j] = i; /* if it's not exchanging, it should stay this round*/ + } + } + } + + if (debug) + { + fprintf(debug, "Replica Exchange Order\n"); + for (i = 0; i < nrepl; i++) + { + fprintf(debug, "Replica %d:", i); + for (j = 0; j < maxswap; j++) + { + if (order[i][j] < 0) + { + break; + } + fprintf(debug, "%2d", order[i][j]); + } + fprintf(debug, "\n"); + } + fflush(debug); + } +} + +static void prepare_to_do_exchange(struct gmx_repl_ex* re, const int replica_id, int* maxswap, gmx_bool* bThisReplicaExchanged) +{ + int i, j; + /* Hold the cyclic decomposition of the (multiple) replica + * exchange. */ + gmx_bool bAnyReplicaExchanged = FALSE; + *bThisReplicaExchanged = FALSE; + + for (i = 0; i < re->nrepl; i++) + { + if (re->destinations[i] != re->ind[i]) + { + /* only mark as exchanged if the index has been shuffled */ + bAnyReplicaExchanged = TRUE; + break; + } + } + if (bAnyReplicaExchanged) + { + /* reinitialize the placeholder arrays */ + for (i = 0; i < re->nrepl; i++) + { + for (j = 0; j < re->nrepl; j++) + { + re->cyclic[i][j] = -1; + re->order[i][j] = -1; + } + } + + /* Identify the cyclic decomposition of the permutation (very + * fast if neighbor replica exchange). */ + cyclic_decomposition(re->destinations, re->cyclic, re->incycle, re->nrepl, maxswap); + + /* Now translate the decomposition into a replica exchange + * order at each step. */ + compute_exchange_order(re->cyclic, re->order, re->nrepl, *maxswap); + + /* Did this replica do any exchange at any point? */ + for (j = 0; j < *maxswap; j++) + { + if (replica_id != re->order[replica_id][j]) + { + *bThisReplicaExchanged = TRUE; + break; + } + } + } +} + +gmx_bool replica_exchange(FILE* fplog, + const t_commrec* cr, + const gmx_multisim_t* ms, + struct gmx_repl_ex* re, + t_state* state, + const gmx_enerdata_t* enerd, + t_state* state_local, + int64_t step, + real time) +{ + int j; + int replica_id = 0; + int exchange_partner; + int maxswap = 0; + /* Number of rounds of exchanges needed to deal with any multiple + * exchanges. */ + /* Where each replica ends up after the exchange attempt(s). */ + /* The order in which multiple exchanges will occur. */ + gmx_bool bThisReplicaExchanged = FALSE; + + /* PLUMED */ + if (PLMD::plumedOutside().active()) + PLMD::plumedOutside().cmd("GREX prepare", nullptr); + /* END PLUMED */ + + if (MAIN(cr)) + { + replica_id = re->repl; + test_for_replica_exchange(fplog, ms, re, enerd, det(state_local->box), step, time); + prepare_to_do_exchange(re, replica_id, &maxswap, &bThisReplicaExchanged); + } + /* Do intra-simulation broadcast so all processors belonging to + * each simulation know whether they need to participate in + * collecting the state. Otherwise, they might as well get on with + * the next thing to do. */ + if (haveDDAtomOrdering(*cr)) + { +#if GMX_MPI + MPI_Bcast(&bThisReplicaExchanged, sizeof(gmx_bool), MPI_BYTE, MAINRANK(cr), cr->mpi_comm_mygroup); +#endif + } + + if (bThisReplicaExchanged) + { + /* Exchange the states */ + /* Collect the global state on the main node */ + if (haveDDAtomOrdering(*cr)) + { + dd_collect_state(cr->dd, state_local, state); + } + else + { + copy_state_serial(state_local, state); + } + + if (MAIN(cr)) + { + /* There will be only one swap cycle with standard replica + * exchange, but there may be multiple swap cycles if we + * allow multiple swaps. */ + + for (j = 0; j < maxswap; j++) + { + exchange_partner = re->order[replica_id][j]; + + if (exchange_partner != replica_id) + { + /* Exchange the global states between the main nodes */ + if (debug) + { + fprintf(debug, "Exchanging %d with %d\n", replica_id, exchange_partner); + } + exchange_state(ms, exchange_partner, state); + } + } + /* For temperature-type replica exchange, we need to scale + * the velocities. */ + if (re->type == ReplicaExchangeType::Temperature || re->type == ReplicaExchangeType::TemperatureLambda) + { + scale_velocities( + state->v, + std::sqrt(re->q[ReplicaExchangeType::Temperature][replica_id] + / re->q[ReplicaExchangeType::Temperature][re->destinations[replica_id]])); + } + } + + /* With domain decomposition the global state is distributed later */ + if (!haveDDAtomOrdering(*cr)) + { + /* Copy the global state to the local state data structure */ + copy_state_serial(state, state_local); + } + } + + return bThisReplicaExchanged; +} + +void print_replica_exchange_statistics(FILE* fplog, struct gmx_repl_ex* re) +{ + int i; + + fprintf(fplog, "\nReplica exchange statistics\n"); + + if (re->nex == 0) + { + fprintf(fplog, + "Repl %d attempts, %d odd, %d even\n", + re->nattempt[0] + re->nattempt[1], + re->nattempt[1], + re->nattempt[0]); + + fprintf(fplog, "Repl average probabilities:\n"); + for (i = 1; i < re->nrepl; i++) + { + if (re->nattempt[i % 2] == 0) + { + re->prob[i] = 0; + } + else + { + re->prob[i] = re->prob_sum[i] / re->nattempt[i % 2]; + } + } + print_ind(fplog, "", re->nrepl, re->ind, nullptr); + print_prob(fplog, "", re->nrepl, re->prob); + + fprintf(fplog, "Repl number of exchanges:\n"); + print_ind(fplog, "", re->nrepl, re->ind, nullptr); + print_count(fplog, "", re->nrepl, re->nexchange); + + fprintf(fplog, "Repl average number of exchanges:\n"); + for (i = 1; i < re->nrepl; i++) + { + if (re->nattempt[i % 2] == 0) + { + re->prob[i] = 0; + } + else + { + re->prob[i] = (static_cast(re->nexchange[i])) / re->nattempt[i % 2]; + } + } + print_ind(fplog, "", re->nrepl, re->ind, nullptr); + print_prob(fplog, "", re->nrepl, re->prob); + + fprintf(fplog, "\n"); + } + /* print the transition matrix */ + print_transition_matrix(fplog, re->nrepl, re->nmoves, re->nattempt); +} + +//! \endcond diff --git a/patches/gromacs-2025.0.diff/src/gromacs/mdrun/replicaexchange.cpp.preplumed b/patches/gromacs-2025.0.diff/src/gromacs/mdrun/replicaexchange.cpp.preplumed new file mode 100644 index 0000000000..6fdee432d2 --- /dev/null +++ b/patches/gromacs-2025.0.diff/src/gromacs/mdrun/replicaexchange.cpp.preplumed @@ -0,0 +1,1424 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 1991- 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 replica exchange routines. + * + * \author David van der Spoel + * \author Mark Abraham + * \ingroup module_mdrun + */ +#include "gmxpre.h" + +#include "replicaexchange.h" + +#include "config.h" + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "gromacs/domdec/collect.h" +#include "gromacs/gmxlib/network.h" +#include "gromacs/math/functions.h" +#include "gromacs/math/units.h" +#include "gromacs/math/vec.h" +#include "gromacs/math/vectypes.h" +#include "gromacs/mdrunutility/multisim.h" +#include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/enerdata.h" +#include "gromacs/mdtypes/inputrec.h" +#include "gromacs/mdtypes/md_enums.h" +#include "gromacs/mdtypes/state.h" +#include "gromacs/random/seed.h" +#include "gromacs/random/threefry.h" +#include "gromacs/random/uniformintdistribution.h" +#include "gromacs/random/uniformrealdistribution.h" +#include "gromacs/topology/ifunc.h" +#include "gromacs/utility/arrayref.h" +#include "gromacs/utility/enumerationhelpers.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/gmxassert.h" +#include "gromacs/utility/pleasecite.h" +#include "gromacs/utility/smalloc.h" + +//! Helps cut off probability values. +constexpr int c_probabilityCutoff = 100; + +/* we don't bother evaluating if events are more rare than exp(-100) = 3.7x10^-44 */ + +//! Rank in the multisimulation +#define MSRANK(ms, nodeid) (nodeid) + +//! Enum for replica exchange flavours +enum class ReplicaExchangeType : int +{ + Temperature, + Lambda, + EndSingle, + TemperatureLambda, + Count +}; +/*! \brief Strings describing replica exchange flavours. + * + * end_single_marker merely notes the end of single variable replica + * exchange. All types higher than it are multiple replica exchange + * methods. + * + * Eventually, should add 'pressure', 'temperature and pressure', + * 'lambda_and_pressure', 'temperature_lambda_pressure'?; Let's wait + * until we feel better about the pressure control methods giving + * exact ensembles. Right now, we assume constant pressure */ +static const char* enumValueToString(ReplicaExchangeType enumValue) +{ + constexpr gmx::EnumerationArray replicateExchangeTypeNames = { + "temperature", "lambda", "end_single_marker", "temperature and lambda" + }; + return replicateExchangeTypeNames[enumValue]; +} + +//! Working data for replica exchange. +struct gmx_repl_ex +{ + //! Replica ID + int repl; + //! Total number of replica + int nrepl; + //! Temperature + real temp; + //! Replica exchange type from ReplicaExchangeType enum + ReplicaExchangeType type; + //! Quantity, e.g. temperature or lambda; first index is ere, second index is replica ID + gmx::EnumerationArray q; + //! Use constant pressure and temperature + gmx_bool bNPT; + //! Replica pressures + real* pres; + //! Replica indices + int* ind; + //! Used for keeping track of all the replica swaps + int* allswaps; + //! Replica exchange interval (number of steps) + int nst; + //! Number of exchanges per interval + int nex; + //! Random seed + int seed; + //! Number of even and odd replica change attempts + int nattempt[2]; + //! Sum of probabilities + real* prob_sum; + //! Number of moves between replicas i and j + int** nmoves; + //! i-th element of the array is the number of exchanges between replica i-1 and i + int* nexchange; + + /*! \brief Helper arrays for replica exchange; allocated here + * so they don't have to be allocated each time */ + //! \{ + int* destinations; + int** cyclic; + int** order; + int* tmpswap; + gmx_bool* incycle; + gmx_bool* bEx; + //! \} + + //! Helper arrays to hold the quantities that are exchanged. + //! \{ + real* prob; + real* Epot; + real* beta; + real* Vol; + real** de; + //! \} +}; + +// TODO We should add Doxygen here some time. +//! \cond + +static gmx_bool repl_quantity(const gmx_multisim_t* ms, struct gmx_repl_ex* re, ReplicaExchangeType ere, real q) +{ + real* qall; + gmx_bool bDiff; + int s; + + snew(qall, ms->numSimulations_); + qall[re->repl] = q; + gmx_sum_sim(ms->numSimulations_, qall, ms); + + bDiff = FALSE; + for (s = 1; s < ms->numSimulations_; s++) + { + if (std::fabs(qall[s] - qall[0]) > 1e-5) /* Each quantity normally has 2-4 significant digits */ + { + bDiff = TRUE; + } + } + + if (bDiff) + { + /* Set the replica exchange type and quantities */ + re->type = ere; + + snew(re->q[ere], re->nrepl); + for (s = 0; s < ms->numSimulations_; s++) + { + re->q[ere][s] = qall[s]; + } + } + sfree(qall); + return bDiff; +} + +gmx_repl_ex_t init_replica_exchange(FILE* fplog, + const gmx_multisim_t* ms, + int numAtomsInSystem, + const t_inputrec* ir, + const ReplicaExchangeParameters& replExParams) +{ + real pres; + int i, j; + struct gmx_repl_ex* re; + gmx_bool bTemp; + gmx_bool bLambda = FALSE; + + fprintf(fplog, "\nInitializing Replica Exchange\n"); + + if (!isMultiSim(ms) || ms->numSimulations_ == 1) + { + gmx_fatal(FARGS, + "Nothing to exchange with only one replica, maybe you forgot to set the " + "-multidir option of mdrun?"); + } + if (replExParams.numExchanges < 0) + { + gmx_fatal(FARGS, "Replica exchange number of exchanges needs to be positive"); + } + + if (!EI_DYNAMICS(ir->eI)) + { + gmx_fatal(FARGS, "Replica exchange is only supported by dynamical simulations"); + /* Note that PAR(cr) is defined by cr->nnodes > 1, which is + * distinct from isMultiSim(ms). A multi-simulation only runs + * with real MPI parallelism, but this does not imply PAR(cr) + * is true! + * + * Since we are using a dynamical integrator, the only + * decomposition is DD, so PAR(cr) and haveDDAtomOrdering(*cr) are + * synonymous. The only way for cr->nnodes > 1 to be true is + * if we are using DD. */ + } + + if (!haveConstantEnsembleTemperature(*ir)) + { + gmx_fatal(FARGS, + "Replica exchange is only supported for systems that have a constant ensemble " + "temperature"); + } + + snew(re, 1); + + re->repl = ms->simulationIndex_; + re->nrepl = ms->numSimulations_; + + fprintf(fplog, "Repl There are %d replicas:\n", re->nrepl); + + /* We only check that the number of atoms in the systms match. + * This, of course, do not guarantee that the systems are the same, + * but it does guarantee that we can perform replica exchange. + */ + check_multi_int(fplog, ms, numAtomsInSystem, "the number of atoms", FALSE); + check_multi_int(fplog, ms, static_cast(ir->eI), "the integrator", FALSE); + check_multi_int64(fplog, ms, ir->init_step + ir->nsteps, "init_step+nsteps", FALSE); + const int nst = replExParams.exchangeInterval; + check_multi_int64(fplog, + ms, + gmx::divideRoundUp(ir->init_step, nst), + "first exchange step: init_step/-replex", + FALSE); + check_multi_int(fplog, ms, static_cast(ir->etc), "the temperature coupling", FALSE); + check_multi_int(fplog, ms, ir->opts.ngtc, "the number of temperature coupling groups", FALSE); + check_multi_int( + fplog, ms, static_cast(ir->pressureCouplingOptions.epc), "the pressure coupling", FALSE); + check_multi_int(fplog, ms, static_cast(ir->efep), "free energy", FALSE); + check_multi_int(fplog, ms, ir->fepvals->n_lambda, "number of lambda states", FALSE); + + re->temp = constantEnsembleTemperature(*ir); + for (i = 1; (i < ir->opts.ngtc); i++) + { + if (ir->opts.ref_t[i] != re->temp) + { + fprintf(fplog, + "\nWARNING: The temperatures of the different temperature coupling groups are " + "not identical\n\n"); + fprintf(stderr, + "\nWARNING: The temperatures of the different temperature coupling groups are " + "not identical\n\n"); + } + } + + re->type = ReplicaExchangeType::Count; + bTemp = repl_quantity(ms, re, ReplicaExchangeType::Temperature, re->temp); + if (ir->efep != FreeEnergyPerturbationType::No) + { + bLambda = repl_quantity( + ms, re, ReplicaExchangeType::Lambda, static_cast(ir->fepvals->init_fep_state)); + } + if (re->type == ReplicaExchangeType::Count) /* nothing was assigned */ + { + gmx_fatal(FARGS, + "The properties of the %d systems are all the same, there is nothing to exchange", + re->nrepl); + } + if (bLambda && bTemp) + { + re->type = ReplicaExchangeType::TemperatureLambda; + } + + if (bTemp) + { + please_cite(fplog, "Sugita1999a"); + if (ir->pressureCouplingOptions.epc != PressureCoupling::No) + { + re->bNPT = TRUE; + fprintf(fplog, "Repl Using Constant Pressure REMD.\n"); + please_cite(fplog, "Okabe2001a"); + } + if (ir->etc == TemperatureCoupling::Berendsen) + { + gmx_fatal(FARGS, + "REMD with the %s thermostat does not produce correct potential energy " + "distributions, consider using the %s thermostat instead", + enumValueToString(ir->etc), + enumValueToString(TemperatureCoupling::VRescale)); + } + } + if (bLambda) + { + if (ir->fepvals->delta_lambda != 0) /* check this? */ + { + gmx_fatal(FARGS, "delta_lambda is not zero"); + } + } + if (re->bNPT) + { + snew(re->pres, re->nrepl); + if (ir->pressureCouplingOptions.epct == PressureCouplingType::SurfaceTension) + { + pres = ir->pressureCouplingOptions.ref_p[ZZ][ZZ]; + } + else + { + pres = 0; + j = 0; + for (i = 0; i < DIM; i++) + { + if (ir->pressureCouplingOptions.compress[i][i] != 0) + { + pres += ir->pressureCouplingOptions.ref_p[i][i]; + j++; + } + } + pres /= j; + } + re->pres[re->repl] = pres; + gmx_sum_sim(re->nrepl, re->pres, ms); + } + + /* Make an index for increasing replica order */ + /* only makes sense if one or the other is varying, not both! + if both are varying, we trust the order the person gave. */ + snew(re->ind, re->nrepl); + for (i = 0; i < re->nrepl; i++) + { + re->ind[i] = i; + } + + if (re->type < ReplicaExchangeType::EndSingle) + { + + for (i = 0; i < re->nrepl; i++) + { + for (j = i + 1; j < re->nrepl; j++) + { + if (re->q[re->type][re->ind[j]] < re->q[re->type][re->ind[i]]) + { + /* Unordered replicas are supposed to work, but there + * is still an issues somewhere. + * Note that at this point still re->ind[i]=i. + */ + gmx_fatal(FARGS, + "Replicas with indices %d < %d have %ss %g > %g, please order your " + "replicas on increasing %s", + i, + j, + enumValueToString(re->type), + re->q[re->type][i], + re->q[re->type][j], + enumValueToString(re->type)); + } + else if (re->q[re->type][re->ind[j]] == re->q[re->type][re->ind[i]]) + { + gmx_fatal(FARGS, "Two replicas have identical %ss", enumValueToString(re->type)); + } + } + } + } + + /* keep track of all the swaps, starting with the initial placement. */ + snew(re->allswaps, re->nrepl); + for (i = 0; i < re->nrepl; i++) + { + re->allswaps[i] = re->ind[i]; + } + + switch (re->type) + { + case ReplicaExchangeType::Temperature: + fprintf(fplog, "\nReplica exchange in temperature\n"); + for (i = 0; i < re->nrepl; i++) + { + fprintf(fplog, " %5.1f", re->q[re->type][re->ind[i]]); + } + fprintf(fplog, "\n"); + break; + case ReplicaExchangeType::Lambda: + fprintf(fplog, "\nReplica exchange in lambda\n"); + for (i = 0; i < re->nrepl; i++) + { + fprintf(fplog, " %3d", static_cast(re->q[re->type][re->ind[i]])); + } + fprintf(fplog, "\n"); + break; + case ReplicaExchangeType::TemperatureLambda: + fprintf(fplog, "\nReplica exchange in temperature and lambda state\n"); + for (i = 0; i < re->nrepl; i++) + { + fprintf(fplog, " %5.1f", re->q[ReplicaExchangeType::Temperature][re->ind[i]]); + } + fprintf(fplog, "\n"); + for (i = 0; i < re->nrepl; i++) + { + fprintf(fplog, " %5d", static_cast(re->q[ReplicaExchangeType::Lambda][re->ind[i]])); + } + fprintf(fplog, "\n"); + break; + default: gmx_incons("Unknown replica exchange quantity"); + } + if (re->bNPT) + { + fprintf(fplog, "\nRepl p"); + for (i = 0; i < re->nrepl; i++) + { + fprintf(fplog, " %5.2f", re->pres[re->ind[i]]); + } + + for (i = 0; i < re->nrepl; i++) + { + if ((i > 0) && (re->pres[re->ind[i]] < re->pres[re->ind[i - 1]])) + { + fprintf(fplog, + "\nWARNING: The reference pressures decrease with increasing " + "temperatures\n\n"); + fprintf(stderr, + "\nWARNING: The reference pressures decrease with increasing " + "temperatures\n\n"); + } + } + } + re->nst = nst; + if (replExParams.randomSeed == -1) + { + if (isMainSim(ms)) + { + re->seed = static_cast(gmx::makeRandomSeed()); + } + else + { + re->seed = 0; + } + gmx_sumi_sim(1, &(re->seed), ms); + } + else + { + re->seed = replExParams.randomSeed; + } + fprintf(fplog, "\nReplica exchange interval: %d\n", re->nst); + fprintf(fplog, "\nReplica random seed: %d\n", re->seed); + + re->nattempt[0] = 0; + re->nattempt[1] = 0; + + snew(re->prob_sum, re->nrepl); + snew(re->nexchange, re->nrepl); + snew(re->nmoves, re->nrepl); + for (i = 0; i < re->nrepl; i++) + { + snew(re->nmoves[i], re->nrepl); + } + fprintf(fplog, "Replica exchange information below: ex and x = exchange, pr = probability\n"); + + /* generate space for the helper functions so we don't have to snew each time */ + + snew(re->destinations, re->nrepl); + snew(re->incycle, re->nrepl); + snew(re->tmpswap, re->nrepl); + snew(re->cyclic, re->nrepl); + snew(re->order, re->nrepl); + for (i = 0; i < re->nrepl; i++) + { + snew(re->cyclic[i], re->nrepl + 1); + snew(re->order[i], re->nrepl); + } + /* allocate space for the functions storing the data for the replicas */ + /* not all of these arrays needed in all cases, but they don't take + up much space, since the max size is nrepl**2 */ + snew(re->prob, re->nrepl); + snew(re->bEx, re->nrepl); + snew(re->beta, re->nrepl); + snew(re->Vol, re->nrepl); + snew(re->Epot, re->nrepl); + snew(re->de, re->nrepl); + for (i = 0; i < re->nrepl; i++) + { + snew(re->de[i], re->nrepl); + } + re->nex = replExParams.numExchanges; + return re; +} + +static void exchange_reals(const gmx_multisim_t gmx_unused* ms, int gmx_unused b, real* v, int n) +{ + real* buf; + int i; + + if (v) + { + snew(buf, n); +#if GMX_MPI + /* + MPI_Sendrecv(v, n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0, + buf,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0, + ms->mainRanksComm_,MPI_STATUS_IGNORE); + */ + { + MPI_Request mpi_req; + + MPI_Isend(v, n * sizeof(real), MPI_BYTE, MSRANK(ms, b), 0, ms->mainRanksComm_, &mpi_req); + MPI_Recv(buf, n * sizeof(real), MPI_BYTE, MSRANK(ms, b), 0, ms->mainRanksComm_, MPI_STATUS_IGNORE); + MPI_Wait(&mpi_req, MPI_STATUS_IGNORE); + } +#endif + for (i = 0; i < n; i++) + { + v[i] = buf[i]; + } + sfree(buf); + } +} + + +static void exchange_doubles(const gmx_multisim_t gmx_unused* ms, int gmx_unused b, double* v, int n) +{ + double* buf; + int i; + + if (v) + { + snew(buf, n); +#if GMX_MPI + /* + MPI_Sendrecv(v, n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0, + buf,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0, + ms->mainRanksComm_,MPI_STATUS_IGNORE); + */ + { + MPI_Request mpi_req; + + MPI_Isend(v, n * sizeof(double), MPI_BYTE, MSRANK(ms, b), 0, ms->mainRanksComm_, &mpi_req); + MPI_Recv(buf, n * sizeof(double), MPI_BYTE, MSRANK(ms, b), 0, ms->mainRanksComm_, MPI_STATUS_IGNORE); + MPI_Wait(&mpi_req, MPI_STATUS_IGNORE); + } +#endif + for (i = 0; i < n; i++) + { + v[i] = buf[i]; + } + sfree(buf); + } +} + +static void exchange_rvecs(const gmx_multisim_t gmx_unused* ms, int gmx_unused b, rvec* v, int n) +{ + rvec* buf; + int i; + + if (v) + { + snew(buf, n); +#if GMX_MPI + /* + MPI_Sendrecv(v[0], n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0, + buf[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0, + ms->mainRanksComm_,MPI_STATUS_IGNORE); + */ + { + MPI_Request mpi_req; + + MPI_Isend(v[0], n * sizeof(rvec), MPI_BYTE, MSRANK(ms, b), 0, ms->mainRanksComm_, &mpi_req); + MPI_Recv(buf[0], n * sizeof(rvec), MPI_BYTE, MSRANK(ms, b), 0, ms->mainRanksComm_, MPI_STATUS_IGNORE); + MPI_Wait(&mpi_req, MPI_STATUS_IGNORE); + } +#endif + for (i = 0; i < n; i++) + { + copy_rvec(buf[i], v[i]); + } + sfree(buf); + } +} + +static void exchange_state(const gmx_multisim_t* ms, int b, t_state* state) +{ + /* When t_state changes, this code should be updated. */ + int ngtc, nnhpres; + ngtc = state->ngtc * state->nhchainlength; + nnhpres = state->nnhpres * state->nhchainlength; + exchange_rvecs(ms, b, state->box, DIM); + exchange_rvecs(ms, b, state->box_rel, DIM); + exchange_rvecs(ms, b, state->boxv, DIM); + exchange_reals(ms, b, &(state->veta), 1); + exchange_reals(ms, b, &(state->vol0), 1); + exchange_rvecs(ms, b, state->svir_prev, DIM); + exchange_rvecs(ms, b, state->fvir_prev, DIM); + exchange_rvecs(ms, b, state->pres_prev, DIM); + exchange_doubles(ms, b, state->nosehoover_xi.data(), ngtc); + exchange_doubles(ms, b, state->nosehoover_vxi.data(), ngtc); + exchange_doubles(ms, b, state->nhpres_xi.data(), nnhpres); + exchange_doubles(ms, b, state->nhpres_vxi.data(), nnhpres); + exchange_doubles(ms, b, state->therm_integral.data(), state->ngtc); + exchange_doubles(ms, b, &state->baros_integral, 1); + exchange_rvecs(ms, b, state->x.rvec_array(), state->numAtoms()); + exchange_rvecs(ms, b, state->v.rvec_array(), state->numAtoms()); +} + +static void copy_state_serial(const t_state* src, t_state* dest) +{ + if (dest != src) + { + /* Currently the local state is always a pointer to the global + * in serial, so we should never end up here. + * TODO: Implement a (trivial) t_state copy once converted to C++. + */ + GMX_RELEASE_ASSERT(false, "State copying is currently not implemented in replica exchange"); + } +} + +static void scale_velocities(gmx::ArrayRef velocities, real fac) +{ + for (auto& v : velocities) + { + v *= fac; + } +} + +static void print_transition_matrix(FILE* fplog, int n, int** nmoves, const int* nattempt) +{ + int i, j, ntot; + float Tprint; + + ntot = nattempt[0] + nattempt[1]; + fprintf(fplog, "\n"); + fprintf(fplog, "Repl"); + for (i = 0; i < n; i++) + { + fprintf(fplog, " "); /* put the title closer to the center */ + } + fprintf(fplog, "Empirical Transition Matrix\n"); + + fprintf(fplog, "Repl"); + for (i = 0; i < n; i++) + { + fprintf(fplog, "%8d", (i + 1)); + } + fprintf(fplog, "\n"); + + for (i = 0; i < n; i++) + { + fprintf(fplog, "Repl"); + for (j = 0; j < n; j++) + { + Tprint = 0.0; + if (nmoves[i][j] > 0) + { + Tprint = nmoves[i][j] / (2.0 * ntot); + } + fprintf(fplog, "%8.4f", Tprint); + } + fprintf(fplog, "%3d\n", i); + } +} + +static void print_ind(FILE* fplog, const char* leg, int n, int* ind, const gmx_bool* bEx) +{ + int i; + + fprintf(fplog, "Repl %2s %2d", leg, ind[0]); + for (i = 1; i < n; i++) + { + fprintf(fplog, " %c %2d", (bEx != nullptr && bEx[i]) ? 'x' : ' ', ind[i]); + } + fprintf(fplog, "\n"); +} + +static void print_allswitchind(FILE* fplog, int n, int* pind, int* allswaps, int* tmpswap) +{ + int i; + + for (i = 0; i < n; i++) + { + tmpswap[i] = allswaps[i]; + } + for (i = 0; i < n; i++) + { + allswaps[i] = tmpswap[pind[i]]; + } + + fprintf(fplog, "\nAccepted Exchanges: "); + for (i = 0; i < n; i++) + { + fprintf(fplog, "%d ", pind[i]); + } + fprintf(fplog, "\n"); + + /* the "Order After Exchange" is the state label corresponding to the configuration that + started in state listed in order, i.e. + + 3 0 1 2 + + means that the: + configuration starting in simulation 3 is now in simulation 0, + configuration starting in simulation 0 is now in simulation 1, + configuration starting in simulation 1 is now in simulation 2, + configuration starting in simulation 2 is now in simulation 3 + */ + fprintf(fplog, "Order After Exchange: "); + for (i = 0; i < n; i++) + { + fprintf(fplog, "%d ", allswaps[i]); + } + fprintf(fplog, "\n\n"); +} + +static void print_prob(FILE* fplog, const char* leg, int n, real* prob) +{ + int i; + char buf[8]; + + fprintf(fplog, "Repl %2s ", leg); + for (i = 1; i < n; i++) + { + if (prob[i] >= 0) + { + sprintf(buf, "%4.2f", prob[i]); + fprintf(fplog, " %3s", buf[0] == '1' ? "1.0" : buf + 1); + } + else + { + fprintf(fplog, " "); + } + } + fprintf(fplog, "\n"); +} + +static void print_count(FILE* fplog, const char* leg, int n, int* count) +{ + int i; + + fprintf(fplog, "Repl %2s ", leg); + for (i = 1; i < n; i++) + { + fprintf(fplog, " %4d", count[i]); + } + fprintf(fplog, "\n"); +} + +static real calc_delta(FILE* fplog, gmx_bool bPrint, struct gmx_repl_ex* re, int a, int b, int ap, int bp) +{ + + real ediff, dpV, delta = 0; + real* Epot = re->Epot; + real* Vol = re->Vol; + real** de = re->de; + real* beta = re->beta; + + /* Two cases; we are permuted and not. In all cases, setting ap = a and bp = b will reduce + to the non permuted case */ + + switch (re->type) + { + case ReplicaExchangeType::Temperature: + /* + * Okabe et. al. Chem. Phys. Lett. 335 (2001) 435-439 + */ + ediff = Epot[b] - Epot[a]; + delta = -(beta[bp] - beta[ap]) * ediff; + break; + case ReplicaExchangeType::Lambda: + /* two cases: when we are permuted, and not. */ + /* non-permuted: + ediff = E_new - E_old + = [H_b(x_a) + H_a(x_b)] - [H_b(x_b) + H_a(x_a)] + = [H_b(x_a) - H_a(x_a)] + [H_a(x_b) - H_b(x_b)] + = de[b][a] + de[a][b] */ + + /* permuted: + ediff = E_new - E_old + = [H_bp(x_a) + H_ap(x_b)] - [H_bp(x_b) + H_ap(x_a)] + = [H_bp(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_bp(x_b)] + = [H_bp(x_a) - H_a(x_a) + H_a(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_b(x_b) + H_b(x_b) - H_bp(x_b)] + = [H_bp(x_a) - H_a(x_a)] - [H_ap(x_a) - H_a(x_a)] + [H_ap(x_b) - H_b(x_b)] - H_bp(x_b) - H_b(x_b)] + = (de[bp][a] - de[ap][a]) + (de[ap][b] - de[bp][b]) */ + /* but, in the current code implementation, we flip configurations, not indices . . . + So let's examine that. + = [H_b(x_ap) - H_a(x_a)] - [H_a(x_ap) - H_a(x_a)] + [H_a(x_bp) - H_b(x_b)] - H_b(x_bp) - H_b(x_b)] + = [H_b(x_ap) - H_a(x_ap)] + [H_a(x_bp) - H_b(x_pb)] + = (de[b][ap] - de[a][ap]) + (de[a][bp] - de[b][bp] + So, if we exchange b<=> bp and a<=> ap, we return to the same result. + So the simple solution is to flip the + position of perturbed and original indices in the tests. + */ + + ediff = (de[bp][a] - de[ap][a]) + (de[ap][b] - de[bp][b]); + delta = ediff * beta[a]; /* assume all same temperature in this case */ + break; + case ReplicaExchangeType::TemperatureLambda: + /* not permuted: */ + /* delta = reduced E_new - reduced E_old + = [beta_b H_b(x_a) + beta_a H_a(x_b)] - [beta_b H_b(x_b) + beta_a H_a(x_a)] + = [beta_b H_b(x_a) - beta_a H_a(x_a)] + [beta_a H_a(x_b) - beta_b H_b(x_b)] + = [beta_b dH_b(x_a) + beta_b H_a(x_a) - beta_a H_a(x_a)] + + [beta_a dH_a(x_b) + beta_a H_b(x_b) - beta_b H_b(x_b)] + = [beta_b dH_b(x_a) + [beta_a dH_a(x_b) + + beta_b (H_a(x_a) - H_b(x_b)]) - beta_a (H_a(x_a) - H_b(x_b)) + = beta_b dH_b(x_a) + beta_a dH_a(x_b) - (beta_b - beta_a)(H_b(x_b) - H_a(x_a) */ + /* delta = beta[b]*de[b][a] + beta[a]*de[a][b] - (beta[b] - beta[a])*(Epot[b] - Epot[a]; */ + /* permuted (big breath!) */ + /* delta = reduced E_new - reduced E_old + = [beta_bp H_bp(x_a) + beta_ap H_ap(x_b)] - [beta_bp H_bp(x_b) + beta_ap H_ap(x_a)] + = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)] + = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)] + - beta_pb H_a(x_a) + beta_ap H_a(x_a) + beta_pb H_a(x_a) - beta_ap H_a(x_a) + - beta_ap H_b(x_b) + beta_bp H_b(x_b) + beta_ap H_b(x_b) - beta_bp H_b(x_b) + = [(beta_bp H_bp(x_a) - beta_bp H_a(x_a)) - (beta_ap H_ap(x_a) - beta_ap H_a(x_a))] + + [(beta_ap H_ap(x_b) - beta_ap H_b(x_b)) - (beta_bp H_bp(x_b) - beta_bp H_b(x_b))] + + beta_pb H_a(x_a) - beta_ap H_a(x_a) + beta_ap H_b(x_b) - beta_bp H_b(x_b) + = [beta_bp (H_bp(x_a) - H_a(x_a)) - beta_ap (H_ap(x_a) - H_a(x_a))] + + [beta_ap (H_ap(x_b) - H_b(x_b)) - beta_bp (H_bp(x_b) - H_b(x_b))] + + beta_pb (H_a(x_a) - H_b(x_b)) - beta_ap (H_a(x_a) - H_b(x_b)) + = ([beta_bp de[bp][a] - beta_ap de[ap][a]) + beta_ap de[ap][b] - beta_bp de[bp][b]) + + (beta_pb-beta_ap)(H_a(x_a) - H_b(x_b)) */ + delta = beta[bp] * (de[bp][a] - de[bp][b]) + beta[ap] * (de[ap][b] - de[ap][a]) + - (beta[bp] - beta[ap]) * (Epot[b] - Epot[a]); + break; + default: gmx_incons("Unknown replica exchange quantity"); + } + if (bPrint) + { + fprintf(fplog, "Repl %d <-> %d dE_term = %10.3e (kT)\n", a, b, delta); + } + if (re->bNPT) + { + /* revist the calculation for 5.0. Might be some improvements. */ + dpV = (beta[ap] * re->pres[ap] - beta[bp] * re->pres[bp]) * (Vol[b] - Vol[a]) / gmx::c_presfac; + if (bPrint) + { + fprintf(fplog, " dpV = %10.3e d = %10.3e\n", dpV, delta + dpV); + } + delta += dpV; + } + return delta; +} + +static void test_for_replica_exchange(FILE* fplog, + const gmx_multisim_t* ms, + struct gmx_repl_ex* re, + const gmx_enerdata_t* enerd, + real vol, + int64_t step, + real time) +{ + int m, i, j, a, b, ap, bp, i0, i1, tmp; + real delta = 0; + gmx_bool bPrint, bMultiEx; + gmx_bool* bEx = re->bEx; + real* prob = re->prob; + int* pind = re->destinations; /* permuted index */ + gmx_bool bEpot = FALSE; + gmx_bool bDLambda = FALSE; + gmx_bool bVol = FALSE; + gmx::ThreeFry2x64<64> rng(re->seed, gmx::RandomDomain::ReplicaExchange); + gmx::UniformRealDistribution uniformRealDist; + gmx::UniformIntDistribution uniformNreplDist(0, re->nrepl - 1); + + bMultiEx = (re->nex > 1); /* multiple exchanges at each state */ + fprintf(fplog, "Replica exchange at step %" PRId64 " time %.5f\n", step, time); + + if (re->bNPT) + { + for (i = 0; i < re->nrepl; i++) + { + re->Vol[i] = 0; + } + bVol = TRUE; + re->Vol[re->repl] = vol; + } + if ((re->type == ReplicaExchangeType::Temperature || re->type == ReplicaExchangeType::TemperatureLambda)) + { + for (i = 0; i < re->nrepl; i++) + { + re->Epot[i] = 0; + } + bEpot = TRUE; + re->Epot[re->repl] = enerd->term[F_EPOT]; + /* temperatures of different states*/ + for (i = 0; i < re->nrepl; i++) + { + re->beta[i] = 1.0 / (re->q[ReplicaExchangeType::Temperature][i] * gmx::c_boltz); + } + } + else + { + for (i = 0; i < re->nrepl; i++) + { + re->beta[i] = 1.0 / (re->temp * gmx::c_boltz); /* we have a single temperature */ + } + } + if (re->type == ReplicaExchangeType::Lambda || re->type == ReplicaExchangeType::TemperatureLambda) + { + bDLambda = TRUE; + /* lambda differences. */ + /* de[i][j] is the energy of the jth simulation in the ith Hamiltonian + minus the energy of the jth simulation in the jth Hamiltonian */ + for (i = 0; i < re->nrepl; i++) + { + for (j = 0; j < re->nrepl; j++) + { + re->de[i][j] = 0; + } + } + for (i = 0; i < re->nrepl; i++) + { + re->de[i][re->repl] = + enerd->foreignLambdaTerms.deltaH(re->q[ReplicaExchangeType::Lambda][i]); + } + } + + /* now actually do the communication */ + if (bVol) + { + gmx_sum_sim(re->nrepl, re->Vol, ms); + } + if (bEpot) + { + gmx_sum_sim(re->nrepl, re->Epot, ms); + } + if (bDLambda) + { + for (i = 0; i < re->nrepl; i++) + { + gmx_sum_sim(re->nrepl, re->de[i], ms); + } + } + + /* make a duplicate set of indices for shuffling */ + for (i = 0; i < re->nrepl; i++) + { + pind[i] = re->ind[i]; + } + + rng.restart(step, 0); + + if (bMultiEx) + { + /* multiple random switch exchange */ + int nself = 0; + + + for (i = 0; i < re->nex + nself; i++) + { + // For now this is superfluous, but just in case we ever add more + // calls in different branches it is safer to always reset the distribution. + uniformNreplDist.reset(); + + /* randomly select a pair */ + /* in theory, could reduce this by identifying only which switches had a nonneglibible + probability of occurring (log p > -100) and only operate on those switches */ + /* find out which state it is from, and what label that state currently has. Likely + more work that useful. */ + i0 = uniformNreplDist(rng); + i1 = uniformNreplDist(rng); + if (i0 == i1) + { + nself++; + continue; /* self-exchange, back up and do it again */ + } + + a = re->ind[i0]; /* what are the indices of these states? */ + b = re->ind[i1]; + ap = pind[i0]; + bp = pind[i1]; + + bPrint = FALSE; /* too noisy */ + /* calculate the energy difference */ + /* if the code changes to flip the STATES, rather than the configurations, + use the commented version of the code */ + /* delta = calc_delta(fplog,bPrint,re,a,b,ap,bp); */ + delta = calc_delta(fplog, bPrint, re, ap, bp, a, b); + + /* we actually only use the first space in the prob and bEx array, + since there are actually many switches between pairs. */ + + if (delta <= 0) + { + /* accepted */ + prob[0] = 1; + bEx[0] = TRUE; + } + else + { + if (delta > c_probabilityCutoff) + { + prob[0] = 0; + } + else + { + prob[0] = std::exp(-delta); + } + // roll a number to determine if accepted. For now it is superfluous to + // reset, but just in case we ever add more calls in different branches + // it is safer to always reset the distribution. + uniformRealDist.reset(); + bEx[0] = uniformRealDist(rng) < prob[0]; + } + re->prob_sum[0] += prob[0]; + + if (bEx[0]) + { + /* swap the states */ + tmp = pind[i0]; + pind[i0] = pind[i1]; + pind[i1] = tmp; + } + } + re->nattempt[0]++; /* keep track of total permutation trials here */ + print_allswitchind(fplog, re->nrepl, pind, re->allswaps, re->tmpswap); + } + else + { + /* standard nearest neighbor replica exchange */ + + m = (step / re->nst) % 2; + for (i = 1; i < re->nrepl; i++) + { + a = re->ind[i - 1]; + b = re->ind[i]; + + bPrint = (re->repl == a || re->repl == b); + if (i % 2 == m) + { + delta = calc_delta(fplog, bPrint, re, a, b, a, b); + if (delta <= 0) + { + /* accepted */ + prob[i] = 1; + bEx[i] = TRUE; + } + else + { + if (delta > c_probabilityCutoff) + { + prob[i] = 0; + } + else + { + prob[i] = std::exp(-delta); + } + // roll a number to determine if accepted. For now it is superfluous to + // reset, but just in case we ever add more calls in different branches + // it is safer to always reset the distribution. + uniformRealDist.reset(); + bEx[i] = uniformRealDist(rng) < prob[i]; + } + re->prob_sum[i] += prob[i]; + + if (bEx[i]) + { + /* swap these two */ + tmp = pind[i - 1]; + pind[i - 1] = pind[i]; + pind[i] = tmp; + re->nexchange[i]++; /* statistics for back compatibility */ + } + } + else + { + prob[i] = -1; + bEx[i] = FALSE; + } + } + /* print some statistics */ + print_ind(fplog, "ex", re->nrepl, re->ind, bEx); + print_prob(fplog, "pr", re->nrepl, prob); + fprintf(fplog, "\n"); + re->nattempt[m]++; + } + + /* record which moves were made and accepted */ + for (i = 0; i < re->nrepl; i++) + { + re->nmoves[re->ind[i]][pind[i]] += 1; + re->nmoves[pind[i]][re->ind[i]] += 1; + } + fflush(fplog); /* make sure we can see what the last exchange was */ +} + +static void cyclic_decomposition(const int* destinations, int** cyclic, gmx_bool* incycle, const int nrepl, int* nswap) +{ + + int i, j, c, p; + int maxlen = 1; + for (i = 0; i < nrepl; i++) + { + incycle[i] = FALSE; + } + for (i = 0; i < nrepl; i++) /* one cycle for each replica */ + { + if (incycle[i]) + { + cyclic[i][0] = -1; + continue; + } + cyclic[i][0] = i; + incycle[i] = TRUE; + c = 1; + p = i; + for (j = 0; j < nrepl; j++) /* potentially all cycles are part, but we will break first */ + { + p = destinations[p]; /* start permuting */ + if (p == i) + { + cyclic[i][c] = -1; + if (c > maxlen) + { + maxlen = c; + } + break; /* we've reached the original element, the cycle is complete, and we marked the end. */ + } + else + { + cyclic[i][c] = p; /* each permutation gives a new member of the cycle */ + incycle[p] = TRUE; + c++; + } + } + } + *nswap = maxlen - 1; + + if (debug) + { + for (i = 0; i < nrepl; i++) + { + fprintf(debug, "Cycle %d:", i); + for (j = 0; j < nrepl; j++) + { + if (cyclic[i][j] < 0) + { + break; + } + fprintf(debug, "%2d", cyclic[i][j]); + } + fprintf(debug, "\n"); + } + fflush(debug); + } +} + +static void compute_exchange_order(int** cyclic, int** order, const int nrepl, const int maxswap) +{ + int i, j; + + for (j = 0; j < maxswap; j++) + { + for (i = 0; i < nrepl; i++) + { + if (cyclic[i][j + 1] >= 0) + { + order[cyclic[i][j + 1]][j] = cyclic[i][j]; + order[cyclic[i][j]][j] = cyclic[i][j + 1]; + } + } + for (i = 0; i < nrepl; i++) + { + if (order[i][j] < 0) + { + order[i][j] = i; /* if it's not exchanging, it should stay this round*/ + } + } + } + + if (debug) + { + fprintf(debug, "Replica Exchange Order\n"); + for (i = 0; i < nrepl; i++) + { + fprintf(debug, "Replica %d:", i); + for (j = 0; j < maxswap; j++) + { + if (order[i][j] < 0) + { + break; + } + fprintf(debug, "%2d", order[i][j]); + } + fprintf(debug, "\n"); + } + fflush(debug); + } +} + +static void prepare_to_do_exchange(struct gmx_repl_ex* re, const int replica_id, int* maxswap, gmx_bool* bThisReplicaExchanged) +{ + int i, j; + /* Hold the cyclic decomposition of the (multiple) replica + * exchange. */ + gmx_bool bAnyReplicaExchanged = FALSE; + *bThisReplicaExchanged = FALSE; + + for (i = 0; i < re->nrepl; i++) + { + if (re->destinations[i] != re->ind[i]) + { + /* only mark as exchanged if the index has been shuffled */ + bAnyReplicaExchanged = TRUE; + break; + } + } + if (bAnyReplicaExchanged) + { + /* reinitialize the placeholder arrays */ + for (i = 0; i < re->nrepl; i++) + { + for (j = 0; j < re->nrepl; j++) + { + re->cyclic[i][j] = -1; + re->order[i][j] = -1; + } + } + + /* Identify the cyclic decomposition of the permutation (very + * fast if neighbor replica exchange). */ + cyclic_decomposition(re->destinations, re->cyclic, re->incycle, re->nrepl, maxswap); + + /* Now translate the decomposition into a replica exchange + * order at each step. */ + compute_exchange_order(re->cyclic, re->order, re->nrepl, *maxswap); + + /* Did this replica do any exchange at any point? */ + for (j = 0; j < *maxswap; j++) + { + if (replica_id != re->order[replica_id][j]) + { + *bThisReplicaExchanged = TRUE; + break; + } + } + } +} + +gmx_bool replica_exchange(FILE* fplog, + const t_commrec* cr, + const gmx_multisim_t* ms, + struct gmx_repl_ex* re, + t_state* state, + const gmx_enerdata_t* enerd, + t_state* state_local, + int64_t step, + real time) +{ + int j; + int replica_id = 0; + int exchange_partner; + int maxswap = 0; + /* Number of rounds of exchanges needed to deal with any multiple + * exchanges. */ + /* Where each replica ends up after the exchange attempt(s). */ + /* The order in which multiple exchanges will occur. */ + gmx_bool bThisReplicaExchanged = FALSE; + + if (MAIN(cr)) + { + replica_id = re->repl; + test_for_replica_exchange(fplog, ms, re, enerd, det(state_local->box), step, time); + prepare_to_do_exchange(re, replica_id, &maxswap, &bThisReplicaExchanged); + } + /* Do intra-simulation broadcast so all processors belonging to + * each simulation know whether they need to participate in + * collecting the state. Otherwise, they might as well get on with + * the next thing to do. */ + if (haveDDAtomOrdering(*cr)) + { +#if GMX_MPI + MPI_Bcast(&bThisReplicaExchanged, sizeof(gmx_bool), MPI_BYTE, MAINRANK(cr), cr->mpi_comm_mygroup); +#endif + } + + if (bThisReplicaExchanged) + { + /* Exchange the states */ + /* Collect the global state on the main node */ + if (haveDDAtomOrdering(*cr)) + { + dd_collect_state(cr->dd, state_local, state); + } + else + { + copy_state_serial(state_local, state); + } + + if (MAIN(cr)) + { + /* There will be only one swap cycle with standard replica + * exchange, but there may be multiple swap cycles if we + * allow multiple swaps. */ + + for (j = 0; j < maxswap; j++) + { + exchange_partner = re->order[replica_id][j]; + + if (exchange_partner != replica_id) + { + /* Exchange the global states between the main nodes */ + if (debug) + { + fprintf(debug, "Exchanging %d with %d\n", replica_id, exchange_partner); + } + exchange_state(ms, exchange_partner, state); + } + } + /* For temperature-type replica exchange, we need to scale + * the velocities. */ + if (re->type == ReplicaExchangeType::Temperature || re->type == ReplicaExchangeType::TemperatureLambda) + { + scale_velocities( + state->v, + std::sqrt(re->q[ReplicaExchangeType::Temperature][replica_id] + / re->q[ReplicaExchangeType::Temperature][re->destinations[replica_id]])); + } + } + + /* With domain decomposition the global state is distributed later */ + if (!haveDDAtomOrdering(*cr)) + { + /* Copy the global state to the local state data structure */ + copy_state_serial(state, state_local); + } + } + + return bThisReplicaExchanged; +} + +void print_replica_exchange_statistics(FILE* fplog, struct gmx_repl_ex* re) +{ + int i; + + fprintf(fplog, "\nReplica exchange statistics\n"); + + if (re->nex == 0) + { + fprintf(fplog, + "Repl %d attempts, %d odd, %d even\n", + re->nattempt[0] + re->nattempt[1], + re->nattempt[1], + re->nattempt[0]); + + fprintf(fplog, "Repl average probabilities:\n"); + for (i = 1; i < re->nrepl; i++) + { + if (re->nattempt[i % 2] == 0) + { + re->prob[i] = 0; + } + else + { + re->prob[i] = re->prob_sum[i] / re->nattempt[i % 2]; + } + } + print_ind(fplog, "", re->nrepl, re->ind, nullptr); + print_prob(fplog, "", re->nrepl, re->prob); + + fprintf(fplog, "Repl number of exchanges:\n"); + print_ind(fplog, "", re->nrepl, re->ind, nullptr); + print_count(fplog, "", re->nrepl, re->nexchange); + + fprintf(fplog, "Repl average number of exchanges:\n"); + for (i = 1; i < re->nrepl; i++) + { + if (re->nattempt[i % 2] == 0) + { + re->prob[i] = 0; + } + else + { + re->prob[i] = (static_cast(re->nexchange[i])) / re->nattempt[i % 2]; + } + } + print_ind(fplog, "", re->nrepl, re->ind, nullptr); + print_prob(fplog, "", re->nrepl, re->prob); + + fprintf(fplog, "\n"); + } + /* print the transition matrix */ + print_transition_matrix(fplog, re->nrepl, re->nmoves, re->nattempt); +} + +//! \endcond diff --git a/patches/patch.sh b/patches/patch.sh index 0a7e4bfa50..aae5a84432 100755 --- a/patches/patch.sh +++ b/patches/patch.sh @@ -288,6 +288,7 @@ case "$action" in if [[ -z $PLUMED_ONLY_PATCH ]]; then #PLUMED_ONLY_PATCH is special for programs that have an embedded plumed integration + #and do not need the Plumed.h, Plumed.inc, and Plumed.cmake files if test -n "$include" ; then test -n "$quiet" || echo "Including Plumed.h, Plumed.inc, and Plumed.cmake ($mode mode)" echo "#include \"$PLUMED_INCLUDEDIR/$PLUMED_PROGRAM_NAME/wrapper/Plumed.h\"" > Plumed.h @@ -356,10 +357,14 @@ case "$action" in fi ;; (save) - if [ ! -e Plumed.h -o ! -e Plumed.inc -o ! -e Plumed.cmake ] - then - echo "ERROR: I cannot find Plumed.h, Plumed.inc, and Plumed.cmake files. You have likely not patched yet." - exit 1 + if [[ -z $PLUMED_ONLY_PATCH ]]; then + #PLUMED_ONLY_PATCH is special for programs that have an embedded plumed integration + #and do not need the Plumed.h, Plumed.inc, and Plumed.cmake files + if [ ! -e Plumed.h ] || [ ! -e Plumed.inc ] || [ ! -e Plumed.cmake ] + then + echo "ERROR: I cannot find Plumed.h, Plumed.inc, and Plumed.cmake files. You have likely not patched yet." + exit 1 + fi fi PREPLUMED=$(find . -name "*.preplumed" | sort) for file in $PLUMED_PREPLUMED_IGNORE; From cf6b758e48bad7f45df078207beaf8fdc192e2dc Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Wed, 29 Jan 2025 17:08:30 +0100 Subject: [PATCH 5/7] Now plumed is forced on --- .../cmake/gmxManagePlumed.cmake | 82 +++++++++++++++++++ .../cmake/gmxManagePlumed.cmake.preplumed | 82 +++++++++++++++++++ 2 files changed, 164 insertions(+) create mode 100644 patches/gromacs-2025.0.diff/cmake/gmxManagePlumed.cmake create mode 100644 patches/gromacs-2025.0.diff/cmake/gmxManagePlumed.cmake.preplumed diff --git a/patches/gromacs-2025.0.diff/cmake/gmxManagePlumed.cmake b/patches/gromacs-2025.0.diff/cmake/gmxManagePlumed.cmake new file mode 100644 index 0000000000..6d092ef4ce --- /dev/null +++ b/patches/gromacs-2025.0.diff/cmake/gmxManagePlumed.cmake @@ -0,0 +1,82 @@ +# +# This file is part of the GROMACS molecular simulation package. +# +# Copyright 2024- 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. + + + + + + + +# Builds the interface to plumed and add the linkage to libPlumed + +gmx_option_multichoice(GMX_USE_PLUMED + "Build the PLUMED wrapper with GROMACS" + ON + ON) +mark_as_advanced(GMX_USE_PLUMED) + +function(gmx_manage_plumed) + # Create a link target, leave it empty if the plumed option is not active + add_library(plumedgmx INTERFACE) + set (GMX_PLUMED_ACTIVE OFF CACHE INTERNAL "Cache entry for PLUMED activation") + if(WIN32) + if(GMX_USE_PLUMED STREQUAL "ON") + message(FATAL_ERROR "PLUMED is not supported on Windows. Reconfigure with -DGMX_USE_PLUMED=OFF.") + endif() + elseif(NOT GMX_USE_PLUMED STREQUAL "OFF") + include(CMakePushCheckState) + cmake_push_check_state(RESET) + set(CMAKE_REQUIRED_LIBRARIES ${CMAKE_DL_LIBS}) + include(CheckFunctionExists) + check_function_exists(dlopen HAVE_DLOPEN) + cmake_pop_check_state() + if(HAVE_DLOPEN) + # Plumed.h compiled in c++ mode creates a fully inlined interface + # So we just need to activate the directory in applied_forces + set(PLUMED_DIR "${CMAKE_SOURCE_DIR}/src/external/plumed") + target_link_libraries( plumedgmx INTERFACE ${CMAKE_DL_LIBS} ) + # The plumedgmx already exists, now we set it up: + target_include_directories(plumedgmx SYSTEM INTERFACE $) + set (GMX_PLUMED_ACTIVE ON CACHE INTERNAL "Cache entry for PLUMED activation") + else() + if(GMX_USE_PLUMED STREQUAL "ON") + message(FATAL_ERROR "PLUMED needs dlopen or anything equivalent. Reconfigure with -DGMX_USE_PLUMED=OFF.") + else() # "AUTO" + message(STATUS "PLUMED needs dlopen or anything equivalent. Disabling support.") + endif() + + endif() + + endif() + +endfunction() diff --git a/patches/gromacs-2025.0.diff/cmake/gmxManagePlumed.cmake.preplumed b/patches/gromacs-2025.0.diff/cmake/gmxManagePlumed.cmake.preplumed new file mode 100644 index 0000000000..beaa31af64 --- /dev/null +++ b/patches/gromacs-2025.0.diff/cmake/gmxManagePlumed.cmake.preplumed @@ -0,0 +1,82 @@ +# +# This file is part of the GROMACS molecular simulation package. +# +# Copyright 2024- 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. + + + + + + + +# Builds the interface to plumed and add the linkage to libPlumed + +gmx_option_multichoice(GMX_USE_PLUMED + "Build the PLUMED wrapper with GROMACS" + AUTO + AUTO ON OFF) +mark_as_advanced(GMX_USE_PLUMED) + +function(gmx_manage_plumed) + # Create a link target, leave it empty if the plumed option is not active + add_library(plumedgmx INTERFACE) + set (GMX_PLUMED_ACTIVE OFF CACHE INTERNAL "Cache entry for PLUMED activation") + if(WIN32) + if(GMX_USE_PLUMED STREQUAL "ON") + message(FATAL_ERROR "PLUMED is not supported on Windows. Reconfigure with -DGMX_USE_PLUMED=OFF.") + endif() + elseif(NOT GMX_USE_PLUMED STREQUAL "OFF") + include(CMakePushCheckState) + cmake_push_check_state(RESET) + set(CMAKE_REQUIRED_LIBRARIES ${CMAKE_DL_LIBS}) + include(CheckFunctionExists) + check_function_exists(dlopen HAVE_DLOPEN) + cmake_pop_check_state() + if(HAVE_DLOPEN) + # Plumed.h compiled in c++ mode creates a fully inlined interface + # So we just need to activate the directory in applied_forces + set(PLUMED_DIR "${CMAKE_SOURCE_DIR}/src/external/plumed") + target_link_libraries( plumedgmx INTERFACE ${CMAKE_DL_LIBS} ) + # The plumedgmx already exists, now we set it up: + target_include_directories(plumedgmx SYSTEM INTERFACE $) + set (GMX_PLUMED_ACTIVE ON CACHE INTERNAL "Cache entry for PLUMED activation") + else() + if(GMX_USE_PLUMED STREQUAL "ON") + message(FATAL_ERROR "PLUMED needs dlopen or anything equivalent. Reconfigure with -DGMX_USE_PLUMED=OFF.") + else() # "AUTO" + message(STATUS "PLUMED needs dlopen or anything equivalent. Disabling support.") + endif() + + endif() + + endif() + +endfunction() From 869a08ab1b8c8730a085ad9779df004e601ede3c Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Thu, 30 Jan 2025 10:57:31 +0100 Subject: [PATCH 6/7] now the patch can move "extra files" --- patches/gromacs-2025.0.config | 4 +- .../external/plumed/PlumedInclude.h.preplumed | 0 .../plumed/PlumedKernel.cpp.preplumed | 0 .../plumed/PlumedOutside.cpp.preplumed | 0 .../plumed/PlumedOutside.h.preplumed | 0 patches/patch.sh | 42 +++++++++++++++---- 6 files changed, 37 insertions(+), 9 deletions(-) delete mode 100644 patches/gromacs-2025.0.diff/src/external/plumed/PlumedInclude.h.preplumed delete mode 100644 patches/gromacs-2025.0.diff/src/external/plumed/PlumedKernel.cpp.preplumed delete mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.cpp.preplumed delete mode 100644 patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.h.preplumed diff --git a/patches/gromacs-2025.0.config b/patches/gromacs-2025.0.config index 9d3bce3bcb..485beae4a4 100644 --- a/patches/gromacs-2025.0.config +++ b/patches/gromacs-2025.0.config @@ -33,4 +33,6 @@ plumed_before_patch(){ plumed_patch_info } -PLUMED_ONLY_PATCH=true +PLUMED_PATCH_NO_INCLUDE=true + +PLUMED_PATCH_EXTRA_FILES="src/external/plumed/PlumedInclude.h src/external/plumed/PlumedKernel.cpp src/gromacs/applied_forces/plumed/PlumedOutside.cpp src/gromacs/applied_forces/plumed/PlumedOutside.h" diff --git a/patches/gromacs-2025.0.diff/src/external/plumed/PlumedInclude.h.preplumed b/patches/gromacs-2025.0.diff/src/external/plumed/PlumedInclude.h.preplumed deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/patches/gromacs-2025.0.diff/src/external/plumed/PlumedKernel.cpp.preplumed b/patches/gromacs-2025.0.diff/src/external/plumed/PlumedKernel.cpp.preplumed deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.cpp.preplumed b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.cpp.preplumed deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.h.preplumed b/patches/gromacs-2025.0.diff/src/gromacs/applied_forces/plumed/PlumedOutside.h.preplumed deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/patches/patch.sh b/patches/patch.sh index aae5a84432..7f9aab5fd8 100755 --- a/patches/patch.sh +++ b/patches/patch.sh @@ -286,8 +286,8 @@ case "$action" in plumed_before_patch fi - if [[ -z $PLUMED_ONLY_PATCH ]]; then - #PLUMED_ONLY_PATCH is special for programs that have an embedded plumed integration + if [[ -z $PLUMED_PATCH_NO_INCLUDE ]]; then + #PLUMED_PATCH_NO_INCLUDE is special for programs that have an embedded plumed integration #and do not need the Plumed.h, Plumed.inc, and Plumed.cmake files if test -n "$include" ; then test -n "$quiet" || echo "Including Plumed.h, Plumed.inc, and Plumed.cmake ($mode mode)" @@ -302,6 +302,11 @@ case "$action" in fi fi + if [[ -n $PLUMED_PATCH_EXTRA_FILES ]]; then + for file in $PLUMED_PATCH_EXTRA_FILES; do + cp "${diff}/$file" "$file" + done + fi if [ -d "$diff" ]; then test -n "$quiet" || echo "Patching with on-the-fly diff from stored originals" PREPLUMED=$(cd "$diff" ; find . -name "*.preplumed" | sort) @@ -357,8 +362,8 @@ case "$action" in fi ;; (save) - if [[ -z $PLUMED_ONLY_PATCH ]]; then - #PLUMED_ONLY_PATCH is special for programs that have an embedded plumed integration + if [[ -z $PLUMED_PATCH_NO_INCLUDE ]]; then + #PLUMED_PATCH_NO_INCLUDE is special for programs that have an embedded plumed integration #and do not need the Plumed.h, Plumed.inc, and Plumed.cmake files if [ ! -e Plumed.h ] || [ ! -e Plumed.inc ] || [ ! -e Plumed.cmake ] then @@ -389,6 +394,11 @@ case "$action" in test -n "$quiet" || echo "Saving your changes to $diff" test -n "$quiet" || echo "Preplumed files:" test -n "$quiet" || echo "$PREPLUMED" + if [[ -n $PLUMED_PATCH_EXTRA_FILES ]]; then + echo "Extra files:" + echo "$PLUMED_PATCH_EXTRA_FILES" + fi + if [ -d "$diff" ] && [ -z "$save_originals" ]; then echo "This patch uses the dir format (originals are saved)" echo "Are you sure you want to save the single diff?" @@ -421,14 +431,24 @@ case "$action" in cp "$file" "$diff/$file" cp "$bckfile" "$diff/$bckfile" else - echo "patch -u -l -b -F 5 -N --suffix=.preplumed \"${file}\" << \\EOF_EOF" >> "$diff" - diff -U 5 "${bckfile}" "$file" --label="$bckfile" --label="$file" >> "$diff" - echo "EOF_EOF" >> "$diff" + { + echo "patch -u -l -b -F 5 -N --suffix=.preplumed \"${file}\" << \\EOF_EOF" + diff -U 5 "${bckfile}" "$file" --label="$bckfile" --label="$file" + echo "EOF_EOF" + } >> "$diff" fi else echo "ERROR: File $file is missing" fi done + if [[ -n $PLUMED_PATCH_EXTRA_FILES ]]; then + test -n "$quiet" || echo "saving extra files" + for file in $PLUMED_PATCH_EXTRA_FILES; do + xx=${diff}/$file + mkdir -p "${xx%/*}" + cp "$file" "${diff}/$file" + done + fi cat </dev/null ; then test -n "$quiet" || echo "Executing plumed_after_revert function" plumed_after_revert From 2fda83870c6851a09025a56d1c621fc449afd4d3 Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Fri, 31 Jan 2025 14:23:37 +0100 Subject: [PATCH 7/7] clang-formatted --- .../applied_forces/plumed/plumedMDModule.cpp | 7 +++---- .../applied_forces/plumed/plumedOptions.h | 18 +++++++++--------- .../plumed/plumedforceprovider.cpp | 3 +-- 3 files changed, 13 insertions(+), 15 deletions(-) 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 c56c259a6a..978fb5bc01 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 @@ -59,7 +59,6 @@ class gmx_multisim_t; namespace gmx { - namespace { @@ -96,10 +95,10 @@ class PlumedMDModule final : public IMDModule { this->options_.setPlumedFile(plumedFilename.plumedFilename_); }); // Retrieve the Multisim options - notifier->simulationSetupNotifier_.subscribe( - [this](const gmx_multisim_t* ms) { this->options_.setMultisim(ms); }); + notifier->simulationSetupNotifier_.subscribe([this](const gmx_multisim_t* ms) + { this->options_.setMultisim(ms); }); + - // Access the temperature if it is constant during the simulation notifier->simulationSetupNotifier_.subscribe( [this](const EnsembleTemperature& ensembleT) 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 097aaa1bd5..cbfaaa1afb 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 @@ -58,14 +58,14 @@ struct EnsembleTemperature; struct PlumedOptions { - std::string plumedFile_; - int natoms_; - const t_commrec* cr_; - const gmx_multisim_t* ms_; - real simulationTimeStep_; - std::optional ensembleTemperature_{}; - StartingBehavior startingBehavior_{}; - bool active_{ false }; + std::string plumedFile_; + int natoms_; + const t_commrec* cr_; + const gmx_multisim_t* ms_; + real simulationTimeStep_; + std::optional ensembleTemperature_{}; + StartingBehavior startingBehavior_{}; + bool active_{ false }; }; class PlumedOptionProvider @@ -99,7 +99,7 @@ class PlumedOptionProvider /*! @brief Sets the address to the communication record object * @param cr the Communication Record object */ - void setComm(const t_commrec& cr); + void setComm(const t_commrec& cr); /*! @brief Sets the address to the multisimulation object * @param ms the address to the multisimulation object */ 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 6b5a0fbae8..85cef6e051 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 @@ -48,10 +48,10 @@ #include "gromacs/math/units.h" #include "gromacs/mdlib/gmx_omp_nthreads.h" #include "gromacs/mdrunutility/handlerestart.h" +#include "gromacs/mdrunutility/multisim.h" #include "gromacs/mdtypes/commrec.h" #include "gromacs/mdtypes/enerdata.h" #include "gromacs/mdtypes/forceoutput.h" -#include "gromacs/mdrunutility/multisim.h" #include "gromacs/utility/exceptions.h" #include "PlumedOutside.h" @@ -119,7 +119,6 @@ try : plumed_(std::make_unique()) } */ - if (isMultiSim(options.ms_)) {