Implement Berendsen and c-rescale pressure coupling for modular simulator
authorPascal Merz <pascal.merz@me.com>
Mon, 24 May 2021 04:27:50 +0000 (04:27 +0000)
committerMark Abraham <mark.j.abraham@gmail.com>
Mon, 24 May 2021 04:27:50 +0000 (04:27 +0000)
src/gromacs/modularsimulator/firstorderpressurecoupling.cpp [new file with mode: 0644]
src/gromacs/modularsimulator/firstorderpressurecoupling.h [new file with mode: 0644]
src/gromacs/modularsimulator/modularsimulator.cpp
src/gromacs/modularsimulator/modularsimulatorinterfaces.h
src/gromacs/modularsimulator/simulatoralgorithm.cpp
src/gromacs/modularsimulator/velocityscalingtemperaturecoupling.h
src/programs/mdrun/tests/simulator.cpp

diff --git a/src/gromacs/modularsimulator/firstorderpressurecoupling.cpp b/src/gromacs/modularsimulator/firstorderpressurecoupling.cpp
new file mode 100644 (file)
index 0000000..c1ecee5
--- /dev/null
@@ -0,0 +1,258 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * 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
+ * http://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 http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief Defines the element performing first-order pressure coupling for the modular simulator
+ *
+ * \author Pascal Merz <pascal.merz@me.com>
+ * \ingroup module_modularsimulator
+ */
+
+#include "gmxpre.h"
+
+#include "firstorderpressurecoupling.h"
+
+#include "gromacs/domdec/domdec_network.h"
+#include "gromacs/mdlib/coupling.h"
+#include "gromacs/mdlib/mdatoms.h"
+#include "gromacs/mdlib/stat.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/mdatom.h"
+#include "gromacs/pbcutil/boxutilities.h"
+
+#include "energydata.h"
+#include "statepropagatordata.h"
+#include "simulatoralgorithm.h"
+
+namespace gmx
+{
+
+template<PressureCoupling pressureCouplingType>
+void FirstOrderPressureCoupling::calculateScalingMatrix(Step step)
+{
+    const auto* pressure         = energyData_->pressure(step);
+    const auto* forceVirial      = energyData_->forceVirial(step);
+    const auto* constraintVirial = energyData_->constraintVirial(step);
+    const auto* box              = statePropagatorData_->constBox();
+
+    previousStepConservedEnergyContribution_ = conservedEnergyContribution_;
+    pressureCouplingCalculateScalingMatrix<pressureCouplingType>(fplog_,
+                                                                 step,
+                                                                 inputrec_,
+                                                                 couplingTimeStep_,
+                                                                 pressure,
+                                                                 box,
+                                                                 forceVirial,
+                                                                 constraintVirial,
+                                                                 boxScalingMatrix_,
+                                                                 &conservedEnergyContribution_);
+    conservedEnergyContributionStep_ = step;
+}
+
+template<PressureCoupling pressureCouplingType>
+void FirstOrderPressureCoupling::scaleBoxAndCoordinates()
+{
+    auto*          box       = statePropagatorData_->box();
+    auto           positions = statePropagatorData_->positionsView().unpaddedArrayRef();
+    ArrayRef<RVec> velocities;
+    if (pressureCouplingType == PressureCoupling::CRescale)
+    {
+        velocities = statePropagatorData_->velocitiesView().unpaddedArrayRef();
+    }
+    // Freeze groups are not implemented
+    ArrayRef<const unsigned short> cFreeze;
+    // Coordinates are always scaled except for GPU update (not implemented currently)
+    const bool scaleCoordinates = true;
+    // Atom range
+    const int startAtom = 0;
+    const int numAtoms  = mdAtoms_->mdatoms()->homenr;
+
+    pressureCouplingScaleBoxAndCoordinates<pressureCouplingType>(
+            inputrec_, boxScalingMatrix_, box, boxRel_, startAtom, numAtoms, positions, velocities, cFreeze, nrnb_, scaleCoordinates);
+}
+
+void FirstOrderPressureCoupling::scheduleTask(Step step, Time /*unused*/, const RegisterRunFunction& registerRunFunction)
+{
+    if (do_per_step(step + couplingFrequency_ + couplingOffset_, couplingFrequency_))
+    {
+        if (pressureCouplingType_ == PressureCoupling::Berendsen)
+        {
+            registerRunFunction([this, step]() {
+                calculateScalingMatrix<PressureCoupling::Berendsen>(step);
+                scaleBoxAndCoordinates<PressureCoupling::Berendsen>();
+            });
+        }
+        else if (pressureCouplingType_ == PressureCoupling::CRescale)
+        {
+            registerRunFunction([this, step]() {
+                calculateScalingMatrix<PressureCoupling::CRescale>(step);
+                scaleBoxAndCoordinates<PressureCoupling::CRescale>();
+            });
+        }
+    }
+}
+
+void FirstOrderPressureCoupling::elementSetup()
+{
+    if (inputrecPreserveShape(inputrec_))
+    {
+        auto*     box  = statePropagatorData_->box();
+        const int ndim = inputrec_->epct == PressureCouplingType::SemiIsotropic ? 2 : 3;
+        do_box_rel(ndim, inputrec_->deform, boxRel_, box, true);
+    }
+}
+
+real FirstOrderPressureCoupling::conservedEnergyContribution(Step step)
+{
+    if (step == conservedEnergyContributionStep_
+        && reportPreviousStepConservedEnergy_ == ReportPreviousStepConservedEnergy::Yes)
+    {
+        return previousStepConservedEnergyContribution_;
+    }
+    return conservedEnergyContribution_;
+}
+
+namespace
+{
+/*!
+ * \brief Enum describing the contents FirstOrderPressureCoupling writes to modular checkpoint
+ *
+ * When changing the checkpoint content, add a new element just above Count, and adjust the
+ * checkpoint functionality.
+ */
+enum class CheckpointVersion
+{
+    Base, //!< First version of modular checkpointing
+    Count //!< Number of entries. Add new versions right above this!
+};
+constexpr auto c_currentVersion = CheckpointVersion(int(CheckpointVersion::Count) - 1);
+} // namespace
+
+template<CheckpointDataOperation operation>
+void FirstOrderPressureCoupling::doCheckpointData(CheckpointData<operation>* checkpointData)
+{
+    checkpointVersion(checkpointData, "FirstOrderPressureCoupling version", c_currentVersion);
+
+    checkpointData->scalar("conserved energy contribution", &conservedEnergyContribution_);
+    checkpointData->scalar("conserved energy step", &conservedEnergyContributionStep_);
+    checkpointData->tensor("relative box vector", boxRel_);
+}
+
+void FirstOrderPressureCoupling::saveCheckpointState(std::optional<WriteCheckpointData> checkpointData,
+                                                     const t_commrec*                   cr)
+{
+    if (MASTER(cr))
+    {
+        doCheckpointData<CheckpointDataOperation::Write>(&checkpointData.value());
+    }
+}
+
+void FirstOrderPressureCoupling::restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData,
+                                                        const t_commrec*                  cr)
+{
+    if (MASTER(cr))
+    {
+        doCheckpointData<CheckpointDataOperation::Read>(&checkpointData.value());
+    }
+    if (DOMAINDECOMP(cr))
+    {
+        dd_bcast(cr->dd, sizeof(conservedEnergyContribution_), &conservedEnergyContribution_);
+        dd_bcast(cr->dd, sizeof(conservedEnergyContributionStep_), &conservedEnergyContributionStep_);
+        dd_bcast(cr->dd, sizeof(boxRel_), boxRel_);
+    }
+}
+
+const std::string& FirstOrderPressureCoupling::clientID()
+{
+    return identifier_;
+}
+
+FirstOrderPressureCoupling::FirstOrderPressureCoupling(int                  couplingFrequency,
+                                                       int                  couplingOffset,
+                                                       real                 couplingTimeStep,
+                                                       StatePropagatorData* statePropagatorData,
+                                                       EnergyData*          energyData,
+                                                       FILE*                fplog,
+                                                       const t_inputrec*    inputrec,
+                                                       const MDAtoms*       mdAtoms,
+                                                       t_nrnb*              nrnb,
+                                                       ReportPreviousStepConservedEnergy reportPreviousStepConservedEnergy) :
+    pressureCouplingType_(inputrec->epc),
+    couplingTimeStep_(couplingTimeStep),
+    couplingFrequency_(couplingFrequency),
+    couplingOffset_(couplingOffset),
+    boxScalingMatrix_{ { 0 } },
+    boxRel_{ { 0 } },
+    conservedEnergyContribution_(0),
+    previousStepConservedEnergyContribution_(0),
+    conservedEnergyContributionStep_(-1),
+    reportPreviousStepConservedEnergy_(reportPreviousStepConservedEnergy),
+    statePropagatorData_(statePropagatorData),
+    energyData_(energyData),
+    fplog_(fplog),
+    inputrec_(inputrec),
+    mdAtoms_(mdAtoms),
+    nrnb_(nrnb),
+    identifier_("FirstOrderPressureCoupling-" + std::string(enumValueToString(pressureCouplingType_)))
+{
+    energyData->addConservedEnergyContribution(
+            [this](Step step, Time /*unused*/) { return conservedEnergyContribution(step); });
+}
+
+ISimulatorElement* FirstOrderPressureCoupling::getElementPointerImpl(
+        LegacySimulatorData*                    legacySimulatorData,
+        ModularSimulatorAlgorithmBuilderHelper* builderHelper,
+        StatePropagatorData*                    statePropagatorData,
+        EnergyData*                             energyData,
+        FreeEnergyPerturbationData gmx_unused* freeEnergyPerturbationData,
+        GlobalCommunicationHelper gmx_unused* globalCommunicationHelper,
+        int                                   offset,
+        ReportPreviousStepConservedEnergy     reportPreviousStepConservedEnergy)
+{
+    return builderHelper->storeElement(std::make_unique<FirstOrderPressureCoupling>(
+            legacySimulatorData->inputrec->nstpcouple,
+            offset,
+            legacySimulatorData->inputrec->delta_t * legacySimulatorData->inputrec->nstpcouple,
+            statePropagatorData,
+            energyData,
+            legacySimulatorData->fplog,
+            legacySimulatorData->inputrec,
+            legacySimulatorData->mdAtoms,
+            legacySimulatorData->nrnb,
+            reportPreviousStepConservedEnergy));
+}
+
+} // namespace gmx
diff --git a/src/gromacs/modularsimulator/firstorderpressurecoupling.h b/src/gromacs/modularsimulator/firstorderpressurecoupling.h
new file mode 100644 (file)
index 0000000..78c9e9a
--- /dev/null
@@ -0,0 +1,178 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * 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
+ * http://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 http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief Declares the element performing first-order pressure coupling for the modular simulator
+ *
+ * \author Pascal Merz <pascal.merz@me.com>
+ * \ingroup module_modularsimulator
+ *
+ * This header is only used within the modular simulator module
+ */
+
+#ifndef GMX_MODULARSIMULATOR_FIRSTORDERPRESSURECOUPLING_H
+#define GMX_MODULARSIMULATOR_FIRSTORDERPRESSURECOUPLING_H
+
+#include "modularsimulatorinterfaces.h"
+
+struct t_inputrec;
+struct t_nrnb;
+
+enum class PressureCoupling;
+
+namespace gmx
+{
+class EnergyData;
+class FreeEnergyPerturbationData;
+class GlobalCommunicationHelper;
+class LegacySimulatorData;
+class MDAtoms;
+class ModularSimulatorAlgorithmBuilderHelper;
+class StatePropagatorData;
+
+/*! \internal
+ * \brief Element performing first-order pressure coupling
+ *
+ * This element implements the first-order pressure coupling algorithms
+ * (Berendesen, c-rescale).
+ */
+class FirstOrderPressureCoupling final : public ISimulatorElement, public ICheckpointHelperClient
+{
+public:
+    //! Constructor
+    FirstOrderPressureCoupling(int                               couplingFrequency,
+                               int                               couplingOffset,
+                               real                              couplingTimeStep,
+                               StatePropagatorData*              statePropagatorData,
+                               EnergyData*                       energyData,
+                               FILE*                             fplog,
+                               const t_inputrec*                 inputrec,
+                               const MDAtoms*                    mdAtoms,
+                               t_nrnb*                           nrnb,
+                               ReportPreviousStepConservedEnergy reportPreviousStepConservedEnergy);
+
+    void scheduleTask(Step step, Time time, const RegisterRunFunction& function) override;
+    //! Setup - initialize relative box matrix
+    void elementSetup() override;
+    //! No teardown needed
+    void elementTeardown() override{};
+
+    //! ICheckpointHelperClient write checkpoint implementation
+    void saveCheckpointState(std::optional<WriteCheckpointData> checkpointData, const t_commrec* cr) override;
+    //! ICheckpointHelperClient read checkpoint implementation
+    void restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData, const t_commrec* cr) override;
+    //! ICheckpointHelperClient key implementation
+    const std::string& clientID() override;
+
+    /*! \brief Factory method implementation
+     *
+     * \param legacySimulatorData  Pointer allowing access to simulator level data
+     * \param builderHelper  ModularSimulatorAlgorithmBuilder helper object
+     * \param statePropagatorData  Pointer to the \c StatePropagatorData object
+     * \param energyData  Pointer to the \c EnergyData object
+     * \param freeEnergyPerturbationData  Pointer to the \c FreeEnergyPerturbationData object
+     * \param globalCommunicationHelper  Pointer to the \c GlobalCommunicationHelper object
+     * \param offset  The step offset at which the barostat is applied
+     * \param reportPreviousStepConservedEnergy  Report the previous or the current step conserved energy
+     *
+     * \return  Pointer to the element to be added. Element needs to have been stored using \c storeElement
+     */
+    static ISimulatorElement*
+    getElementPointerImpl(LegacySimulatorData*                    legacySimulatorData,
+                          ModularSimulatorAlgorithmBuilderHelper* builderHelper,
+                          StatePropagatorData*                    statePropagatorData,
+                          EnergyData*                             energyData,
+                          FreeEnergyPerturbationData gmx_unused* freeEnergyPerturbationData,
+                          GlobalCommunicationHelper gmx_unused* globalCommunicationHelper,
+                          int                                   offset,
+                          ReportPreviousStepConservedEnergy     reportPreviousStepConservedEnergy);
+
+private:
+    //! Calculate the scaling matrix
+    template<PressureCoupling pressureCouplingType>
+    void calculateScalingMatrix(Step step);
+    //! Scale the box and coordinates according to the current scaling matrix
+    template<PressureCoupling pressureCouplingType>
+    void scaleBoxAndCoordinates();
+    //! Helper function returning the conserved energy contribution
+    real conservedEnergyContribution(Step step);
+
+    //! The pressure coupling type required
+    const PressureCoupling pressureCouplingType_;
+    //! The coupling time step (simulation time step x coupling frequency)
+    const real couplingTimeStep_;
+    //! The frequency at which pressure coupling is applied
+    const int couplingFrequency_;
+    //! The offset at which pressure coupling is applied
+    const int couplingOffset_;
+
+    //! The current box scaling matrix
+    tensor boxScalingMatrix_;
+    //! Relative box shape
+    tensor boxRel_;
+    //! Contribution to the conserved energy
+    double conservedEnergyContribution_;
+    //! Contribution to the conserved energy
+    double previousStepConservedEnergyContribution_;
+    //! Current step of the conserved energy contribution
+    Step conservedEnergyContributionStep_;
+    //! Whether we're reporting current step or previous step conserved energy
+    const ReportPreviousStepConservedEnergy reportPreviousStepConservedEnergy_;
+
+    // TODO: Clarify relationship to data objects and find a more robust alternative to raw pointers (#3583)
+    //! Pointer to the micro state
+    StatePropagatorData* statePropagatorData_;
+    //! Pointer to the energy data
+    EnergyData* energyData_;
+
+    // Access to ISimulator data
+    //! Handles logging.
+    FILE* fplog_;
+    //! Contains user input mdp options.
+    const t_inputrec* inputrec_;
+    //! Atom parameters for this domain.
+    const MDAtoms* mdAtoms_;
+    //! Manages flop accounting.
+    t_nrnb* nrnb_;
+
+    //! CheckpointHelper identifier
+    const std::string identifier_;
+    //! Helper function to read from / write to CheckpointData
+    template<CheckpointDataOperation operation>
+    void doCheckpointData(CheckpointData<operation>* checkpointData);
+};
+
+} // namespace gmx
+
+#endif // GMX_MODULARSIMULATOR_FIRSTORDERPRESSURECOUPLING_H
index 2431e18aa2e759101ada4a27a1c968354bfe00b2..f1fdd17d8a3e941323f0789c77fb6eb7caf89f21 100644 (file)
@@ -78,6 +78,7 @@
 
 #include "computeglobalselement.h"
 #include "constraintelement.h"
