Implement Berendsen thermostat for modular simulator

Refs #3423

diff --git a/src/gromacs/modularsimulator/modularsimulator.cpp b/src/gromacs/modularsimulator/modularsimulator.cpp
index 5e4aa552cd3f1401d3822368906bc756af64a7f9..bc49d4665cb1145941880b23a426d2c988f728c5 100644 (file)
--- a/src/gromacs/modularsimulator/modularsimulator.cpp
+++ b/src/gromacs/modularsimulator/modularsimulator.cpp
@@ -108,7 +108,8 @@ void ModularSimulator::addIntegrationElements(ModularSimulatorAlgorithmBuilder*
         // The leap frog integration algorithm
         builder->add<ForceElement>();
         builder->add<StatePropagatorData::Element>();
-        if (legacySimulatorData_->inputrec->etc == etcVRESCALE)
+        if (legacySimulatorData_->inputrec->etc == etcVRESCALE
+            || legacySimulatorData_->inputrec->etc == etcBERENDSEN)
         {
             builder->add<VelocityScalingTemperatureCoupling>(-1, UseFullStepKE::No,
                                                              ReportPreviousStepConservedEnergy::No);
@@ -140,7 +141,8 @@ void ModularSimulator::addIntegrationElements(ModularSimulatorAlgorithmBuilder*
         }
         builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>();
         builder->add<StatePropagatorData::Element>();
-        if (legacySimulatorData_->inputrec->etc == etcVRESCALE)
+        if (legacySimulatorData_->inputrec->etc == etcVRESCALE
+            || legacySimulatorData_->inputrec->etc == etcBERENDSEN)
         {
             builder->add<VelocityScalingTemperatureCoupling>(
                     0, UseFullStepKE::Yes, ReportPreviousStepConservedEnergy::Yes);
@@ -219,11 +221,11 @@ bool ModularSimulator::isInputCompatible(bool                             exitOn
     isInputCompatible =
             isInputCompatible
             && conditionalAssert(!doRerun, "Rerun is not supported by the modular simulator.");
-    isInputCompatible =
-            isInputCompatible
-            && conditionalAssert(
-                       inputrec->etc == etcNO || inputrec->etc == etcVRESCALE,
-                       "Only v-rescale thermostat is supported by the modular simulator.");
+    isInputCompatible = isInputCompatible
+                        && conditionalAssert(inputrec->etc == etcNO || inputrec->etc == etcVRESCALE
+                                                     || inputrec->etc == etcBERENDSEN,
+                                             "Only v-rescale and Berendsen thermostat are "
+                                             "supported by the modular simulator.");
     isInputCompatible =
             isInputCompatible
             && conditionalAssert(

diff --git a/src/gromacs/modularsimulator/velocityscalingtemperaturecoupling.cpp b/src/gromacs/modularsimulator/velocityscalingtemperaturecoupling.cpp
index 8ceaf49255bddcaf3a4e44623df76648cee6a3f6..e6633898272e50ee12cd518909387746b84be9db 100644 (file)
--- a/src/gromacs/modularsimulator/velocityscalingtemperaturecoupling.cpp
+++ b/src/gromacs/modularsimulator/velocityscalingtemperaturecoupling.cpp
@@ -199,6 +199,68 @@ private:
     ArrayRef<real> lambdaStartVelocities_;
 };
 
+/*! \internal
+ * \brief Implements Berendsen temperature coupling
+ */
+class BerendsenTemperatureCoupling final : public ITemperatureCouplingImpl
+{
+public:
+    //! Apply the v-rescale temperature control
+    real apply(Step gmx_unused                step,
+               int                            temperatureGroup,
+               real                           currentKineticEnergy,
+               real                           currentTemperature,
+               const TemperatureCouplingData& temperatureCouplingData) override
+    {
+        if (!(temperatureCouplingData.couplingTime[temperatureGroup] >= 0
+              && temperatureCouplingData.numDegreesOfFreedom[temperatureGroup] > 0
+              && currentKineticEnergy > 0))
+        {
+            lambdaStartVelocities_[temperatureGroup] = 1.0;
+            return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup];
+        }
+
+        real lambda =
+                std::sqrt(1.0
+                          + (temperatureCouplingData.couplingTimeStep
+                             / temperatureCouplingData.couplingTime[temperatureGroup])
+                                    * (temperatureCouplingData.referenceTemperature[temperatureGroup] / currentTemperature
+                                       - 1.0));
+        lambdaStartVelocities_[temperatureGroup] = std::max<real>(std::min<real>(lambda, 1.25), 0.8);
+        if (debug)
+        {
+            fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n", temperatureGroup,
+                    currentTemperature, lambdaStartVelocities_[temperatureGroup]);
+        }
+        return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup]
+               - (lambdaStartVelocities_[temperatureGroup] * lambdaStartVelocities_[temperatureGroup]
+                  - 1) * currentKineticEnergy;
+    }
+
+    //! Connect with propagator - Berendsen only scales start step velocities
+    void connectWithPropagator(const PropagatorThermostatConnection& connectionData,
+                               int                                   numTemperatureGroups) override
+    {
+        connectionData.setNumVelocityScalingVariables(numTemperatureGroups);
+        lambdaStartVelocities_ = connectionData.getViewOnVelocityScaling();
+    }
+
+    //! No data to write to checkpoint
+    void writeCheckpoint(std::optional<WriteCheckpointData> gmx_unused checkpointData,
+                         const t_commrec gmx_unused* cr) override
+    {
+    }
+    //! No data to read from checkpoints
+    void readCheckpoint(std::optional<ReadCheckpointData> gmx_unused checkpointData,
+                        const t_commrec gmx_unused* cr) override
+    {
+    }
+
+private:
+    //! View on the scaling factor of the propagator (pre-step velocities)
+    ArrayRef<real> lambdaStartVelocities_;
+};
+
 VelocityScalingTemperatureCoupling::VelocityScalingTemperatureCoupling(
         int                               nstcouple,
         int                               offset,
@@ -233,6 +295,10 @@ VelocityScalingTemperatureCoupling::VelocityScalingTemperatureCoupling(
     {
         temperatureCouplingImpl_ = std::make_unique<VRescaleTemperatureCoupling>(seed);
     }
+    else if (couplingType == etcBERENDSEN)
+    {
+        temperatureCouplingImpl_ = std::make_unique<BerendsenTemperatureCoupling>();
+    }
     else
     {
         throw NotImplementedError("Temperature coupling " + std::string(ETCOUPLTYPE(couplingType))

diff --git a/src/programs/mdrun/tests/exactcontinuation.cpp b/src/programs/mdrun/tests/exactcontinuation.cpp
index f3ff5f6279196825c11d38ab26eb71e1b3ea253f..737d98a38946e40ba25e8df3b7925cecc229d4f4 100644 (file)
--- a/src/programs/mdrun/tests/exactcontinuation.cpp
+++ b/src/programs/mdrun/tests/exactcontinuation.cpp
@@ -359,7 +359,8 @@ TEST_P(MdrunNoAppendContinuationIsExact, WithinTolerances)
     // TODO: Update this as modular simulator gains functionality
     const bool isModularSimulatorExplicitlyDisabled = (getenv("GMX_DISABLE_MODULAR_SIMULATOR") != nullptr);
     const bool isTCouplingCompatibleWithModularSimulator =
-            (temperatureCoupling == "no" || temperatureCoupling == "v-rescale");
+            (temperatureCoupling == "no" || temperatureCoupling == "v-rescale"
+             || temperatureCoupling == "berendsen");
     if (integrator == "md-vv" && pressureCoupling == "parrinello-rahman"
         && (isModularSimulatorExplicitlyDisabled || !isTCouplingCompatibleWithModularSimulator))
     {

diff --git a/src/programs/mdrun/tests/simulator.cpp b/src/programs/mdrun/tests/simulator.cpp
index 004b31e2afb6e1f92f0d015ab87ee2a273d29b70..6ed4f5f3c27531ced441ad871c7d9472b8df4a23 100644 (file)
--- a/src/programs/mdrun/tests/simulator.cpp
+++ b/src/programs/mdrun/tests/simulator.cpp
@@ -222,35 +222,37 @@ TEST_P(SimulatorComparisonTest, WithinTolerances)
 // These tests are very sensitive, so we only run them in double precision.
 // As we change call ordering, they might actually become too strict to be useful.
 #if !GMX_GPU_OPENCL && GMX_DOUBLE
-INSTANTIATE_TEST_CASE_P(SimulatorsAreEquivalentDefaultModular,
-                        SimulatorComparisonTest,
-                        ::testing::Combine(::testing::Combine(::testing::Values("argon12", "tip3p5"),
-                                                              ::testing::Values("md-vv"),
-                                                              ::testing::Values("no", "v-rescale"),
-                                                              ::testing::Values("no")),
-                                           ::testing::Values("GMX_DISABLE_MODULAR_SIMULATOR")));
+INSTANTIATE_TEST_CASE_P(
+        SimulatorsAreEquivalentDefaultModular,
+        SimulatorComparisonTest,
+        ::testing::Combine(::testing::Combine(::testing::Values("argon12", "tip3p5"),
+                                              ::testing::Values("md-vv"),
+                                              ::testing::Values("no", "v-rescale", "berendsen"),
+                                              ::testing::Values("no")),
+                           ::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"),
+                                              ::testing::Values("no", "v-rescale", "berendsen"),
                                               ::testing::Values("no", "Parrinello-Rahman")),
                            ::testing::Values("GMX_USE_MODULAR_SIMULATOR")));
 #else
-INSTANTIATE_TEST_CASE_P(DISABLED_SimulatorsAreEquivalentDefaultModular,
-                        SimulatorComparisonTest,
-                        ::testing::Combine(::testing::Combine(::testing::Values("argon12", "tip3p5"),
-                                                              ::testing::Values("md-vv"),
-                                                              ::testing::Values("no", "v-rescale"),
-                                                              ::testing::Values("no")),
-                                           ::testing::Values("GMX_DISABLE_MODULAR_SIMULATOR")));
+INSTANTIATE_TEST_CASE_P(
+        DISABLED_SimulatorsAreEquivalentDefaultModular,
+        SimulatorComparisonTest,
+        ::testing::Combine(::testing::Combine(::testing::Values("argon12", "tip3p5"),
+                                              ::testing::Values("md-vv"),
+                                              ::testing::Values("no", "v-rescale", "berendsen"),
+                                              ::testing::Values("no")),
+                           ::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"),
+                                              ::testing::Values("no", "v-rescale", "berendsen"),
                                               ::testing::Values("no", "Parrinello-Rahman")),
                            ::testing::Values("GMX_USE_MODULAR_SIMULATOR")));
 #endif