Refactor temperature annealing

This refactors temperature annealing by dividing it into two subfunctions,
one which returns a new temperature, keeping the inputrec constant and
being independent of the Update object, and one which updates the inputrec
and the Update object.

This separation allows most of the code to be reused by modular simulator,
and it is also a first step towards a design which does not need to modify
the input record (#3854)

Refs #3417, #3854

diff --git a/src/gromacs/mdlib/coupling.cpp b/src/gromacs/mdlib/coupling.cpp
index e8c2badf1f0a2bb11e118c328d622b731fe45f82..4c24f6350cc8ce0428ac4fc7139fa0f8c5d4ca50 100644 (file)
--- a/src/gromacs/mdlib/coupling.cpp
+++ b/src/gromacs/mdlib/coupling.cpp
@@ -2142,74 +2142,80 @@ bool initSimulatedAnnealing(t_inputrec* ir, gmx::Update* upd)
     return doSimAnnealing;
 }
 
-/* set target temperatures if we are annealing */
-void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd)
+real computeAnnealingTargetTemperature(const t_inputrec& inputrec, int temperatureGroup, real time)
 {
-    for (int temperatureGroup = 0; temperatureGroup < ir->opts.ngtc; temperatureGroup++)
+    GMX_RELEASE_ASSERT(temperatureGroup >= 0 && temperatureGroup < inputrec.opts.ngtc,
+                       "Invalid temperature group.");
+    if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::No)
     {
-        const auto& inputrec = *ir;
-        if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::No)
-        {
-            continue;
-        }
-        GMX_RELEASE_ASSERT(
-                inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Single
-                        || inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Periodic,
-                gmx::formatString("Unknown simulated annealing algorithm for temperature group %d", temperatureGroup)
-                        .c_str());
-        real       thist   = 0;
-        const auto npoints = inputrec.opts.anneal_npoints[temperatureGroup];
-        if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Periodic)
-        {
-            /* calculate time modulo the period */
-            const auto pert = inputrec.opts.anneal_time[temperatureGroup][npoints - 1];
-            const auto n    = static_cast<int>(t / pert);
-            thist           = t - n * pert; /* modulo time */
-            /* Make sure rounding didn't get us outside the interval */
-            if (std::fabs(thist - pert) < GMX_REAL_EPS * 100)
-            {
-                thist = 0;
-            }
-        }
-        else if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Single)
-        {
-            thist = t;
-        }
-        /* We are doing annealing for this group if we got here,
-         * and we have the (relative) time as thist.
-         * calculate target temp */
-        int j = 0;
-        while ((j < npoints - 1) && (thist > (inputrec.opts.anneal_time[temperatureGroup][j + 1])))
+        // No change of temperature, return current reference temperature
+        return inputrec.opts.ref_t[temperatureGroup];
+    }
+    GMX_RELEASE_ASSERT(
+            inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Single
+                    || inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Periodic,
+            gmx::formatString("Unknown simulated annealing algorithm for temperature group %d", temperatureGroup)
+                    .c_str());
+    real       thist   = 0;
+    const auto npoints = inputrec.opts.anneal_npoints[temperatureGroup];
+    if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Periodic)
+    {
+        /* calculate time modulo the period */
+        const auto pert = inputrec.opts.anneal_time[temperatureGroup][npoints - 1];
+        const auto n    = static_cast<int>(time / pert);
+        thist           = time - n * pert; /* modulo time */
+        /* Make sure rounding didn't get us outside the interval */
+        if (std::fabs(thist - pert) < GMX_REAL_EPS * 100)
         {
-            j++;
+            thist = 0;
         }
-        if (j < npoints - 1)
+    }
+    else if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Single)
+    {
+        thist = time;
+    }
+    /* We are doing annealing for this group if we got here,
+     * and we have the (relative) time as thist.
+     * calculate target temp */
+    int j = 0;
+    while ((j < npoints - 1) && (thist > (inputrec.opts.anneal_time[temperatureGroup][j + 1])))
+    {
+        j++;
+    }
+    if (j < npoints - 1)
+    {
+        /* Found our position between points j and j+1.
+         * Interpolate: x is the amount from j+1, (1-x) from point j
+         * First treat possible jumps in temperature as a special case.
+         */
+        if ((inputrec.opts.anneal_time[temperatureGroup][j + 1]
+             - inputrec.opts.anneal_time[temperatureGroup][j])
+            < GMX_REAL_EPS * 100)
         {
-            /* Found our position between points j and j+1.
-             * Interpolate: x is the amount from j+1, (1-x) from point j
-             * First treat possible jumps in temperature as a special case.
-             */
-            if ((inputrec.opts.anneal_time[temperatureGroup][j + 1]
-                 - inputrec.opts.anneal_time[temperatureGroup][j])
-                < GMX_REAL_EPS * 100)
-            {
-                ir->opts.ref_t[temperatureGroup] = inputrec.opts.anneal_temp[temperatureGroup][j + 1];
-            }
-            else
-            {
-                const real x = ((thist - inputrec.opts.anneal_time[temperatureGroup][j])
-                                / (inputrec.opts.anneal_time[temperatureGroup][j + 1]
-                                   - inputrec.opts.anneal_time[temperatureGroup][j]));
-                ir->opts.ref_t[temperatureGroup] =
-                        x * inputrec.opts.anneal_temp[temperatureGroup][j + 1]
-                        + (1 - x) * inputrec.opts.anneal_temp[temperatureGroup][j];
-            }
+            return inputrec.opts.anneal_temp[temperatureGroup][j + 1];
         }
         else
         {
-            ir->opts.ref_t[temperatureGroup] = inputrec.opts.anneal_temp[temperatureGroup][npoints - 1];
+            const real x = ((thist - inputrec.opts.anneal_time[temperatureGroup][j])
+                            / (inputrec.opts.anneal_time[temperatureGroup][j + 1]
+                               - inputrec.opts.anneal_time[temperatureGroup][j]));
+            return x * inputrec.opts.anneal_temp[temperatureGroup][j + 1]
+                   + (1 - x) * inputrec.opts.anneal_temp[temperatureGroup][j];
         }
     }