+#include "firstorderpressurecoupling.h"
 #include "forceelement.h"
 #include "parrinellorahmanbarostat.h"
 #include "simulatoralgorithm.h"
@@ -131,6 +132,11 @@ void ModularSimulator::addIntegrationElements(ModularSimulatorAlgorithmBuilder*
         {
             builder->add<ParrinelloRahmanBarostat>(Offset(-1), PropagatorTag("LeapFrogPropagator"));
         }
+        else if (legacySimulatorData_->inputrec->epc == PressureCoupling::Berendsen
+                 || legacySimulatorData_->inputrec->epc == PressureCoupling::CRescale)
+        {
+            builder->add<FirstOrderPressureCoupling>(0, ReportPreviousStepConservedEnergy::No);
+        }
     }
     else if (legacySimulatorData_->inputrec->eI == IntegrationAlgorithm::VV)
     {
@@ -165,6 +171,11 @@ void ModularSimulator::addIntegrationElements(ModularSimulatorAlgorithmBuilder*
         {
             builder->add<ParrinelloRahmanBarostat>(Offset(-1), PropagatorTag("VelocityHalfStep"));
         }
+        else if (legacySimulatorData_->inputrec->epc == PressureCoupling::Berendsen
+                 || legacySimulatorData_->inputrec->epc == PressureCoupling::CRescale)
+        {
+            builder->add<FirstOrderPressureCoupling>(0, ReportPreviousStepConservedEnergy::Yes);
+        }
     }
     else
     {
@@ -237,11 +248,13 @@ bool ModularSimulator::isInputCompatible(bool                             exitOn
                                                          && inputrec->eI == IntegrationAlgorithm::MD),
                                              "Only v-rescale, Berendsen and Nose-Hoover "
                                              "thermostats are supported by the modular simulator.");
-    isInputCompatible =
-            isInputCompatible
-            && conditionalAssert(
-                    inputrec->epc == PressureCoupling::No || inputrec->epc == PressureCoupling::ParrinelloRahman,
-                    "Only Parrinello-Rahman barostat is supported by the modular simulator.");
+    isInputCompatible = isInputCompatible
+                        && conditionalAssert(inputrec->epc == PressureCoupling::No
+                                                     || inputrec->epc == PressureCoupling::ParrinelloRahman
+                                                     || inputrec->epc == PressureCoupling::Berendsen
+                                                     || inputrec->epc == PressureCoupling::CRescale,
+                                             "Only Parrinello-Rahman, Berendsen and C-Rescale "
+                                             "barostats are supported by the modular simulator.");
     isInputCompatible =
             isInputCompatible
             && conditionalAssert(
index 3dcbb6b03a2ec27c2962087699de1d62431a6301..e6b4cd6eabeb72aef432a0a422dc04560c9f8c8c 100644 (file)
@@ -553,6 +553,14 @@ struct PropagatorConnection
     std::function<PropagatorCallback()> getPRScalingCallback;
 };
 
+//! Enum describing whether an element is reporting conserved energy from the previous step
+enum class ReportPreviousStepConservedEnergy
+{
+    Yes,
+    No,
+    Count
+};
+
 //! /}
 } // namespace gmx
 
