Prepare simulated annealing refactoring

Prepare simulated annealing refactoring by moving variable declarations
to where they are first used, rename a non-descriptive integer `i`, use
a const reference to the inputrec where possible, replace a gmx_fatal
by an assert, and replace a switch / case construct by an if / else.

This is pure refactoring, no change in functionality, but will make
follow-up changes significantly easier to review.

diff --git a/src/gromacs/mdlib/coupling.cpp b/src/gromacs/mdlib/coupling.cpp
index d0f06c484845c2a968ebdbf4fb38cdce52c23ceb..e8c2badf1f0a2bb11e118c328d622b731fe45f82 100644 (file)
--- a/src/gromacs/mdlib/coupling.cpp
+++ b/src/gromacs/mdlib/coupling.cpp
@@ -2145,39 +2145,41 @@ bool initSimulatedAnnealing(t_inputrec* ir, gmx::Update* upd)
 /* set target temperatures if we are annealing */
 void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd)
 {
-    int  i, j, n, npoints;
-    real pert, thist = 0, x;
-
-    for (i = 0; i < ir->opts.ngtc; i++)
-    {
-        npoints = ir->opts.anneal_npoints[i];
-        switch (ir->opts.annealing[i])
-        {
-            case SimulatedAnnealing::No: continue;
-            case SimulatedAnnealing::Periodic:
-                /* calculate time modulo the period */
-                pert  = ir->opts.anneal_time[i][npoints - 1];
-                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;
-                }
-                break;
-            case SimulatedAnnealing::Single: thist = t; break;
-            default:
-                gmx_fatal(FARGS,
-                          "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",
-                          i,
-                          ir->opts.ngtc,
-                          npoints);
+    for (int temperatureGroup = 0; temperatureGroup < ir->opts.ngtc; temperatureGroup++)
+    {
+        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 */
-        j = 0;
-        while ((j < npoints - 1) && (thist > (ir->opts.anneal_time[i][j + 1])))
+        int j = 0;
+        while ((j < npoints - 1) && (thist > (inputrec.opts.anneal_time[temperatureGroup][j + 1])))
         {
             j++;
         }
@@ -2187,21 +2189,25 @@ void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd)
              * Interpolate: x is the amount from j+1, (1-x) from point j
              * First treat possible jumps in temperature as a special case.
              */
-            if ((ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]) < GMX_REAL_EPS * 100)
+            if ((inputrec.opts.anneal_time[temperatureGroup][j + 1]
+                 - inputrec.opts.anneal_time[temperatureGroup][j])
+                < GMX_REAL_EPS * 100)
             {
-                ir->opts.ref_t[i] = ir->opts.anneal_temp[i][j + 1];
+                ir->opts.ref_t[temperatureGroup] = inputrec.opts.anneal_temp[temperatureGroup][j + 1];
             }
             else
             {
-                x = ((thist - ir->opts.anneal_time[i][j])
-                     / (ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]));
-                ir->opts.ref_t[i] =
-                        x * ir->opts.anneal_temp[i][j + 1] + (1 - x) * ir->opts.anneal_temp[i][j];
+                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];
             }
         }
         else
         {
-            ir->opts.ref_t[i] = ir->opts.anneal_temp[i][npoints - 1];
+            ir->opts.ref_t[temperatureGroup] = inputrec.opts.anneal_temp[temperatureGroup][npoints - 1];
         }
     }