+    else
+    {
+        return inputrec.opts.anneal_temp[temperatureGroup][npoints - 1];
+    }
+}
+
+/* set target temperatures if we are annealing */
+void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd)
+{
+    for (int temperatureGroup = 0; temperatureGroup < ir->opts.ngtc; temperatureGroup++)
+    {
+        ir->opts.ref_t[temperatureGroup] = computeAnnealingTargetTemperature(*ir, temperatureGroup, t);
+    }
 
     upd->update_temperature_constants(*ir);
 }

diff --git a/src/gromacs/mdlib/coupling.h b/src/gromacs/mdlib/coupling.h
index 15211f5a83b9ed1227bd17f97fb4a65830560786..4e8857c3b0310cdd6f16637b12532e7622b49484 100644 (file)
--- a/src/gromacs/mdlib/coupling.h
+++ b/src/gromacs/mdlib/coupling.h
@@ -175,6 +175,16 @@ void rescale_velocities(const gmx_ekindata_t*               ekind,
                         gmx::ArrayRef<gmx::RVec>            v);
 /* Rescale the velocities with the scaling factor in ekind */
 
+/*!
+ * \brief Compute the new annealing temperature for a temperature group
+ *
+ * \param inputrec          The input record
+ * \param temperatureGroup  The temperature group
+ * \param time              The current time
+ * \return  The new reference temperature for the group
+ */
+real computeAnnealingTargetTemperature(const t_inputrec& inputrec, int temperatureGroup, real time);
+
 //! Check whether we do simulated annealing.
 bool doSimulatedAnnealing(const t_inputrec* ir);