index 6628e57009cda0084515daacc70b787123437a56..9808bbe530b84996bbdc9c7cda348cba5316fe76 100644 (file)
@@ -77,6 +77,7 @@
 #include "checkpointhelper.h"
 #include "domdechelper.h"
 #include "energydata.h"
+#include "firstorderpressurecoupling.h"
 #include "freeenergyperturbationdata.h"
 #include "modularsimulator.h"
 #include "pmeloadbalancehelper.h"
index 489c9262adc09042ad4353af3ad57099ebdb8f60..2cd8f58541e62134d3b3a2a50c5792e8257ba5a9 100644 (file)
@@ -67,14 +67,6 @@ enum class UseFullStepKE
     Count
 };
 
-//! Enum describing whether the thermostat is reporting conserved energy from the previous step
-enum class ReportPreviousStepConservedEnergy
-{
-    Yes,
-    No,
-    Count
-};
-
 /*! \internal
  * \ingroup module_modularsimulator
  * \brief Element implementing the a velocity-scaling thermostat
index a1889c6c35274645fd96e4c3c22212f3f4ea3611..8e3c4791d7214fbff56ec9c947365b32c50028f7 100644 (file)
@@ -90,6 +90,8 @@ TEST_P(SimulatorComparisonTest, WithinTolerances)
     const auto& pcoupling           = std::get<3>(mdpParams);
     const auto& environmentVariable = std::get<1>(params);
 
+    int maxNumWarnings = 0;
+
     // TODO At some point we should also test PME-only ranks.
     const int numRanksAvailable = getNumberOfTestMpiRanks();
     if (!isNumberOfPpRanksSupported(simulationName, numRanksAvailable))
@@ -110,6 +112,20 @@ TEST_P(SimulatorComparisonTest, WithinTolerances)
         return;
     }
 
+    if (tcoupling == "nose-hoover" && pcoupling == "berendsen")
+    {
+        if (integrator == "md-vv")
+        {
+            // Combination not allowed by legacy do_md.
+            return;
+        }
+        else
+        {
+            // "Using Berendsen pressure coupling invalidates the true ensemble for the thermostat"
+            maxNumWarnings++;
+        }
+    }
+
     const std::string envVariableModSimOn  = "GMX_USE_MODULAR_SIMULATOR";
     const std::string envVariableModSimOff = "GMX_DISABLE_MODULAR_SIMULATOR";
 
@@ -191,7 +207,7 @@ TEST_P(SimulatorComparisonTest, WithinTolerances)
     runner_.tprFileName_ = fileManager_.getTemporaryFilePath("sim.tpr");
     runner_.useTopGroAndNdxFromDatabase(simulationName);
     runner_.useStringAsMdpFile(prepareMdpFileContents(mdpFieldValues));
-    runGrompp(&runner_);
+    runGrompp(&runner_, { SimulationOptionTuple("-maxwarn", std::to_string(maxNumWarnings)) });
 
     // Backup current state of both environment variables and unset them
     const char* environmentVariableBackupOn  = getenv(envVariableModSimOn.c_str());
@@ -242,16 +258,17 @@ INSTANTIATE_TEST_CASE_P(
         ::testing::Combine(::testing::Combine(::testing::Values("argon12", "tip3p5"),
                                               ::testing::Values("md-vv"),
                                               ::testing::Values("no", "v-rescale", "berendsen"),
-                                              ::testing::Values("no")),
+                                              ::testing::Values("no", "berendsen", "c-rescale")),
                            ::testing::Values("GMX_DISABLE_MODULAR_SIMULATOR")));
 INSTANTIATE_TEST_CASE_P(
         SimulatorsAreEquivalentDefaultLegacy,
         SimulatorComparisonTest,
         ::testing::Combine(
-                ::testing::Combine(::testing::Values("argon12", "tip3p5"),
-                                   ::testing::Values("md"),
-                                   ::testing::Values("no", "v-rescale", "berendsen", "nose-hoover"),
-                                   ::testing::Values("no", "Parrinello-Rahman")),
+                ::testing::Combine(
+                        ::testing::Values("argon12", "tip3p5"),
+                        ::testing::Values("md"),
+                        ::testing::Values("no", "v-rescale", "berendsen", "nose-hoover"),
+                        ::testing::Values("no", "Parrinello-Rahman", "berendsen", "c-rescale")),
                 ::testing::Values("GMX_USE_MODULAR_SIMULATOR")));
 #else
 INSTANTIATE_TEST_CASE_P(
@@ -260,16 +277,17 @@ INSTANTIATE_TEST_CASE_P(
         ::testing::Combine(::testing::Combine(::testing::Values("argon12", "tip3p5"),
                                               ::testing::Values("md-vv"),
                                               ::testing::Values("no", "v-rescale", "berendsen"),
-                                              ::testing::Values("no")),
+                                              ::testing::Values("no", "berendsen", "c-rescale")),
                            ::testing::Values("GMX_DISABLE_MODULAR_SIMULATOR")));
 INSTANTIATE_TEST_CASE_P(
         DISABLED_SimulatorsAreEquivalentDefaultLegacy,
         SimulatorComparisonTest,
         ::testing::Combine(
-                ::testing::Combine(::testing::Values("argon12", "tip3p5"),
-                                   ::testing::Values("md"),
-                                   ::testing::Values("no", "v-rescale", "berendsen", "nose-hoover"),
-                                   ::testing::Values("no", "Parrinello-Rahman")),
+                ::testing::Combine(
+                        ::testing::Values("argon12", "tip3p5"),
+                        ::testing::Values("md"),
+                        ::testing::Values("no", "v-rescale", "berendsen", "nose-hoover"),
+                        ::testing::Values("no", "Parrinello-Rahman", "berendsen", "c-rescale")),
                 ::testing::Values("GMX_USE_MODULAR_SIMULATOR")));
 #endif