Reimplement constant acceleration groups
[alexxy/gromacs.git] / src / gromacs / mdtypes / inputrec.cpp
index e42cfe2385f763b346edb0ec5ccc58855b969b16..d93e399c4f3ce703839a108807dddca21ec29baa 100644 (file)
@@ -3,7 +3,8 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2010, The GROMACS development team.
- * Copyright (c) 2012,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2012,2014,2015,2016,2017 by the GROMACS development team.
+ * Copyright (c) 2018,2019,2020,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.
@@ -36,6 +37,7 @@
  */
 #include "gmxpre.h"
 
+#include "gromacs/utility/enumerationhelpers.h"
 #include "inputrec.h"
 
 #include <cstdio>
 #include <cstring>
 
 #include <algorithm>
+#include <memory>
+#include <numeric>
 
+#include "gromacs/applied_forces/awh/read_params.h"
 #include "gromacs/math/veccompare.h"
 #include "gromacs/math/vecdump.h"
-#include "gromacs/mdtypes/awh-params.h"
+#include "gromacs/mdlib/vcm.h"
+#include "gromacs/mdtypes/awh_params.h"
 #include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/mdtypes/pull-params.h"
+#include "gromacs/mdtypes/multipletimestepping.h"
+#include "gromacs/mdtypes/pull_params.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/utility/compare.h"
 #include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/snprintf.h"
 #include "gromacs/utility/strconvert.h"
+#include "gromacs/utility/stringutil.h"
 #include "gromacs/utility/textwriter.h"
 #include "gromacs/utility/txtdump.h"
 
 //! Macro to select a bool name
-#define EBOOL(e)       gmx::boolToString(e)
+#define EBOOL(e) gmx::boolToString(e)
+
+/* Default values for nstcalcenergy, used when the are no other restrictions. */
+constexpr int c_defaultNstCalcEnergy = 10;
 
 /* The minimum number of integration steps required for reasonably accurate
  * integration of first and second order coupling algorithms.
  */
-const int nstmin_berendsen_tcouple =  5;
+const int nstmin_berendsen_tcouple = 5;
 const int nstmin_berendsen_pcouple = 10;
 const int nstmin_harmonic          = 20;
 
-t_inputrec::t_inputrec()
+/* Default values for T- and P- coupling intervals, used when the are no other
+ * restrictions.
+ */
+constexpr int c_defaultNstTCouple = 10;
+constexpr int c_defaultNstPCouple = 10;
+
+t_inputrec::t_inputrec() :
+    fepvals(std::make_unique<t_lambda>()),
+    simtempvals(std::make_unique<t_simtemp>()),
+    expandedvals(std::make_unique<t_expanded>())
 {
-    std::memset(this, 0, sizeof(*this));
-    snew(fepvals, 1);
-    snew(expandedvals, 1);
-    snew(simtempvals, 1);
 }
 
 t_inputrec::~t_inputrec()
@@ -83,56 +99,50 @@ t_inputrec::~t_inputrec()
     done_inputrec(this);
 }
 
-static int nst_wanted(const t_inputrec *ir)
+int ir_optimal_nstcalcenergy(const t_inputrec* ir)
 {
+    int nst;
+
     if (ir->nstlist > 0)
     {
-        return ir->nstlist;
+        nst = ir->nstlist;
     }
     else
     {
-        return 10;
+        nst = c_defaultNstCalcEnergy;
     }
-}
 
-int ir_optimal_nstcalcenergy(const t_inputrec *ir)
-{
-    return nst_wanted(ir);
+    if (ir->useMts)
+    {
+        nst = std::lcm(nst, ir->mtsLevels.back().stepFactor);
+    }
+
+    return nst;
 }
 
-int tcouple_min_integration_steps(int etc)
+int tcouple_min_integration_steps(TemperatureCoupling etc)
 {
     int n;
 
     switch (etc)
     {
-        case etcNO:
-            n = 0;
-            break;
-        case etcBERENDSEN:
-        case etcYES:
-            n = nstmin_berendsen_tcouple;
-            break;
-        case etcVRESCALE:
+        case TemperatureCoupling::No: n = 0; break;
+        case TemperatureCoupling::Berendsen:
+        case TemperatureCoupling::Yes: n = nstmin_berendsen_tcouple; break;
+        case TemperatureCoupling::VRescale:
             /* V-rescale supports instantaneous rescaling */
             n = 0;
             break;
-        case etcNOSEHOOVER:
-            n = nstmin_harmonic;
-            break;
-        case etcANDERSEN:
-        case etcANDERSENMASSIVE:
-            n = 1;
-            break;
-        default:
-            gmx_incons("Unknown etc value");
-            n = 0;
+        case TemperatureCoupling::NoseHoover: n = nstmin_harmonic; break;
+        case TemperatureCoupling::Andersen:
+        case TemperatureCoupling::AndersenMassive: n = 1; break;
+        default: gmx_incons("Unknown etc value");
     }
 
     return n;
 }
 
-int ir_optimal_nsttcouple(const t_inputrec *ir)
+int ir_optimal_nsttcouple(const t_inputrecir)
 {
     int  nmin, nwanted, n;
     real tau_min;
@@ -140,10 +150,10 @@ int ir_optimal_nsttcouple(const t_inputrec *ir)
 
     nmin = tcouple_min_integration_steps(ir->etc);
 
-    nwanted = nst_wanted(ir);
+    nwanted = c_defaultNstTCouple;
 
     tau_min = 1e20;
-    if (ir->etc != etcNO)
+    if (ir->etc != TemperatureCoupling::No)
     {
         for (g = 0; g < ir->opts.ngtc; g++)
         {
@@ -154,13 +164,13 @@ int ir_optimal_nsttcouple(const t_inputrec *ir)
         }
     }
 
-    if (nmin == 0 || ir->delta_t*nwanted <= tau_min)
+    if (nmin == 0 || ir->delta_t * nwanted <= tau_min)
     {
         n = nwanted;
     }
     else
     {
-        n = (int)(tau_min/(ir->delta_t*nmin) + 0.001);
+        n = static_cast<int>(tau_min / (ir->delta_t * nmin) + 0.001);
         if (n < 1)
         {
             n = 1;
@@ -174,193 +184,162 @@ int ir_optimal_nsttcouple(const t_inputrec *ir)
     return n;
 }
 
-int pcouple_min_integration_steps(int epc)
+int pcouple_min_integration_steps(PressureCoupling epc)
 {
     int n;
 
     switch (epc)
     {
-        case epcNO:
-            n = 0;
-            break;
-        case etcBERENDSEN:
-        case epcISOTROPIC:
-            n = nstmin_berendsen_pcouple;
-            break;
-        case epcPARRINELLORAHMAN:
-        case epcMTTK:
-            n = nstmin_harmonic;
-            break;
-        default:
-            gmx_incons("Unknown epc value");
-            n = 0;
+        case PressureCoupling::No: n = 0; break;
+        case PressureCoupling::Berendsen:
+        case PressureCoupling::CRescale:
+        case PressureCoupling::Isotropic: n = nstmin_berendsen_pcouple; break;
+        case PressureCoupling::ParrinelloRahman:
+        case PressureCoupling::Mttk: n = nstmin_harmonic; break;
+        default: gmx_incons("Unknown epc value");
     }
 
     return n;
 }
 
-int ir_optimal_nstpcouple(const t_inputrec *ir)
+int ir_optimal_nstpcouple(const t_inputrecir)
 {
-    int  nmin, nwanted, n;
+    const int minIntegrationSteps = pcouple_min_integration_steps(ir->epc);
 
-    nmin = pcouple_min_integration_steps(ir->epc);
+    const int nwanted = c_defaultNstPCouple;
 
-    nwanted = nst_wanted(ir);
+    // With multiple time-stepping we can only compute the pressure at slowest steps
+    const int minNstPCouple = (ir->useMts ? ir->mtsLevels.back().stepFactor : 1);
 
-    if (nmin == 0 || ir->delta_t*nwanted <= ir->tau_p)
+    int n;
+    if (minIntegrationSteps == 0 || ir->delta_t * nwanted <= ir->tau_p)
     {
         n = nwanted;
     }
     else
     {
-        n = static_cast<int>(ir->tau_p/(ir->delta_t*nmin) + 0.001);
-        if (n < 1)
+        n = static_cast<int>(ir->tau_p / (ir->delta_t * minIntegrationSteps) + 0.001);
+        if (n < minNstPCouple)
         {
-            n = 1;
+            n = minNstPCouple;
         }
-        while (nwanted % n != 0)
+        // Without MTS we try to make nstpcouple a "nice" number
+        if (!ir->useMts)
         {
-            n--;
+            while (nwanted % n != 0)
+            {
+                n--;
+            }
         }
     }
 
+    // With MTS, nstpcouple should be a multiple of the slowest MTS interval
+    if (ir->useMts)
+    {
+        n = n - (n % minNstPCouple);
+    }
+
     return n;
 }
 
-gmx_bool ir_coulomb_switched(const t_inputrec *ir)
+gmx_bool ir_coulomb_switched(const t_inputrecir)
 {
-    return (ir->coulombtype == eelSWITCH ||
-            ir->coulombtype == eelSHIFT ||
-            ir->coulombtype == eelENCADSHIFT ||
-            ir->coulombtype == eelPMESWITCH ||
-            ir->coulombtype == eelPMEUSERSWITCH ||
-            ir->coulomb_modifier == eintmodPOTSWITCH ||
-            ir->coulomb_modifier == eintmodFORCESWITCH);
+    return (ir->coulombtype == CoulombInteractionType::Switch
+            || ir->coulombtype == CoulombInteractionType::Shift
+            || ir->coulombtype == CoulombInteractionType::PmeSwitch
+            || ir->coulombtype == CoulombInteractionType::PmeUserSwitch
+            || ir->coulomb_modifier == InteractionModifiers::PotSwitch
+            || ir->coulomb_modifier == InteractionModifiers::ForceSwitch);
 }
 
-gmx_bool ir_coulomb_is_zero_at_cutoff(const t_inputrec *ir)
+gmx_bool ir_coulomb_is_zero_at_cutoff(const t_inputrecir)
 {
-    return (ir->cutoff_scheme == ecutsVERLET ||
-            ir_coulomb_switched(ir) || ir->coulomb_modifier != eintmodNONE ||
-            ir->coulombtype == eelRF_ZERO);
+    return (ir->cutoff_scheme == CutoffScheme::Verlet || ir_coulomb_switched(ir)
+            || ir->coulomb_modifier != InteractionModifiers::None
+            || ir->coulombtype == CoulombInteractionType::RFZero);
 }
 
-gmx_bool ir_coulomb_might_be_zero_at_cutoff(const t_inputrec *ir)
+gmx_bool ir_coulomb_might_be_zero_at_cutoff(const t_inputrecir)
 {
-    return (ir_coulomb_is_zero_at_cutoff(ir) || ir->coulombtype == eelUSER || ir->coulombtype == eelPMEUSER);
+    return (ir_coulomb_is_zero_at_cutoff(ir) || ir->coulombtype == CoulombInteractionType::User
+            || ir->coulombtype == CoulombInteractionType::PmeUser);
 }
 
-gmx_bool ir_vdw_switched(const t_inputrec *ir)
+gmx_bool ir_vdw_switched(const t_inputrecir)
 {
-    return (ir->vdwtype == evdwSWITCH ||
-            ir->vdwtype == evdwSHIFT ||
-            ir->vdwtype == evdwENCADSHIFT ||
-            ir->vdw_modifier == eintmodPOTSWITCH ||
-            ir->vdw_modifier == eintmodFORCESWITCH);
+    return (ir->vdwtype == VanDerWaalsType::Switch || ir->vdwtype == VanDerWaalsType::Shift
+            || ir->vdw_modifier == InteractionModifiers::PotSwitch
+            || ir->vdw_modifier == InteractionModifiers::ForceSwitch);
 }
 
-gmx_bool ir_vdw_is_zero_at_cutoff(const t_inputrec *ir)
+gmx_bool ir_vdw_is_zero_at_cutoff(const t_inputrecir)
 {
-    return (ir->cutoff_scheme == ecutsVERLET ||
-            ir_vdw_switched(ir) || ir->vdw_modifier != eintmodNONE);
+    return (ir->cutoff_scheme == CutoffScheme::Verlet || ir_vdw_switched(ir)
+            || ir->vdw_modifier != InteractionModifiers::None);
 }
 
-gmx_bool ir_vdw_might_be_zero_at_cutoff(const t_inputrec *ir)
+gmx_bool ir_vdw_might_be_zero_at_cutoff(const t_inputrecir)
 {
-    return (ir_vdw_is_zero_at_cutoff(ir) || ir->vdwtype == evdwUSER);
+    return (ir_vdw_is_zero_at_cutoff(ir) || ir->vdwtype == VanDerWaalsType::User);
 }
 
-static void done_pull_group(t_pull_group *pgrp)
+static void done_t_rot(t_rot* rot)
 {
-    if (pgrp->nat > 0)
+    if (rot == nullptr)
     {
-        sfree(pgrp->ind);
-        sfree(pgrp->weight);
+        return;
     }
-}
-
-static void done_pull_params(pull_params_t *pull)
-{
-    int i;
-
-    for (i = 0; i < pull->ngroup+1; i++)
+    if (rot->grp != nullptr)
     {
-        done_pull_group(pull->group);
+        for (int i = 0; i < rot->ngrp; i++)
+        {
+            sfree(rot->grp[i].ind);
+            sfree(rot->grp[i].x_ref);
+        }
+        sfree(rot->grp);
     }
-
-    sfree(pull->group);
-    sfree(pull->coord);
+    sfree(rot);
 }
 
-static void done_lambdas(t_lambda *fep)
+static void done_t_swapCoords(t_swapcoords* swapCoords)
 {
-    if (fep->n_lambda > 0)
+    if (swapCoords == nullptr)
     {
-        for (int i = 0; i < efptNR; i++)
-        {
-            sfree(fep->all_lambda[i]);
-        }
+        return;
+    }
+    for (int i = 0; i < swapCoords->ngrp; i++)
+    {
+        sfree(swapCoords->grp[i].ind);
+        sfree(swapCoords->grp[i].molname);
     }
-    sfree(fep->all_lambda);
+    sfree(swapCoords->grp);
+    sfree(swapCoords);
 }
 
-void done_inputrec(t_inputrec *ir)
+void done_inputrec(t_inputrecir)
 {
     sfree(ir->opts.nrdf);
     sfree(ir->opts.ref_t);
+    for (int i = 0; i < ir->opts.ngtc; i++)
+    {
+        sfree(ir->opts.anneal_time[i]);
+        sfree(ir->opts.anneal_temp[i]);
+    }
     sfree(ir->opts.annealing);
     sfree(ir->opts.anneal_npoints);
     sfree(ir->opts.anneal_time);
     sfree(ir->opts.anneal_temp);
     sfree(ir->opts.tau_t);
-    sfree(ir->opts.acc);
+    sfree(ir->opts.acceleration);
     sfree(ir->opts.nFreeze);
-    sfree(ir->opts.QMmethod);
-    sfree(ir->opts.QMbasis);
-    sfree(ir->opts.QMcharge);
-    sfree(ir->opts.QMmult);
-    sfree(ir->opts.bSH);
-    sfree(ir->opts.CASorbitals);
-    sfree(ir->opts.CASelectrons);
-    sfree(ir->opts.SAon);
-    sfree(ir->opts.SAoff);
-    sfree(ir->opts.SAsteps);
     sfree(ir->opts.egp_flags);
-    done_lambdas(ir->fepvals);
-    sfree(ir->fepvals);
-    sfree(ir->expandedvals);
-    sfree(ir->simtempvals);
 
-    if (ir->pull)
-    {
-        done_pull_params(ir->pull);
-        sfree(ir->pull);
-    }
+    done_t_swapCoords(ir->swap);
+    done_t_rot(ir->rot);
     delete ir->params;
 }
 
-static void pr_qm_opts(FILE *fp, int indent, const char *title, const t_grpopts *opts)
-{
-    fprintf(fp, "%s:\n", title);
-
-    pr_int(fp, indent, "ngQM", opts->ngQM);
-    if (opts->ngQM > 0)
-    {
-        pr_ivec(fp, indent, "QMmethod", opts->QMmethod, opts->ngQM, FALSE);
-        pr_ivec(fp, indent, "QMbasis", opts->QMbasis, opts->ngQM, FALSE);
-        pr_ivec(fp, indent, "QMcharge", opts->QMcharge, opts->ngQM, FALSE);
-        pr_ivec(fp, indent, "QMmult", opts->QMmult, opts->ngQM, FALSE);
-        pr_bvec(fp, indent, "SH", opts->bSH, opts->ngQM, FALSE);
-        pr_ivec(fp, indent, "CASorbitals", opts->CASorbitals, opts->ngQM, FALSE);
-        pr_ivec(fp, indent, "CASelectrons", opts->CASelectrons, opts->ngQM, FALSE);
-        pr_rvec(fp, indent, "SAon", opts->SAon, opts->ngQM, FALSE);
-        pr_rvec(fp, indent, "SAoff", opts->SAoff, opts->ngQM, FALSE);
-        pr_ivec(fp, indent, "SAsteps", opts->SAsteps, opts->ngQM, FALSE);
-    }
-}
-
-static void pr_grp_opts(FILE *out, int indent, const char *title, const t_grpopts *opts,
-                        gmx_bool bMDPformat)
+static void pr_grp_opts(FILE* out, int indent, const char* title, const t_grpopts* opts, gmx_bool bMDPformat)
 {
     int i, m, j;
 
@@ -397,7 +376,7 @@ static void pr_grp_opts(FILE *out, int indent, const char *title, const t_grpopt
     fprintf(out, "annealing%s", bMDPformat ? " = " : ":");
     for (i = 0; (i < opts->ngtc); i++)
     {
-        fprintf(out, "  %10s", EANNEAL(opts->annealing[i]));
+        fprintf(out, "  %10s", enumValueToString(opts->annealing[i]));
     }
     fprintf(out, "\n");
 
@@ -433,7 +412,7 @@ static void pr_grp_opts(FILE *out, int indent, const char *title, const t_grpopt
     {
         for (m = 0; (m < DIM); m++)
         {
-            fprintf(out, "  %10g", opts->acc[i][m]);
+            fprintf(out, "  %10g", opts->acceleration[i][m]);
         }
     }
     fprintf(out, "\n");
@@ -456,7 +435,7 @@ static void pr_grp_opts(FILE *out, int indent, const char *title, const t_grpopt
         fprintf(out, "energygrp-flags[%3d]:", i);
         for (m = 0; (m < opts->ngener); m++)
         {
-            fprintf(out, " %d", opts->egp_flags[opts->ngener*i+m]);
+            fprintf(out, " %d", opts->egp_flags[opts->ngener * i + m]);
         }
         fprintf(out, "\n");
     }
@@ -464,13 +443,19 @@ static void pr_grp_opts(FILE *out, int indent, const char *title, const t_grpopt
     fflush(out);
 }
 
-static void pr_matrix(FILE *fp, int indent, const char *title, const rvec *m,
-                      gmx_bool bMDPformat)
+static void pr_matrix(FILE* fp, int indent, const char* title, const rvec* m, gmx_bool bMDPformat)
 {
     if (bMDPformat)
     {
-        fprintf(fp, "%-10s    = %g %g %g %g %g %g\n", title,
-                m[XX][XX], m[YY][YY], m[ZZ][ZZ], m[XX][YY], m[XX][ZZ], m[YY][ZZ]);
+        fprintf(fp,
+                "%-10s    = %g %g %g %g %g %g\n",
+                title,
+                m[XX][XX],
+                m[YY][YY],
+                m[ZZ][ZZ],
+                m[XX][YY],
+                m[XX][ZZ],
+                m[YY][ZZ]);
     }
     else
     {
@@ -484,34 +469,30 @@ static void pr_matrix(FILE *fp, int indent, const char *title, const rvec *m,
 #define PR(t, s) pr_real(fp, indent, t, s)
 #define PD(t, s) pr_double(fp, indent, t, s)
 
-static void pr_pull_group(FILE *fp, int indent, int g, const t_pull_group *pgrp)
+static void pr_pull_group(FILE* fp, int indent, int g, const t_pull_group* pgrp)
 {
     pr_indent(fp, indent);
     fprintf(fp, "pull-group %d:\n", g);
     indent += 2;
-    pr_ivec_block(fp, indent, "atom", pgrp->ind, pgrp->nat, TRUE);
-    pr_rvec(fp, indent, "weight", pgrp->weight, pgrp->nweight, TRUE);
+    pr_ivec_block(fp, indent, "atom", pgrp->ind.data(), pgrp->ind.size(), TRUE);
+    pr_rvec(fp, indent, "weight", pgrp->weight.data(), pgrp->weight.size(), TRUE);
     PI("pbcatom", pgrp->pbcatom);
 }
 
-static void pr_pull_coord(FILE *fp, int indent, int c, const t_pull_coord *pcrd)
+static void pr_pull_coord(FILE* fp, int indent, int c, const t_pull_coord* pcrd)
 {
-    int g;
-
     pr_indent(fp, indent);
     fprintf(fp, "pull-coord %d:\n", c);
-    PS("type", EPULLTYPE(pcrd->eType));
-    if (pcrd->eType == epullEXTERNAL)
+    PS("type", enumValueToString(pcrd->eType));
+    if (pcrd->eType == PullingAlgorithm::External)
     {
-        PS("potential-provider", pcrd->externalPotentialProvider);
+        PS("potential-provider", pcrd->externalPotentialProvider.c_str());
     }
-    PS("geometry", EPULLGEOM(pcrd->eGeom));
-    for (g = 0; g < pcrd->ngroup; g++)
+    PS("geometry", enumValueToString(pcrd->eGeom));
+    for (int g = 0; g < pcrd->ngroup; g++)
     {
-        char buf[10];
-
-        sprintf(buf, "group[%d]", g);
-        PI(buf, pcrd->group[g]);
+        std::string buffer = gmx::formatString("group[%d]", g);
+        PI(buffer.c_str(), pcrd->group[g]);
     }
     pr_ivec(fp, indent, "dim", pcrd->dim, DIM, TRUE);
     pr_rvec(fp, indent, "origin", pcrd->origin, DIM, TRUE);
@@ -523,38 +504,38 @@ static void pr_pull_coord(FILE *fp, int indent, int c, const t_pull_coord *pcrd)
     PR("kB", pcrd->kB);
 }
 
-static void pr_simtempvals(FILE *fp, int indent, const t_simtemp *simtemp, int n_lambda)
+static void pr_simtempvals(FILE* fp, int indent, const t_simtemp* simtemp, int n_lambda)
 {
-    PS("simulated-tempering-scaling", ESIMTEMP(simtemp->eSimTempScale));
+    PS("simulated-tempering-scaling", enumValueToString(simtemp->eSimTempScale));
     PR("sim-temp-low", simtemp->simtemp_low);
     PR("sim-temp-high", simtemp->simtemp_high);
-    pr_rvec(fp, indent, "simulated tempering temperatures", simtemp->temperatures, n_lambda, TRUE);
+    pr_rvec(fp, indent, "simulated tempering temperatures", simtemp->temperatures.data(), n_lambda, TRUE);
 }
 
-static void pr_expandedvals(FILE *fp, int indent, const t_expanded *expand, int n_lambda)
+static void pr_expandedvals(FILE* fp, int indent, const t_expanded* expand, int n_lambda)
 {
 
     PI("nstexpanded", expand->nstexpanded);
-    PS("lmc-stats", elamstats_names[expand->elamstats]);
-    PS("lmc-move", elmcmove_names[expand->elmcmove]);
-    PS("lmc-weights-equil", elmceq_names[expand->elmceq]);
-    if (expand->elmceq == elmceqNUMATLAM)
+    PS("lmc-stats", enumValueToString(expand->elamstats));
+    PS("lmc-move", enumValueToString(expand->elmcmove));
+    PS("lmc-weights-equil", enumValueToString(expand->elmceq));
+    if (expand->elmceq == LambdaWeightWillReachEquilibrium::NumAtLambda)
     {
         PI("weight-equil-number-all-lambda", expand->equil_n_at_lam);
     }
-    if (expand->elmceq == elmceqSAMPLES)
+    if (expand->elmceq == LambdaWeightWillReachEquilibrium::Samples)
     {
         PI("weight-equil-number-samples", expand->equil_samples);
     }
-    if (expand->elmceq == elmceqSTEPS)
+    if (expand->elmceq == LambdaWeightWillReachEquilibrium::Steps)
     {
         PI("weight-equil-number-steps", expand->equil_steps);
     }
-    if (expand->elmceq == elmceqWLDELTA)
+    if (expand->elmceq == LambdaWeightWillReachEquilibrium::WLDelta)
     {
         PR("weight-equil-wl-delta", expand->equil_wl_delta);
     }
-    if (expand->elmceq == elmceqRATIO)
+    if (expand->elmceq == LambdaWeightWillReachEquilibrium::Ratio)
     {
         PR("weight-equil-count-ratio", expand->equil_ratio);
     }
@@ -573,13 +554,13 @@ static void pr_expandedvals(FILE *fp, int indent, const t_expanded *expand, int
     PS("wl-oneovert", EBOOL(expand->bWLoneovert));
 
     pr_indent(fp, indent);
-    pr_rvec(fp, indent, "init-lambda-weights", expand->init_lambda_weights, n_lambda, TRUE);
+    pr_rvec(fp, indent, "init-lambda-weights", expand->init_lambda_weights.data(), n_lambda, TRUE);
     PS("init-weights", EBOOL(expand->bInit_weights));
 }
 
-static void pr_fepvals(FILE *fp, int indent, const t_lambda *fep, gmx_bool bMDPformat)
+static void pr_fepvals(FILE* fp, int indent, const t_lambda* fep, gmx_bool bMDPformat)
 {
-    int i, j;
+    int j;
 
     PR("init-lambda", fep->init_lambda);
     PI("init-lambda-state", fep->init_fep_state);
@@ -594,9 +575,9 @@ static void pr_fepvals(FILE *fp, int indent, const t_lambda *fep, gmx_bool bMDPf
     {
         pr_indent(fp, indent);
         fprintf(fp, "separate-dvdl%s\n", bMDPformat ? " = " : ":");
-        for (i = 0; i < efptNR; i++)
+        for (auto i : gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, bool>::keys())
         {
-            fprintf(fp, "%18s = ", efpt_names[i]);
+            fprintf(fp, "%18s = ", enumValueToString(i));
             if (fep->separate_dvdl[i])
             {
                 fprintf(fp, "  TRUE");
@@ -608,9 +589,10 @@ static void pr_fepvals(FILE *fp, int indent, const t_lambda *fep, gmx_bool bMDPf
             fprintf(fp, "\n");
         }
         fprintf(fp, "all-lambdas%s\n", bMDPformat ? " = " : ":");
-        for (i = 0; i < efptNR; i++)
+        for (auto key : gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, bool>::keys())
         {
-            fprintf(fp, "%18s = ", efpt_names[i]);
+            fprintf(fp, "%18s = ", enumValueToString(key));
+            int i = static_cast<int>(key);
             for (j = 0; j < fep->n_lambda; j++)
             {
                 fprintf(fp, "  %10g", fep->all_lambda[i][j]);
@@ -619,7 +601,7 @@ static void pr_fepvals(FILE *fp, int indent, const t_lambda *fep, gmx_bool bMDPf
         }
     }
     PI("calc-lambda-neighbors", fep->lambda_neighbors);
-    PS("dhdl-print-energy", edHdLPrintEnergy_names[fep->edHdLPrintEnergy]);
+    PS("dhdl-print-energy", enumValueToString(fep->edHdLPrintEnergy));
     PR("sc-alpha", fep->sc_alpha);
     PI("sc-power", fep->sc_power);
     PR("sc-r-power", fep->sc_r_power);
@@ -628,132 +610,129 @@ static void pr_fepvals(FILE *fp, int indent, const t_lambda *fep, gmx_bool bMDPf
     PS("sc-coul", EBOOL(fep->bScCoul));
     PI("dh-hist-size", fep->dh_hist_size);
     PD("dh-hist-spacing", fep->dh_hist_spacing);
-    PS("separate-dhdl-file", SEPDHDLFILETYPE(fep->separate_dhdl_file));
-    PS("dhdl-derivatives", DHDLDERIVATIVESTYPE(fep->dhdl_derivatives));
+    PS("separate-dhdl-file", enumValueToString(fep->separate_dhdl_file));
+    PS("dhdl-derivatives", enumValueToString(fep->dhdl_derivatives));
+    PS("sc-function", enumValueToString(fep->softcoreFunction));
+    PR("sc-gapsys-scale-linpoint-lj", fep->scGapsysScaleLinpointLJ);
+    PR("sc-gapsys-scale-linpoint-q", fep->scGapsysScaleLinpointQ);
+    PR("sc-gapsys-sigma-lj", fep->scGapsysSigmaLJ);
 };
 
-static void pr_pull(FILE *fp, int indent, const pull_params_t *pull)
+static void pr_pull(FILE* fp, int indent, const pull_params_t& pull)
 {
     int g;
 
-    PR("pull-cylinder-r", pull->cylinder_r);
-    PR("pull-constr-tol", pull->constr_tol);
-    PS("pull-print-COM", EBOOL(pull->bPrintCOM));
-    PS("pull-print-ref-value", EBOOL(pull->bPrintRefValue));
-    PS("pull-print-components", EBOOL(pull->bPrintComp));
-    PI("pull-nstxout", pull->nstxout);
-    PI("pull-nstfout", pull->nstfout);
-    PI("pull-ngroups", pull->ngroup);
-    for (g = 0; g < pull->ngroup; g++)
+    PR("pull-cylinder-r", pull.cylinder_r);
+    PR("pull-constr-tol", pull.constr_tol);
+    PS("pull-print-COM", EBOOL(pull.bPrintCOM));
+    PS("pull-print-ref-value", EBOOL(pull.bPrintRefValue));
+    PS("pull-print-components", EBOOL(pull.bPrintComp));
+    PI("pull-nstxout", pull.nstxout);
+    PI("pull-nstfout", pull.nstfout);
+    PS("pull-pbc-ref-prev-step-com", EBOOL(pull.bSetPbcRefToPrevStepCOM));
+    PS("pull-xout-average", EBOOL(pull.bXOutAverage));
+    PS("pull-fout-average", EBOOL(pull.bFOutAverage));
+    PI("pull-ngroups", pull.ngroup);
+    for (g = 0; g < pull.ngroup; g++)
     {
-        pr_pull_group(fp, indent, g, &pull->group[g]);
+        pr_pull_group(fp, indent, g, &pull.group[g]);
     }
-    PI("pull-ncoords", pull->ncoord);
-    for (g = 0; g < pull->ncoord; g++)
+    PI("pull-ncoords", pull.ncoord);
+    for (g = 0; g < pull.ncoord; g++)
     {
-        pr_pull_coord(fp, indent, g, &pull->coord[g]);
+        pr_pull_coord(fp, indent, g, &pull.coord[g]);
     }
 }
 
-static void pr_awh_bias_dim(FILE *fp, int indent, gmx::AwhDimParams *awhDimParams, char *prefix)
+static void pr_awh_bias_dim(FILE* fp, int indent, const gmx::AwhDimParams& awhDimParams, const char* prefix)
 {
     pr_indent(fp, indent);
     indent++;
-    fprintf(fp,  "%s:\n", prefix);
-    PS("coord-provider", EAWHCOORDPROVIDER(awhDimParams->eCoordProvider));
-    PI("coord-index", awhDimParams->coordIndex + 1);
-    PR("start", awhDimParams->origin);
-    PR("end", awhDimParams->end);
-    PR("period", awhDimParams->period);
-    PR("force-constant", awhDimParams->forceConstant);
-    PR("diffusion", awhDimParams->diffusion);
-    PR("start", awhDimParams->origin);
-    PR("end", awhDimParams->end);
-    PR("cover-diameter", awhDimParams->coverDiameter);
+    fprintf(fp, "%s:\n", prefix);
+    PS("coord-provider", enumValueToString(awhDimParams.coordinateProvider()));
+    PI("coord-index", awhDimParams.coordinateIndex() + 1);
+    PR("start", awhDimParams.origin());
+    PR("end", awhDimParams.end());
+    PR("period", awhDimParams.period());
+    PR("force-constant", awhDimParams.forceConstant());
+    PR("diffusion", awhDimParams.diffusion());
+    PR("cover-diameter", awhDimParams.coverDiameter());
 }
 
-static void pr_awh_bias(FILE *fp, int indent, gmx::AwhBiasParams *awhBiasParams, char *prefix)
+static void pr_awh_bias(FILE* fp, int indent, const gmx::AwhBiasParams& awhBiasParams, const char* prefix)
 {
     char opt[STRLEN];
 
     sprintf(opt, "%s-error-init", prefix);
-    PR(opt, awhBiasParams->errorInitial);
+    PR(opt, awhBiasParams.initialErrorEstimate());
     sprintf(opt, "%s-growth", prefix);
-    PS(opt, EAWHGROWTH(awhBiasParams->eGrowth));
+    PS(opt, enumValueToString(awhBiasParams.growthType()));
     sprintf(opt, "%s-target", prefix);
-    PS(opt, EAWHTARGET(awhBiasParams->eTarget));
+    PS(opt, enumValueToString(awhBiasParams.targetDistribution()));
     sprintf(opt, "%s-target-beta-scalng", prefix);
-    PR(opt, awhBiasParams->targetBetaScaling);
+    PR(opt, awhBiasParams.targetBetaScaling());
     sprintf(opt, "%s-target-cutoff", prefix);
-    PR(opt, awhBiasParams->targetCutoff);
+    PR(opt, awhBiasParams.targetCutoff());
     sprintf(opt, "%s-user-data", prefix);
-    PS(opt, EBOOL(awhBiasParams->bUserData));
+    PS(opt, EBOOL(awhBiasParams.userPMFEstimate()));
     sprintf(opt, "%s-share-group", prefix);
-    PI(opt, awhBiasParams->shareGroup);
+    PI(opt, awhBiasParams.shareGroup());
     sprintf(opt, "%s-equilibrate-histogram", prefix);
-    PS(opt, EBOOL(awhBiasParams->equilibrateHistogram));
+    PS(opt, EBOOL(awhBiasParams.equilibrateHistogram()));
     sprintf(opt, "%s-ndim", prefix);
-    PI(opt, awhBiasParams->ndim);
+    PI(opt, awhBiasParams.ndim());
 
-    for (int d = 0; d < awhBiasParams->ndim; d++)
+    int d = 0;
+    for (const auto& dimParam : awhBiasParams.dimParams())
     {
         char prefixdim[STRLEN];
         sprintf(prefixdim, "%s-dim%d", prefix, d + 1);
-        pr_awh_bias_dim(fp, indent, &awhBiasParams->dimParams[d], prefixdim);
+        pr_awh_bias_dim(fp, indent, dimParam, prefixdim);
+        d++;
     }
 }
 
-static void pr_awh(FILE *fp, int indent, gmx::AwhParams *awhParams)
+static void pr_awh(FILE* fp, int indent, gmx::AwhParams* awhParams)
 {
-    int  k;
-    char opt[STRLEN], prefix[STRLEN];
-
-    sprintf(prefix, "%s", "awh");
-
-    sprintf(opt, "%s-potential", prefix);
-    PS(opt, EAWHPOTENTIAL(awhParams->ePotential));
-    sprintf(opt, "%s-seed", prefix);
-    PI(opt, awhParams->seed);
-    sprintf(opt, "%s-nstout", prefix);
-    PI(opt, awhParams->nstOut);
-    sprintf(opt, "%s-nstsample", prefix);
-    PI(opt, awhParams->nstSampleCoord);
-    sprintf(opt, "%s-nsamples-update", prefix);
-    PI(opt, awhParams->numSamplesUpdateFreeEnergy);
-    sprintf(opt, "%s-share-bias-multisim", prefix);
-    PS(opt, EBOOL(awhParams->shareBiasMultisim));
-    sprintf(opt, "%s-nbias", prefix);
-    PI(opt, awhParams->numBias);
-
-    for (k = 0; k < awhParams->numBias; k++)
-    {
-        sprintf(prefix, "awh%d", k + 1);
-        pr_awh_bias(fp, indent, &awhParams->awhBiasParams[k], prefix);
+    PS("awh-potential", enumValueToString(awhParams->potential()));
+    PI("awh-seed", awhParams->seed());
+    PI("awh-nstout", awhParams->nstout());
+    PI("awh-nstsample", awhParams->nstSampleCoord());
+    PI("awh-nsamples-update", awhParams->numSamplesUpdateFreeEnergy());
+    PS("awh-share-bias-multisim", EBOOL(awhParams->shareBiasMultisim()));
+    PI("awh-nbias", awhParams->numBias());
+
+    int k = 0;
+    for (const auto& awhBiasParam : awhParams->awhBiasParams())
+    {
+        auto prefix = gmx::formatString("awh%d", k + 1);
+        pr_awh_bias(fp, indent, awhBiasParam, prefix.c_str());
+        k++;
     }
 }
 
-static void pr_rotgrp(FILE *fp, int indent, int g, const t_rotgrp *rotg)
+static void pr_rotgrp(FILE* fp, int indent, int g, const t_rotgrp* rotg)
 {
     pr_indent(fp, indent);
     fprintf(fp, "rot-group %d:\n", g);
     indent += 2;
-    PS("rot-type", EROTGEOM(rotg->eType));
+    PS("rot-type", enumValueToString(rotg->eType));
     PS("rot-massw", EBOOL(rotg->bMassW));
     pr_ivec_block(fp, indent, "atom", rotg->ind, rotg->nat, TRUE);
     pr_rvecs(fp, indent, "x-ref", rotg->x_ref, rotg->nat);
-    pr_rvec(fp, indent, "rot-vec", rotg->vec, DIM, TRUE);
+    pr_rvec(fp, indent, "rot-vec", rotg->inputVec, DIM, TRUE);
     pr_rvec(fp, indent, "rot-pivot", rotg->pivot, DIM, TRUE);
     PR("rot-rate", rotg->rate);
     PR("rot-k", rotg->k);
     PR("rot-slab-dist", rotg->slab_dist);
     PR("rot-min-gauss", rotg->min_gaussian);
     PR("rot-eps", rotg->eps);
-    PS("rot-fit-method", EROTFIT(rotg->eFittype));
+    PS("rot-fit-method", enumValueToString(rotg->eFittype));
     PI("rot-potfit-nstep", rotg->PotAngle_nstep);
     PR("rot-potfit-step", rotg->PotAngle_step);
 }
 
-static void pr_rot(FILE *fp, int indent, const t_rot *rot)
+static void pr_rot(FILE* fp, int indent, const t_rot* rot)
 {
     int g;
 
@@ -767,13 +746,15 @@ static void pr_rot(FILE *fp, int indent, const t_rot *rot)
 }
 
 
-static void pr_swap(FILE *fp, int indent, const t_swapcoords *swap)
+static void pr_swap(FILE* fp, int indent, const t_swapcoords* swap)
 {
     char str[STRLEN];
 
     /* Enums for better readability of the code */
-    enum {
-        eCompA = 0, eCompB
+    enum
+    {
+        eCompA = 0,
+        eCompB
     };
 
 
@@ -789,11 +770,19 @@ static void pr_swap(FILE *fp, int indent, const t_swapcoords *swap)
     }
 
     /* The solvent group */
-    snprintf(str, STRLEN, "solvent group %s", swap->grp[eGrpSolvent].molname);
-    pr_ivec_block(fp, indent, str, swap->grp[eGrpSolvent].ind, swap->grp[eGrpSolvent].nat, TRUE);
+    snprintf(str,
+             STRLEN,
+             "solvent group %s",
+             swap->grp[static_cast<int>(SwapGroupSplittingType::Solvent)].molname);
+    pr_ivec_block(fp,
+                  indent,
+                  str,
+                  swap->grp[static_cast<int>(SwapGroupSplittingType::Solvent)].ind,
+                  swap->grp[static_cast<int>(SwapGroupSplittingType::Solvent)].nat,
+                  TRUE);
 
     /* Now print the indices for all the ion groups: */
-    for (int ig = eSwapFixedGrpNR; ig < swap->ngrp; ig++)
+    for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < swap->ngrp; ig++)
     {
         snprintf(str, STRLEN, "ion group %s", swap->grp[ig].molname);
         pr_ivec_block(fp, indent, str, swap->grp[ig].ind, swap->grp[ig].nat, TRUE);
@@ -810,9 +799,9 @@ static void pr_swap(FILE *fp, int indent, const t_swapcoords *swap)
     /* Print the requested ion counts for both compartments */
     for (int ic = eCompA; ic <= eCompB; ic++)
     {
-        for (int ig = eSwapFixedGrpNR; ig < swap->ngrp; ig++)
+        for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < swap->ngrp; ig++)
         {
-            snprintf(str, STRLEN, "%s-in-%c", swap->grp[ig].molname, 'A'+ic);
+            snprintf(str, STRLEN, "%s-in-%c", swap->grp[ig].molname, 'A' + ic);
             PI(str, swap->grp[ig].nmolReq[ic]);
         }
     }
@@ -823,17 +812,16 @@ static void pr_swap(FILE *fp, int indent, const t_swapcoords *swap)
 }
 
 
-static void pr_imd(FILE *fp, int indent, const t_IMD *imd)
+static void pr_imd(FILE* fp, int indent, const t_IMD* imd)
 {
     PI("IMD-atoms", imd->nat);
     pr_ivec_block(fp, indent, "atom", imd->ind, imd->nat, TRUE);
 }
 
 
-void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
-                 gmx_bool bMDPformat)
+void pr_inputrec(FILE* fp, int indent, const char* title, const t_inputrec* ir, gmx_bool bMDPformat)
 {
-    const char *infbuf = "inf";
+    const charinfbuf = "inf";
 
     if (available(fp, ir, indent, title))
     {
@@ -845,13 +833,37 @@ void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
          * options are written in the default mdout.mdp, and with
          * the same user-exposed names to facilitate debugging.
          */
-        PS("integrator", EI(ir->eI));
+        PS("integrator", enumValueToString(ir->eI));
         PR("tinit", ir->init_t);
         PR("dt", ir->delta_t);
         PSTEP("nsteps", ir->nsteps);
         PSTEP("init-step", ir->init_step);
         PI("simulation-part", ir->simulation_part);
-        PS("comm-mode", ECOM(ir->comm_mode));
+        PS("mts", EBOOL(ir->useMts));
+        if (ir->useMts)
+        {
+            for (int mtsIndex = 1; mtsIndex < static_cast<int>(ir->mtsLevels.size()); mtsIndex++)
+            {
+                const auto&       mtsLevel = ir->mtsLevels[mtsIndex];
+                const std::string forceKey = gmx::formatString("mts-level%d-forces", mtsIndex + 1);
+                std::string       forceGroups;
+                for (int i = 0; i < static_cast<int>(gmx::MtsForceGroups::Count); i++)
+                {
+                    if (mtsLevel.forceGroups[i])
+                    {
+                        if (!forceGroups.empty())
+                        {
+                            forceGroups += " ";
+                        }
+                        forceGroups += gmx::mtsForceGroupNames[i];
+                    }
+                }
+                PS(forceKey.c_str(), forceGroups.c_str());
+                const std::string factorKey = gmx::formatString("mts-level%d-factor", mtsIndex + 1);
+                PI(factorKey.c_str(), mtsLevel.stepFactor);
+            }
+        }
+        PS("comm-mode", enumValueToString(ir->comm_mode));
         PI("nstcomm", ir->nstcomm);
 
         /* Langevin dynamics */
@@ -880,17 +892,16 @@ void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
         PR("compressed-x-precision", ir->x_compression_precision);
 
         /* Neighborsearching parameters */
-        PS("cutoff-scheme", ECUTSCHEME(ir->cutoff_scheme));
+        PS("cutoff-scheme", enumValueToString(ir->cutoff_scheme));
         PI("nstlist", ir->nstlist);
-        PS("ns-type", ENS(ir->ns_type));
-        PS("pbc", epbc_names[ir->ePBC]);
+        PS("pbc", c_pbcTypeNames[ir->pbcType].c_str());
         PS("periodic-molecules", EBOOL(ir->bPeriodicMols));
         PR("verlet-buffer-tolerance", ir->verletbuf_tol);
         PR("rlist", ir->rlist);
 
         /* Options for electrostatics and VdW */
-        PS("coulombtype", EELTYPE(ir->coulombtype));
-        PS("coulomb-modifier", INTMODIFIER(ir->coulomb_modifier));
+        PS("coulombtype", enumValueToString(ir->coulombtype));
+        PS("coulomb-modifier", enumValueToString(ir->coulomb_modifier));
         PR("rcoulomb-switch", ir->rcoulomb_switch);
         PR("rcoulomb", ir->rcoulomb);
         if (ir->epsilon_r != 0)
@@ -909,11 +920,11 @@ void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
         {
             PS("epsilon-rf", infbuf);
         }
-        PS("vdw-type", EVDWTYPE(ir->vdwtype));
-        PS("vdw-modifier", INTMODIFIER(ir->vdw_modifier));
+        PS("vdw-type", enumValueToString(ir->vdwtype));
+        PS("vdw-modifier", enumValueToString(ir->vdw_modifier));
         PR("rvdw-switch", ir->rvdw_switch);
         PR("rvdw", ir->rvdw);
-        PS("DispCorr", EDISPCORR(ir->eDispCorr));
+        PS("DispCorr", enumValueToString(ir->eDispCorr));
         PR("table-extension", ir->tabext);
 
         PR("fourierspacing", ir->fourier_spacing);
@@ -923,30 +934,28 @@ void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
         PI("pme-order", ir->pme_order);
         PR("ewald-rtol", ir->ewald_rtol);
         PR("ewald-rtol-lj", ir->ewald_rtol_lj);
-        PS("lj-pme-comb-rule", ELJPMECOMBNAMES(ir->ljpme_combination_rule));
-        PR("ewald-geometry", ir->ewald_geometry);
+        PS("lj-pme-comb-rule", enumValueToString(ir->ljpme_combination_rule));
+        PS("ewald-geometry", enumValueToString(ir->ewald_geometry));
         PR("epsilon-surface", ir->epsilon_surface);
 
         /* Options for weak coupling algorithms */
-        PS("tcoupl", ETCOUPLTYPE(ir->etc));
+        PS("tcoupl", enumValueToString(ir->etc));
         PI("nsttcouple", ir->nsttcouple);
         PI("nh-chain-length", ir->opts.nhchainlength);
         PS("print-nose-hoover-chain-variables", EBOOL(ir->bPrintNHChains));
 
-        PS("pcoupl", EPCOUPLTYPE(ir->epc));
-        PS("pcoupltype", EPCOUPLTYPETYPE(ir->epct));
+        PS("pcoupl", enumValueToString(ir->epc));
+        PS("pcoupltype", enumValueToString(ir->epct));
         PI("nstpcouple", ir->nstpcouple);
         PR("tau-p", ir->tau_p);
         pr_matrix(fp, indent, "compressibility", ir->compress, bMDPformat);
         pr_matrix(fp, indent, "ref-p", ir->ref_p, bMDPformat);
-        PS("refcoord-scaling", EREFSCALINGTYPE(ir->refcoord_scaling));
+        PS("refcoord-scaling", enumValueToString(ir->refcoord_scaling));
 
         if (bMDPformat)
         {
-            fprintf(fp, "posres-com  = %g %g %g\n", ir->posres_com[XX],
-                    ir->posres_com[YY], ir->posres_com[ZZ]);
-            fprintf(fp, "posres-comB = %g %g %g\n", ir->posres_comB[XX],
-                    ir->posres_comB[YY], ir->posres_comB[ZZ]);
+            fprintf(fp, "posres-com  = %g %g %g\n", ir->posres_com[XX], ir->posres_com[YY], ir->posres_com[ZZ]);
+            fprintf(fp, "posres-comB = %g %g %g\n", ir->posres_comB[XX], ir->posres_comB[YY], ir->posres_comB[ZZ]);
         }
         else
         {
@@ -956,13 +965,11 @@ void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
 
         /* QMMM */
         PS("QMMM", EBOOL(ir->bQMMM));
-        PI("QMconstraints", ir->QMconstraints);
-        PI("QMMMscheme", ir->QMMMscheme);
-        PR("MMChargeScaleFactor", ir->scalefactor);
-        pr_qm_opts(fp, indent, "qm-opts", &(ir->opts));
+        fprintf(fp, "%s:\n", "qm-opts");
+        pr_int(fp, indent, "ngQM", ir->opts.ngQM);
 
         /* CONSTRAINT OPTIONS */
-        PS("constraint-algorithm", ECONSTRTYPE(ir->eConstrAlg));
+        PS("constraint-algorithm", enumValueToString(ir->eConstrAlg));
         PS("continuation", EBOOL(ir->bContinuation));
 
         PS("Shake-SOR", EBOOL(ir->bShakeSOR));
@@ -973,7 +980,7 @@ void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
 
         /* Walls */
         PI("nwall", ir->nwall);
-        PS("wall-type", EWALLTYPE(ir->wall_type));
+        PS("wall-type", enumValueToString(ir->wall_type));
         PR("wall-r-linpot", ir->wall_r_linpot);
         /* wall-atomtype */
         PI("wall-atomtype[0]", ir->wall_atomtype[0]);
@@ -987,14 +994,14 @@ void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
         PS("pull", EBOOL(ir->bPull));
         if (ir->bPull)
         {
-            pr_pull(fp, indent, ir->pull);
+            pr_pull(fp, indent, *ir->pull);
         }
 
         /* AWH BIASING */
         PS("awh", EBOOL(ir->bDoAwh));
         if (ir->bDoAwh)
         {
-            pr_awh(fp, indent, ir->awhParams);
+            pr_awh(fp, indent, ir->awhParams.get());
         }
 
         /* ENFORCED ROTATION */
@@ -1012,8 +1019,8 @@ void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
         }
 
         /* NMR refinement stuff */
-        PS("disre", EDISRETYPE(ir->eDisre));
-        PS("disre-weighting", EDISREWEIGHTING(ir->eDisreWeighting));
+        PS("disre", enumValueToString(ir->eDisre));
+        PS("disre-weighting", enumValueToString(ir->eDisreWeighting));
         PS("disre-mixed", EBOOL(ir->bDisreMixed));
         PR("dr-fc", ir->dr_fc);
         PR("dr-tau", ir->dr_tau);
@@ -1024,14 +1031,14 @@ void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
         PR("nstorireout", ir->nstorireout);
 
         /* FREE ENERGY VARIABLES */
-        PS("free-energy", EFEPTYPE(ir->efep));
-        if (ir->efep != efepNO || ir->bSimTemp)
+        PS("free-energy", enumValueToString(ir->efep));
+        if (ir->efep != FreeEnergyPerturbationType::No || ir->bSimTemp)
         {
-            pr_fepvals(fp, indent, ir->fepvals, bMDPformat);
+            pr_fepvals(fp, indent, ir->fepvals.get(), bMDPformat);
         }
         if (ir->bExpanded)
         {
-            pr_expandedvals(fp, indent, ir->expandedvals, ir->fepvals->n_lambda);
+            pr_expandedvals(fp, indent, ir->expandedvals.get(), ir->fepvals->n_lambda);
         }
 
         /* NON-equilibrium MD stuff */
@@ -1042,12 +1049,12 @@ void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
         PS("simulated-tempering", EBOOL(ir->bSimTemp));
         if (ir->bSimTemp)
         {
-            pr_simtempvals(fp, indent, ir->simtempvals, ir->fepvals->n_lambda);
+            pr_simtempvals(fp, indent, ir->simtempvals.get(), ir->fepvals->n_lambda);
         }
 
         /* ION/WATER SWAPPING FOR COMPUTATIONAL ELECTROPHYSIOLOGY */
-        PS("swapcoords", ESWAPTYPE(ir->eSwapCoords));
-        if (ir->eSwapCoords != eswapNO)
+        PS("swapcoords", enumValueToString(ir->eSwapCoords));
+        if (ir->eSwapCoords != SwapType::No)
         {
             pr_swap(fp, indent, ir->swap);
         }
@@ -1076,12 +1083,12 @@ void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
 #undef PR
 #undef PI
 
-static void cmp_grpopts(FILE *fp, const t_grpopts *opt1, const t_grpopts *opt2, real ftol, real abstol)
+static void cmp_grpopts(FILE* fp, const t_grpopts* opt1, const t_grpopts* opt2, real ftol, real abstol)
 {
     int  i, j;
     char buf1[256], buf2[256];
 
-    cmp_int(fp, "inputrec->grpopts.ngtc", -1,  opt1->ngtc, opt2->ngtc);
+    cmp_int(fp, "inputrec->grpopts.ngtc", -1, opt1->ngtc, opt2->ngtc);
     cmp_int(fp, "inputrec->grpopts.ngacc", -1, opt1->ngacc, opt2->ngacc);
     cmp_int(fp, "inputrec->grpopts.ngfrz", -1, opt1->ngfrz, opt2->ngfrz);
     cmp_int(fp, "inputrec->grpopts.ngener", -1, opt1->ngener, opt2->ngener);
@@ -1090,9 +1097,8 @@ static void cmp_grpopts(FILE *fp, const t_grpopts *opt1, const t_grpopts *opt2,
         cmp_real(fp, "inputrec->grpopts.nrdf", i, opt1->nrdf[i], opt2->nrdf[i], ftol, abstol);
         cmp_real(fp, "inputrec->grpopts.ref_t", i, opt1->ref_t[i], opt2->ref_t[i], ftol, abstol);
         cmp_real(fp, "inputrec->grpopts.tau_t", i, opt1->tau_t[i], opt2->tau_t[i], ftol, abstol);
-        cmp_int(fp, "inputrec->grpopts.annealing", i, opt1->annealing[i], opt2->annealing[i]);
-        cmp_int(fp, "inputrec->grpopts.anneal_npoints", i,
-                opt1->anneal_npoints[i], opt2->anneal_npoints[i]);
+        cmpEnum(fp, "inputrec->grpopts.annealing", opt1->annealing[i], opt2->annealing[i]);
+        cmp_int(fp, "inputrec->grpopts.anneal_npoints", i, opt1->anneal_npoints[i], opt2->anneal_npoints[i]);
         if (opt1->anneal_npoints[i] == opt2->anneal_npoints[i])
         {
             sprintf(buf1, "inputrec->grpopts.anneal_time[%d]", i);
@@ -1111,15 +1117,13 @@ static void cmp_grpopts(FILE *fp, const t_grpopts *opt1, const t_grpopts *opt2,
             for (j = i; j < opt1->ngener; j++)
             {
                 sprintf(buf1, "inputrec->grpopts.egp_flags[%d]", i);
-                cmp_int(fp, buf1, j,
-                        opt1->egp_flags[opt1->ngener*i+j],
-                        opt2->egp_flags[opt1->ngener*i+j]);
+                cmp_int(fp, buf1, j, opt1->egp_flags[opt1->ngener * i + j], opt2->egp_flags[opt1->ngener * i + j]);
             }
         }
     }
     for (i = 0; (i < std::min(opt1->ngacc, opt2->ngacc)); i++)
     {
-        cmp_rvec(fp, "inputrec->grpopts.acc", i, opt1->acc[i], opt2->acc[i], ftol, abstol);
+        cmp_rvec(fp, "inputrec->grpopts.acceleration", i, opt1->acceleration[i], opt2->acceleration[i], ftol, abstol);
     }
     for (i = 0; (i < std::min(opt1->ngfrz, opt2->ngfrz)); i++)
     {
@@ -1127,74 +1131,153 @@ static void cmp_grpopts(FILE *fp, const t_grpopts *opt1, const t_grpopts *opt2,
     }
 }
 
-static void cmp_pull(FILE *fp)
+static void cmp_pull(FILEfp)
 {
-    fprintf(fp, "WARNING: Both files use COM pulling, but comparing of the pull struct is not implemented (yet). The pull parameters could be the same or different.\n");
+    fprintf(fp,
+            "WARNING: Both files use COM pulling, but comparing of the pull struct is not "
+            "implemented (yet). The pull parameters could be the same or different.\n");
 }
 
-static void cmp_awhDimParams(FILE *fp, const gmx::AwhDimParams *dimp1, const gmx::AwhDimParams *dimp2, int dimIndex, real ftol, real abstol)
+static void cmp_awhDimParams(FILE*                    fp,
+                             const gmx::AwhDimParams& dimp1,
+                             const gmx::AwhDimParams& dimp2,
+                             int                      dimIndex,
+                             real                     ftol,
+                             real                     abstol)
 {
     /* Note that we have double index here, but the compare functions only
      * support one index, so here we only print the dim index and not the bias.
      */
-    cmp_int(fp, "inputrec.awhParams->bias?->dim->coord_index", dimIndex, dimp1->coordIndex, dimp2->coordIndex);
-    cmp_double(fp, "inputrec->awhParams->bias?->dim->period", dimIndex, dimp1->period, dimp2->period, ftol, abstol);
-    cmp_double(fp, "inputrec->awhParams->bias?->dim->diffusion", dimIndex, dimp1->diffusion, dimp2->diffusion, ftol, abstol);
-    cmp_double(fp, "inputrec->awhParams->bias?->dim->origin", dimIndex, dimp1->origin, dimp2->origin, ftol, abstol);
-    cmp_double(fp, "inputrec->awhParams->bias?->dim->end", dimIndex, dimp1->end, dimp2->end, ftol, abstol);
-    cmp_double(fp, "inputrec->awhParams->bias?->dim->coord_value_init", dimIndex, dimp1->coordValueInit, dimp2->coordValueInit, ftol, abstol);
-    cmp_double(fp, "inputrec->awhParams->bias?->dim->coverDiameter", dimIndex, dimp1->coverDiameter, dimp2->coverDiameter, ftol, abstol);
+    cmp_int(fp,
+            "inputrec.awhParams->bias?->dim->coord_index",
+            dimIndex,
+            dimp1.coordinateIndex(),
+            dimp2.coordinateIndex());
+    cmp_double(fp, "inputrec->awhParams->bias?->dim->period", dimIndex, dimp1.period(), dimp2.period(), ftol, abstol);
+    cmp_double(fp,
+               "inputrec->awhParams->bias?->dim->diffusion",
+               dimIndex,
+               dimp1.diffusion(),
+               dimp2.diffusion(),
+               ftol,
+               abstol);
+    cmp_double(fp, "inputrec->awhParams->bias?->dim->origin", dimIndex, dimp1.origin(), dimp2.origin(), ftol, abstol);
+    cmp_double(fp, "inputrec->awhParams->bias?->dim->end", dimIndex, dimp1.end(), dimp2.end(), ftol, abstol);
+    cmp_double(fp,
+               "inputrec->awhParams->bias?->dim->coord_value_init",
+               dimIndex,
+               dimp1.initialCoordinate(),
+               dimp2.initialCoordinate(),
+               ftol,
+               abstol);
+    cmp_double(fp,
+               "inputrec->awhParams->bias?->dim->coverDiameter",
+               dimIndex,
+               dimp1.coverDiameter(),
+               dimp2.coverDiameter(),
+               ftol,
+               abstol);
 }
 
-static void cmp_awhBiasParams(FILE *fp, const gmx::AwhBiasParams *bias1, const gmx::AwhBiasParams *bias2, int biasIndex, real ftol, real abstol)
+static void cmp_awhBiasParams(FILE*                     fp,
+                              const gmx::AwhBiasParams& bias1,
+                              const gmx::AwhBiasParams& bias2,
+                              int                       biasIndex,
+                              real                      ftol,
+                              real                      abstol)
 {
-    cmp_int(fp, "inputrec->awhParams->ndim", biasIndex, bias1->ndim, bias2->ndim);
-    cmp_int(fp, "inputrec->awhParams->biaseTarget", biasIndex, bias1->eTarget, bias2->eTarget);
-    cmp_double(fp, "inputrec->awhParams->biastargetBetaScaling", biasIndex, bias1->targetBetaScaling, bias2->targetBetaScaling, ftol, abstol);
-    cmp_double(fp, "inputrec->awhParams->biastargetCutoff", biasIndex, bias1->targetCutoff, bias2->targetCutoff, ftol, abstol);
-    cmp_int(fp, "inputrec->awhParams->biaseGrowth", biasIndex, bias1->eGrowth, bias2->eGrowth);
-    cmp_bool(fp, "inputrec->awhParams->biasbUserData", biasIndex, bias1->bUserData, bias2->bUserData);
-    cmp_double(fp, "inputrec->awhParams->biaserror_initial", biasIndex, bias1->errorInitial, bias2->errorInitial, ftol, abstol);
-    cmp_int(fp, "inputrec->awhParams->biasShareGroup", biasIndex, bias1->shareGroup, bias2->shareGroup);
-
-    for (int dim = 0; dim < std::min(bias1->ndim, bias2->ndim); dim++)
-    {
-        cmp_awhDimParams(fp, &bias1->dimParams[dim], &bias2->dimParams[dim], dim, ftol, abstol);
+    cmp_int(fp, "inputrec->awhParams->ndim", biasIndex, bias1.ndim(), bias2.ndim());
+    cmpEnum<gmx::AwhTargetType>(
+            fp, "inputrec->awhParams->biaseTarget", bias1.targetDistribution(), bias2.targetDistribution());
+    cmp_double(fp,
+               "inputrec->awhParams->biastargetBetaScaling",
+               biasIndex,
+               bias1.targetBetaScaling(),
+               bias2.targetBetaScaling(),
+               ftol,
+               abstol);
+    cmp_double(fp,
+               "inputrec->awhParams->biastargetCutoff",
+               biasIndex,
+               bias1.targetCutoff(),
+               bias2.targetCutoff(),
+               ftol,
+               abstol);
+    cmpEnum<gmx::AwhHistogramGrowthType>(
+            fp, "inputrec->awhParams->biaseGrowth", bias1.growthType(), bias2.growthType());
+    cmp_bool(fp, "inputrec->awhParams->biasbUserData", biasIndex, bias1.userPMFEstimate(), bias2.userPMFEstimate());
+    cmp_double(fp,
+               "inputrec->awhParams->biaserror_initial",
+               biasIndex,
+               bias1.initialErrorEstimate(),
+               bias2.initialErrorEstimate(),
+               ftol,
+               abstol);
+    cmp_int(fp, "inputrec->awhParams->biasShareGroup", biasIndex, bias1.shareGroup(), bias2.shareGroup());
+
+    const auto dimParams1 = bias1.dimParams();
+    const auto dimParams2 = bias2.dimParams();
+    for (int dim = 0; dim < std::min(bias1.ndim(), bias2.ndim()); dim++)
+    {
+        cmp_awhDimParams(fp, dimParams1[dim], dimParams2[dim], dim, ftol, abstol);
     }
 }
 
-static void cmp_awhParams(FILE *fp, const gmx::AwhParams *awh1, const gmx::AwhParams *awh2, real ftol, real abstol)
+static void cmp_awhParams(FILE* fp, const gmx::AwhParams& awh1, const gmx::AwhParams& awh2, real ftol, real abstol)
 {
-    cmp_int(fp, "inputrec->awhParams->nbias", -1, awh1->numBias, awh2->numBias);
-    cmp_int64(fp, "inputrec->awhParams->seed", awh1->seed, awh2->seed);
-    cmp_int(fp, "inputrec->awhParams->nstout", -1, awh1->nstOut, awh2->nstOut);
-    cmp_int(fp, "inputrec->awhParams->nstsample_coord", -1, awh1->nstSampleCoord, awh2->nstSampleCoord);
-    cmp_int(fp, "inputrec->awhParams->nsamples_update_free_energy", -1, awh1->numSamplesUpdateFreeEnergy, awh2->numSamplesUpdateFreeEnergy);
-    cmp_int(fp, "inputrec->awhParams->ePotential", -1, awh1->ePotential, awh2->ePotential);
-    cmp_bool(fp, "inputrec->awhParams->shareBiasMultisim", -1, awh1->shareBiasMultisim, awh2->shareBiasMultisim);
-
-    if (awh1->numBias == awh2->numBias)
-    {
-        for (int bias = 0; bias < awh1->numBias; bias++)
+    cmp_int(fp, "inputrec->awhParams->nbias", -1, awh1.numBias(), awh2.numBias());
+    cmp_int64(fp, "inputrec->awhParams->seed", awh1.seed(), awh2.seed());
+    cmp_int(fp, "inputrec->awhParams->nstout", -1, awh1.nstout(), awh2.nstout());
+    cmp_int(fp, "inputrec->awhParams->nstsample_coord", -1, awh1.nstSampleCoord(), awh2.nstSampleCoord());
+    cmp_int(fp,
+            "inputrec->awhParams->nsamples_update_free_energy",
+            -1,
+            awh1.numSamplesUpdateFreeEnergy(),
+            awh2.numSamplesUpdateFreeEnergy());
+    cmpEnum<gmx::AwhPotentialType>(
+            fp, "inputrec->awhParams->ePotential", awh1.potential(), awh2.potential());
+    cmp_bool(fp, "inputrec->awhParams->shareBiasMultisim", -1, awh1.shareBiasMultisim(), awh2.shareBiasMultisim());
+
+    if (awh1.numBias() == awh2.numBias())
+    {
+        const auto awhBiasParams1 = awh1.awhBiasParams();
+        const auto awhBiasParams2 = awh2.awhBiasParams();
+        for (int bias = 0; bias < awh1.numBias(); bias++)
         {
-            cmp_awhBiasParams(fp, &awh1->awhBiasParams[bias], &awh2->awhBiasParams[bias], bias, ftol, abstol);
+            cmp_awhBiasParams(fp, awhBiasParams1[bias], awhBiasParams2[bias], bias, ftol, abstol);
         }
     }
 }
 
-static void cmp_simtempvals(FILE *fp, const t_simtemp *simtemp1, const t_simtemp *simtemp2, int n_lambda, real ftol, real abstol)
+static void cmp_simtempvals(FILE*            fp,
+                            const t_simtemp* simtemp1,
+                            const t_simtemp* simtemp2,
+                            int              n_lambda,
+                            real             ftol,
+                            real             abstol)
 {
     int i;
-    cmp_int(fp, "inputrec->simtempvals->eSimTempScale", -1, simtemp1->eSimTempScale, simtemp2->eSimTempScale);
+    cmpEnum(fp, "inputrec->simtempvals->eSimTempScale", simtemp1->eSimTempScale, simtemp2->eSimTempScale);
     cmp_real(fp, "inputrec->simtempvals->simtemp_high", -1, simtemp1->simtemp_high, simtemp2->simtemp_high, ftol, abstol);
     cmp_real(fp, "inputrec->simtempvals->simtemp_low", -1, simtemp1->simtemp_low, simtemp2->simtemp_low, ftol, abstol);
     for (i = 0; i < n_lambda; i++)
     {
-        cmp_real(fp, "inputrec->simtempvals->temperatures", -1, simtemp1->temperatures[i], simtemp2->temperatures[i], ftol, abstol);
+        cmp_real(fp,
+                 "inputrec->simtempvals->temperatures",
+                 -1,
+                 simtemp1->temperatures[i],
+                 simtemp2->temperatures[i],
+                 ftol,
+                 abstol);
     }
 }
 
-static void cmp_expandedvals(FILE *fp, const t_expanded *expand1, const t_expanded *expand2, int n_lambda, real ftol, real abstol)
+static void cmp_expandedvals(FILE*             fp,
+                             const t_expanded* expand1,
+                             const t_expanded* expand2,
+                             int               n_lambda,
+                             real              ftol,
+                             real              abstol)
 {
     int i;
 
@@ -1203,25 +1286,50 @@ static void cmp_expandedvals(FILE *fp, const t_expanded *expand1, const t_expand
 
     for (i = 0; i < n_lambda; i++)
     {
-        cmp_real(fp, "inputrec->expandedvals->init_lambda_weights", -1,
-                 expand1->init_lambda_weights[i], expand2->init_lambda_weights[i], ftol, abstol);
+        cmp_real(fp,
+                 "inputrec->expandedvals->init_lambda_weights",
+                 -1,
+                 expand1->init_lambda_weights[i],
+                 expand2->init_lambda_weights[i],
+                 ftol,
+                 abstol);
     }
 
-    cmp_int(fp, "inputrec->expandedvals->lambda-stats", -1, expand1->elamstats, expand2->elamstats);
-    cmp_int(fp, "inputrec->expandedvals->lambda-mc-move", -1, expand1->elmcmove, expand2->elmcmove);
+    cmpEnum(fp, "inputrec->expandedvals->lambda-stats", expand1->elamstats, expand2->elamstats);
+    cmpEnum(fp, "inputrec->expandedvals->lambda-mc-move", expand1->elmcmove, expand2->elmcmove);
     cmp_int(fp, "inputrec->expandedvals->lmc-repeats", -1, expand1->lmc_repeats, expand2->lmc_repeats);
     cmp_int(fp, "inputrec->expandedvals->lmc-gibbsdelta", -1, expand1->gibbsdeltalam, expand2->gibbsdeltalam);
     cmp_int(fp, "inputrec->expandedvals->lmc-forced-nstart", -1, expand1->lmc_forced_nstart, expand2->lmc_forced_nstart);
-    cmp_int(fp, "inputrec->expandedvals->lambda-weights-equil", -1, expand1->elmceq, expand2->elmceq);
-    cmp_int(fp, "inputrec->expandedvals->,weight-equil-number-all-lambda", -1, expand1->equil_n_at_lam, expand2->equil_n_at_lam);
+    cmpEnum(fp, "inputrec->expandedvals->lambda-weights-equil", expand1->elmceq, expand2->elmceq);
+    cmp_int(fp,
+            "inputrec->expandedvals->,weight-equil-number-all-lambda",
+            -1,
+            expand1->equil_n_at_lam,
+            expand2->equil_n_at_lam);
     cmp_int(fp, "inputrec->expandedvals->weight-equil-number-samples", -1, expand1->equil_samples, expand2->equil_samples);
     cmp_int(fp, "inputrec->expandedvals->weight-equil-number-steps", -1, expand1->equil_steps, expand2->equil_steps);
-    cmp_real(fp, "inputrec->expandedvals->weight-equil-wl-delta", -1, expand1->equil_wl_delta, expand2->equil_wl_delta, ftol, abstol);
-    cmp_real(fp, "inputrec->expandedvals->weight-equil-count-ratio", -1, expand1->equil_ratio, expand2->equil_ratio, ftol, abstol);
-    cmp_bool(fp, "inputrec->expandedvals->symmetrized-transition-matrix", -1, expand1->bSymmetrizedTMatrix, expand2->bSymmetrizedTMatrix);
+    cmp_real(fp,
+             "inputrec->expandedvals->weight-equil-wl-delta",
+             -1,
+             expand1->equil_wl_delta,
+             expand2->equil_wl_delta,
+             ftol,
+             abstol);
+    cmp_real(fp,
+             "inputrec->expandedvals->weight-equil-count-ratio",
+             -1,
+             expand1->equil_ratio,
+             expand2->equil_ratio,
+             ftol,
+             abstol);
+    cmp_bool(fp,
+             "inputrec->expandedvals->symmetrized-transition-matrix",
+             -1,
+             expand1->bSymmetrizedTMatrix,
+             expand2->bSymmetrizedTMatrix);
     cmp_int(fp, "inputrec->expandedvals->nstTij", -1, expand1->nstTij, expand2->nstTij);
     cmp_int(fp, "inputrec->expandedvals->mininum-var-min", -1, expand1->minvarmin, expand2->minvarmin); /*default is reasonable */
-    cmp_int(fp, "inputrec->expandedvals->weight-c-range", -1, expand1->c_range, expand2->c_range);      /* default is just C=0 */
+    cmp_int(fp, "inputrec->expandedvals->weight-c-range", -1, expand1->c_range, expand2->c_range); /* default is just C=0 */
     cmp_real(fp, "inputrec->expandedvals->wl-scale", -1, expand1->wl_scale, expand2->wl_scale, ftol, abstol);
     cmp_real(fp, "inputrec->expandedvals->init-wl-delta", -1, expand1->init_wl_delta, expand2->init_wl_delta, ftol, abstol);
     cmp_real(fp, "inputrec->expandedvals->wl-ratio", -1, expand1->wl_ratio, expand2->wl_ratio, ftol, abstol);
@@ -1230,35 +1338,56 @@ static void cmp_expandedvals(FILE *fp, const t_expanded *expand1, const t_expand
     cmp_real(fp, "inputrec->expandedvals->mc-temperature", -1, expand1->mc_temp, expand2->mc_temp, ftol, abstol);
 }
 
-static void cmp_fepvals(FILE *fp, const t_lambda *fep1, const t_lambda *fep2, real ftol, real abstol)
+static void cmp_fepvals(FILE* fp, const t_lambda* fep1, const t_lambda* fep2, real ftol, real abstol)
 {
     int i, j;
     cmp_int(fp, "inputrec->nstdhdl", -1, fep1->nstdhdl, fep2->nstdhdl);
     cmp_double(fp, "inputrec->fepvals->init_fep_state", -1, fep1->init_fep_state, fep2->init_fep_state, ftol, abstol);
     cmp_double(fp, "inputrec->fepvals->delta_lambda", -1, fep1->delta_lambda, fep2->delta_lambda, ftol, abstol);
     cmp_int(fp, "inputrec->fepvals->n_lambda", -1, fep1->n_lambda, fep2->n_lambda);
-    for (i = 0; i < efptNR; i++)
+    for (i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
     {
         for (j = 0; j < std::min(fep1->n_lambda, fep2->n_lambda); j++)
         {
-            cmp_double(fp, "inputrec->fepvals->all_lambda", -1, fep1->all_lambda[i][j], fep2->all_lambda[i][j], ftol, abstol);
+            cmp_double(fp,
+                       "inputrec->fepvals->all_lambda",
+                       -1,
+                       fep1->all_lambda[i][j],
+                       fep2->all_lambda[i][j],
+                       ftol,
+                       abstol);
         }
     }
-    cmp_int(fp, "inputrec->fepvals->lambda_neighbors", 1, fep1->lambda_neighbors,
-            fep2->lambda_neighbors);
+    cmp_int(fp, "inputrec->fepvals->lambda_neighbors", 1, fep1->lambda_neighbors, fep2->lambda_neighbors);
     cmp_real(fp, "inputrec->fepvals->sc_alpha", -1, fep1->sc_alpha, fep2->sc_alpha, ftol, abstol);
     cmp_int(fp, "inputrec->fepvals->sc_power", -1, fep1->sc_power, fep2->sc_power);
     cmp_real(fp, "inputrec->fepvals->sc_r_power", -1, fep1->sc_r_power, fep2->sc_r_power, ftol, abstol);
     cmp_real(fp, "inputrec->fepvals->sc_sigma", -1, fep1->sc_sigma, fep2->sc_sigma, ftol, abstol);
-    cmp_int(fp, "inputrec->fepvals->edHdLPrintEnergy", -1, fep1->edHdLPrintEnergy, fep1->edHdLPrintEnergy);
+    cmpEnum(fp, "inputrec->fepvals->edHdLPrintEnergy", fep1->edHdLPrintEnergy, fep1->edHdLPrintEnergy);
     cmp_bool(fp, "inputrec->fepvals->bScCoul", -1, fep1->bScCoul, fep1->bScCoul);
-    cmp_int(fp, "inputrec->separate_dhdl_file", -1, fep1->separate_dhdl_file, fep2->separate_dhdl_file);
-    cmp_int(fp, "inputrec->dhdl_derivatives", -1, fep1->dhdl_derivatives, fep2->dhdl_derivatives);
+    cmpEnum(fp, "inputrec->separate_dhdl_file", fep1->separate_dhdl_file, fep2->separate_dhdl_file);
+    cmpEnum(fp, "inputrec->dhdl_derivatives", fep1->dhdl_derivatives, fep2->dhdl_derivatives);
     cmp_int(fp, "inputrec->dh_hist_size", -1, fep1->dh_hist_size, fep2->dh_hist_size);
     cmp_double(fp, "inputrec->dh_hist_spacing", -1, fep1->dh_hist_spacing, fep2->dh_hist_spacing, ftol, abstol);
+    cmpEnum(fp, "inputrec->fepvals->softcoreFunction", fep1->softcoreFunction, fep2->softcoreFunction);
+    cmp_real(fp,
+             "inputrec->fepvals->scGapsysScaleLinpointLJ",
+             -1,
+             fep1->scGapsysScaleLinpointLJ,
+             fep2->scGapsysScaleLinpointLJ,
+             ftol,
+             abstol);
+    cmp_real(fp,
+             "inputrec->fepvals->scGapsysScaleLinpointQ",
+             -1,
+             fep1->scGapsysScaleLinpointQ,
+             fep2->scGapsysScaleLinpointQ,
+             ftol,
+             abstol);
+    cmp_real(fp, "inputrec->fepvals->scGapsysSigmaLJ", -1, fep1->scGapsysSigmaLJ, fep2->scGapsysSigmaLJ, ftol, abstol);
 }
 
-void cmp_inputrec(FILE *fp, const t_inputrec *ir1, const t_inputrec *ir2, real ftol, real abstol)
+void cmp_inputrec(FILE* fp, const t_inputrec* ir1, const t_inputrec* ir2, real ftol, real abstol)
 {
     fprintf(fp, "comparing inputrec\n");
 
@@ -1269,17 +1398,31 @@ void cmp_inputrec(FILE *fp, const t_inputrec *ir1, const t_inputrec *ir2, real f
      * #define CII(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
      * #define CIR(s) cmp_real(fp,"inputrec->"#s,0,ir1->##s,ir2->##s,ftol)
      */
-    cmp_int(fp, "inputrec->eI", -1, ir1->eI, ir2->eI);
+    cmpEnum(fp, "inputrec->eI", ir1->eI, ir2->eI);
     cmp_int64(fp, "inputrec->nsteps", ir1->nsteps, ir2->nsteps);
     cmp_int64(fp, "inputrec->init_step", ir1->init_step, ir2->init_step);
     cmp_int(fp, "inputrec->simulation_part", -1, ir1->simulation_part, ir2->simulation_part);
-    cmp_int(fp, "inputrec->ePBC", -1, ir1->ePBC, ir2->ePBC);
-    cmp_int(fp, "inputrec->bPeriodicMols", -1, ir1->bPeriodicMols, ir2->bPeriodicMols);
-    cmp_int(fp, "inputrec->cutoff_scheme", -1, ir1->cutoff_scheme, ir2->cutoff_scheme);
-    cmp_int(fp, "inputrec->ns_type", -1, ir1->ns_type, ir2->ns_type);
+    cmp_int(fp, "inputrec->mts", -1, static_cast<int>(ir1->useMts), static_cast<int>(ir2->useMts));
+    if (ir1->useMts && ir2->useMts)
+    {
+        cmp_int(fp, "inputrec->mts-levels", -1, ir1->mtsLevels.size(), ir2->mtsLevels.size());
+        cmp_int(fp,
+                "inputrec->mts-level2-forces",
+                -1,
+                ir1->mtsLevels[1].forceGroups.to_ulong(),
+                ir2->mtsLevels[1].forceGroups.to_ulong());
+        cmp_int(fp,
+                "inputrec->mts-level2-factor",
+                -1,
+                ir1->mtsLevels[1].stepFactor,
+                ir2->mtsLevels[1].stepFactor);
+    }
+    cmp_int(fp, "inputrec->pbcType", -1, static_cast<int>(ir1->pbcType), static_cast<int>(ir2->pbcType));
+    cmp_bool(fp, "inputrec->bPeriodicMols", -1, ir1->bPeriodicMols, ir2->bPeriodicMols);
+    cmpEnum(fp, "inputrec->cutoff_scheme", ir1->cutoff_scheme, ir2->cutoff_scheme);
     cmp_int(fp, "inputrec->nstlist", -1, ir1->nstlist, ir2->nstlist);
     cmp_int(fp, "inputrec->nstcomm", -1, ir1->nstcomm, ir2->nstcomm);
-    cmp_int(fp, "inputrec->comm_mode", -1, ir1->comm_mode, ir2->comm_mode);
+    cmpEnum(fp, "inputrec->comm_mode", ir1->comm_mode, ir2->comm_mode);
     cmp_int(fp, "inputrec->nstlog", -1, ir1->nstlog, ir2->nstlog);
     cmp_int(fp, "inputrec->nstxout", -1, ir1->nstxout, ir2->nstxout);
     cmp_int(fp, "inputrec->nstvout", -1, ir1->nstvout, ir2->nstvout);
@@ -1289,21 +1432,35 @@ void cmp_inputrec(FILE *fp, const t_inputrec *ir1, const t_inputrec *ir2, real f
     cmp_int(fp, "inputrec->nstxout_compressed", -1, ir1->nstxout_compressed, ir2->nstxout_compressed);
     cmp_double(fp, "inputrec->init_t", -1, ir1->init_t, ir2->init_t, ftol, abstol);
     cmp_double(fp, "inputrec->delta_t", -1, ir1->delta_t, ir2->delta_t, ftol, abstol);
-    cmp_real(fp, "inputrec->x_compression_precision", -1, ir1->x_compression_precision, ir2->x_compression_precision, ftol, abstol);
+    cmp_real(fp,
+             "inputrec->x_compression_precision",
+             -1,
+             ir1->x_compression_precision,
+             ir2->x_compression_precision,
+             ftol,
+             abstol);
     cmp_real(fp, "inputrec->fourierspacing", -1, ir1->fourier_spacing, ir2->fourier_spacing, ftol, abstol);
     cmp_int(fp, "inputrec->nkx", -1, ir1->nkx, ir2->nkx);
     cmp_int(fp, "inputrec->nky", -1, ir1->nky, ir2->nky);
     cmp_int(fp, "inputrec->nkz", -1, ir1->nkz, ir2->nkz);
     cmp_int(fp, "inputrec->pme_order", -1, ir1->pme_order, ir2->pme_order);
     cmp_real(fp, "inputrec->ewald_rtol", -1, ir1->ewald_rtol, ir2->ewald_rtol, ftol, abstol);
-    cmp_int(fp, "inputrec->ewald_geometry", -1, ir1->ewald_geometry, ir2->ewald_geometry);
+    cmpEnum(fp, "inputrec->ewald_geometry", ir1->ewald_geometry, ir2->ewald_geometry);
     cmp_real(fp, "inputrec->epsilon_surface", -1, ir1->epsilon_surface, ir2->epsilon_surface, ftol, abstol);
-    cmp_int(fp, "inputrec->bContinuation", -1, ir1->bContinuation, ir2->bContinuation);
-    cmp_int(fp, "inputrec->bShakeSOR", -1, ir1->bShakeSOR, ir2->bShakeSOR);
-    cmp_int(fp, "inputrec->etc", -1, ir1->etc, ir2->etc);
-    cmp_int(fp, "inputrec->bPrintNHChains", -1, ir1->bPrintNHChains, ir2->bPrintNHChains);
-    cmp_int(fp, "inputrec->epc", -1, ir1->epc, ir2->epc);
-    cmp_int(fp, "inputrec->epct", -1, ir1->epct, ir2->epct);
+    cmp_int(fp,
+            "inputrec->bContinuation",
+            -1,
+            static_cast<int>(ir1->bContinuation),
+            static_cast<int>(ir2->bContinuation));
+    cmp_int(fp, "inputrec->bShakeSOR", -1, static_cast<int>(ir1->bShakeSOR), static_cast<int>(ir2->bShakeSOR));
+    cmpEnum(fp, "inputrec->etc", ir1->etc, ir2->etc);
+    cmp_int(fp,
+            "inputrec->bPrintNHChains",
+            -1,
+            static_cast<int>(ir1->bPrintNHChains),
+            static_cast<int>(ir2->bPrintNHChains));
+    cmpEnum(fp, "inputrec->epc", ir1->epc, ir2->epc);
+    cmpEnum(fp, "inputrec->epct", ir1->epct, ir2->epct);
     cmp_real(fp, "inputrec->tau_p", -1, ir1->tau_p, ir2->tau_p, ftol, abstol);
     cmp_rvec(fp, "inputrec->ref_p(x)", -1, ir1->ref_p[XX], ir2->ref_p[XX], ftol, abstol);
     cmp_rvec(fp, "inputrec->ref_p(y)", -1, ir1->ref_p[YY], ir2->ref_p[YY], ftol, abstol);
@@ -1311,39 +1468,50 @@ void cmp_inputrec(FILE *fp, const t_inputrec *ir1, const t_inputrec *ir2, real f
     cmp_rvec(fp, "inputrec->compress(x)", -1, ir1->compress[XX], ir2->compress[XX], ftol, abstol);
     cmp_rvec(fp, "inputrec->compress(y)", -1, ir1->compress[YY], ir2->compress[YY], ftol, abstol);
     cmp_rvec(fp, "inputrec->compress(z)", -1, ir1->compress[ZZ], ir2->compress[ZZ], ftol, abstol);
-    cmp_int(fp, "refcoord_scaling", -1, ir1->refcoord_scaling, ir2->refcoord_scaling);
+    cmpEnum(fp, "refcoord_scaling", ir1->refcoord_scaling, ir2->refcoord_scaling);
     cmp_rvec(fp, "inputrec->posres_com", -1, ir1->posres_com, ir2->posres_com, ftol, abstol);
     cmp_rvec(fp, "inputrec->posres_comB", -1, ir1->posres_comB, ir2->posres_comB, ftol, abstol);
     cmp_real(fp, "inputrec->verletbuf_tol", -1, ir1->verletbuf_tol, ir2->verletbuf_tol, ftol, abstol);
     cmp_real(fp, "inputrec->rlist", -1, ir1->rlist, ir2->rlist, ftol, abstol);
     cmp_real(fp, "inputrec->rtpi", -1, ir1->rtpi, ir2->rtpi, ftol, abstol);
-    cmp_int(fp, "inputrec->coulombtype", -1, ir1->coulombtype, ir2->coulombtype);
-    cmp_int(fp, "inputrec->coulomb_modifier", -1, ir1->coulomb_modifier, ir2->coulomb_modifier);
+    cmpEnum(fp, "inputrec->coulombtype", ir1->coulombtype, ir2->coulombtype);
+    cmpEnum(fp, "inputrec->coulomb_modifier", ir1->coulomb_modifier, ir2->coulomb_modifier);
     cmp_real(fp, "inputrec->rcoulomb_switch", -1, ir1->rcoulomb_switch, ir2->rcoulomb_switch, ftol, abstol);
     cmp_real(fp, "inputrec->rcoulomb", -1, ir1->rcoulomb, ir2->rcoulomb, ftol, abstol);
-    cmp_int(fp, "inputrec->vdwtype", -1, ir1->vdwtype, ir2->vdwtype);
-    cmp_int(fp, "inputrec->vdw_modifier", -1, ir1->vdw_modifier, ir2->vdw_modifier);  cmp_real(fp, "inputrec->rvdw_switch", -1, ir1->rvdw_switch, ir2->rvdw_switch, ftol, abstol);
+    cmpEnum(fp, "inputrec->vdwtype", ir1->vdwtype, ir2->vdwtype);
+    cmpEnum(fp, "inputrec->vdw_modifier", ir1->vdw_modifier, ir2->vdw_modifier);
+    cmp_real(fp, "inputrec->rvdw_switch", -1, ir1->rvdw_switch, ir2->rvdw_switch, ftol, abstol);
     cmp_real(fp, "inputrec->rvdw", -1, ir1->rvdw, ir2->rvdw, ftol, abstol);
     cmp_real(fp, "inputrec->epsilon_r", -1, ir1->epsilon_r, ir2->epsilon_r, ftol, abstol);
     cmp_real(fp, "inputrec->epsilon_rf", -1, ir1->epsilon_rf, ir2->epsilon_rf, ftol, abstol);
     cmp_real(fp, "inputrec->tabext", -1, ir1->tabext, ir2->tabext, ftol, abstol);
 
-    cmp_int(fp, "inputrec->eDispCorr", -1, ir1->eDispCorr, ir2->eDispCorr);
+    cmpEnum(fp, "inputrec->eDispCorr", ir1->eDispCorr, ir2->eDispCorr);
     cmp_real(fp, "inputrec->shake_tol", -1, ir1->shake_tol, ir2->shake_tol, ftol, abstol);
-    cmp_int(fp, "inputrec->efep", -1, ir1->efep, ir2->efep);
-    cmp_fepvals(fp, ir1->fepvals, ir2->fepvals, ftol, abstol);
-    cmp_int(fp, "inputrec->bSimTemp", -1, ir1->bSimTemp, ir2->bSimTemp);
+    cmpEnum(fp, "inputrec->efep", ir1->efep, ir2->efep);
+    cmp_fepvals(fp, ir1->fepvals.get(), ir2->fepvals.get(), ftol, abstol);
+    cmp_int(fp, "inputrec->bSimTemp", -1, static_cast<int>(ir1->bSimTemp), static_cast<int>(ir2->bSimTemp));
     if ((ir1->bSimTemp == ir2->bSimTemp) && (ir1->bSimTemp))
     {
-        cmp_simtempvals(fp, ir1->simtempvals, ir2->simtempvals, std::min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda), ftol, abstol);
+        cmp_simtempvals(fp,
+                        ir1->simtempvals.get(),
+                        ir2->simtempvals.get(),
+                        std::min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda),
+                        ftol,
+                        abstol);
     }
-    cmp_int(fp, "inputrec->bExpanded", -1, ir1->bExpanded, ir2->bExpanded);
+    cmp_int(fp, "inputrec->bExpanded", -1, static_cast<int>(ir1->bExpanded), static_cast<int>(ir2->bExpanded));
     if ((ir1->bExpanded == ir2->bExpanded) && (ir1->bExpanded))
     {
-        cmp_expandedvals(fp, ir1->expandedvals, ir2->expandedvals, std::min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda), ftol, abstol);
+        cmp_expandedvals(fp,
+                         ir1->expandedvals.get(),
+                         ir2->expandedvals.get(),
+                         std::min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda),
+                         ftol,
+                         abstol);
     }
     cmp_int(fp, "inputrec->nwall", -1, ir1->nwall, ir2->nwall);
-    cmp_int(fp, "inputrec->wall_type", -1, ir1->wall_type, ir2->wall_type);
+    cmpEnum(fp, "inputrec->wall_type", ir1->wall_type, ir2->wall_type);
     cmp_int(fp, "inputrec->wall_atomtype[0]", -1, ir1->wall_atomtype[0], ir2->wall_atomtype[0]);
     cmp_int(fp, "inputrec->wall_atomtype[1]", -1, ir1->wall_atomtype[1], ir2->wall_atomtype[1]);
     cmp_real(fp, "inputrec->wall_density[0]", -1, ir1->wall_density[0], ir2->wall_density[0], ftol, abstol);
@@ -1359,13 +1527,13 @@ void cmp_inputrec(FILE *fp, const t_inputrec *ir1, const t_inputrec *ir2, real f
     cmp_bool(fp, "inputrec->bDoAwh", -1, ir1->bDoAwh, ir2->bDoAwh);
     if (ir1->bDoAwh && ir2->bDoAwh)
     {
-        cmp_awhParams(fp, ir1->awhParams, ir2->awhParams, ftol, abstol);
+        cmp_awhParams(fp, *ir1->awhParams, *ir2->awhParams, ftol, abstol);
     }
 
-    cmp_int(fp, "inputrec->eDisre", -1, ir1->eDisre, ir2->eDisre);
+    cmpEnum(fp, "inputrec->eDisre", ir1->eDisre, ir2->eDisre);
     cmp_real(fp, "inputrec->dr_fc", -1, ir1->dr_fc, ir2->dr_fc, ftol, abstol);
-    cmp_int(fp, "inputrec->eDisreWeighting", -1, ir1->eDisreWeighting, ir2->eDisreWeighting);
-    cmp_int(fp, "inputrec->bDisreMixed", -1, ir1->bDisreMixed, ir2->bDisreMixed);
+    cmpEnum(fp, "inputrec->eDisreWeighting", ir1->eDisreWeighting, ir2->eDisreWeighting);
+    cmp_int(fp, "inputrec->bDisreMixed", -1, static_cast<int>(ir1->bDisreMixed), static_cast<int>(ir2->bDisreMixed));
     cmp_int(fp, "inputrec->nstdisreout", -1, ir1->nstdisreout, ir2->nstdisreout);
     cmp_real(fp, "inputrec->dr_tau", -1, ir1->dr_tau, ir2->dr_tau, ftol, abstol);
     cmp_real(fp, "inputrec->orires_fc", -1, ir1->orires_fc, ir2->orires_fc, ftol, abstol);
@@ -1377,7 +1545,7 @@ void cmp_inputrec(FILE *fp, const t_inputrec *ir1, const t_inputrec *ir2, real f
     cmp_real(fp, "inputrec->fc_stepsize", -1, ir1->fc_stepsize, ir2->fc_stepsize, ftol, abstol);
     cmp_int(fp, "inputrec->nstcgsteep", -1, ir1->nstcgsteep, ir2->nstcgsteep);
     cmp_int(fp, "inputrec->nbfgscorr", 0, ir1->nbfgscorr, ir2->nbfgscorr);
-    cmp_int(fp, "inputrec->eConstrAlg", -1, ir1->eConstrAlg, ir2->eConstrAlg);
+    cmpEnum(fp, "inputrec->eConstrAlg", ir1->eConstrAlg, ir2->eConstrAlg);
     cmp_int(fp, "inputrec->nProjOrder", -1, ir1->nProjOrder, ir2->nProjOrder);
     cmp_real(fp, "inputrec->LincsWarnAngle", -1, ir1->LincsWarnAngle, ir2->LincsWarnAngle, ftol, abstol);
     cmp_int(fp, "inputrec->nLincsIter", -1, ir1->nLincsIter, ir2->nLincsIter);
@@ -1402,76 +1570,82 @@ void cmp_inputrec(FILE *fp, const t_inputrec *ir1, const t_inputrec *ir2, real f
     gmx::compareKeyValueTrees(&writer, *ir1->params, *ir2->params, ftol, abstol);
 }
 
-void comp_pull_AB(FILE *fp, pull_params_t *pull, real ftol, real abstol)
+void comp_pull_AB(FILE* fp, const pull_params_t& pull, real ftol, real abstol)
 {
-    int i;
-
-    for (i = 0; i < pull->ncoord; i++)
+    for (int i = 0; i < pull.ncoord; i++)
     {
         fprintf(fp, "comparing pull coord %d\n", i);
-        cmp_real(fp, "pull-coord->k", -1, pull->coord[i].k, pull->coord[i].kB, ftol, abstol);
+        cmp_real(fp, "pull-coord->k", -1, pull.coord[i].k, pull.coord[i].kB, ftol, abstol);
     }
 }
 
-gmx_bool inputrecDeform(const t_inputrec *ir)
+gmx_bool inputrecDeform(const t_inputrecir)
 {
-    return (ir->deform[XX][XX] != 0 || ir->deform[YY][YY] != 0 || ir->deform[ZZ][ZZ] != 0 ||
-            ir->deform[YY][XX] != 0 || ir->deform[ZZ][XX] != 0 || ir->deform[ZZ][YY] != 0);
+    return (ir->deform[XX][XX] != 0 || ir->deform[YY][YY] != 0 || ir->deform[ZZ][ZZ] != 0
+            || ir->deform[YY][XX] != 0 || ir->deform[ZZ][XX] != 0 || ir->deform[ZZ][YY] != 0);
 }
 
-gmx_bool inputrecDynamicBox(const t_inputrec *ir)
+gmx_bool inputrecDynamicBox(const t_inputrecir)
 {
-    return (ir->epc != epcNO || ir->eI == eiTPI || inputrecDeform(ir));
+    return (ir->epc != PressureCoupling::No || ir->eI == IntegrationAlgorithm::TPI || inputrecDeform(ir));
 }
 
-gmx_bool inputrecPreserveShape(const t_inputrec *ir)
+gmx_bool inputrecPreserveShape(const t_inputrecir)
 {
-    return  (ir->epc != epcNO && ir->deform[XX][XX] == 0 &&
-             (ir->epct == epctISOTROPIC || ir->epct == epctSEMIISOTROPIC));
+    return (ir->epc != PressureCoupling::No && ir->deform[XX][XX] == 0
+            && (ir->epct == PressureCouplingType::Isotropic
+                || ir->epct == PressureCouplingType::SemiIsotropic));
 }
 
-gmx_bool inputrecNeedMutot(const t_inputrec *ir)
+gmx_bool inputrecNeedMutot(const t_inputrecir)
 {
-    return ((ir->coulombtype == eelEWALD || EEL_PME(ir->coulombtype)) &&
-            (ir->ewald_geometry == eewg3DC || ir->epsilon_surface != 0));
+    return ((ir->coulombtype == CoulombInteractionType::Ewald || EEL_PME(ir->coulombtype))
+            && (ir->ewald_geometry == EwaldGeometry::ThreeDC || ir->epsilon_surface != 0));
 }
 
-gmx_bool inputrecExclForces(const t_inputrec *ir)
+gmx_bool inputrecExclForces(const t_inputrecir)
 {
     return (EEL_FULL(ir->coulombtype) || (EEL_RF(ir->coulombtype)));
 }
 
-gmx_bool inputrecNptTrotter(const t_inputrec *ir)
+gmx_bool inputrecNptTrotter(const t_inputrecir)
 {
-    return ( ( (ir->eI == eiVV) || (ir->eI == eiVVAK) ) &&
-             (ir->epc == epcMTTK) && (ir->etc == etcNOSEHOOVER) );
+    return (((ir->eI == IntegrationAlgorithm::VV) || (ir->eI == IntegrationAlgorithm::VVAK))
+            && (ir->epc == PressureCoupling::Mttk) && (ir->etc == TemperatureCoupling::NoseHoover));
 }
 
-gmx_bool inputrecNvtTrotter(const t_inputrec *ir)
+gmx_bool inputrecNvtTrotter(const t_inputrecir)
 {
-    return ( ( (ir->eI == eiVV) || (ir->eI == eiVVAK) ) &&
-             (ir->epc != epcMTTK) && (ir->etc == etcNOSEHOOVER) );
+    return (((ir->eI == IntegrationAlgorithm::VV) || (ir->eI == IntegrationAlgorithm::VVAK))
+            && (ir->epc != PressureCoupling::Mttk) && (ir->etc == TemperatureCoupling::NoseHoover));
 }
 
-gmx_bool inputrecNphTrotter(const t_inputrec *ir)
+gmx_bool inputrecNphTrotter(const t_inputrecir)
 {
-    return ( ( (ir->eI == eiVV) || (ir->eI == eiVVAK) ) &&
-             (ir->epc == epcMTTK) && (ir->etc != etcNOSEHOOVER) );
+    return (((ir->eI == IntegrationAlgorithm::VV) || (ir->eI == IntegrationAlgorithm::VVAK))
+            && (ir->epc == PressureCoupling::Mttk) && (ir->etc != TemperatureCoupling::NoseHoover));
 }
 
-bool inputrecPbcXY2Walls(const t_inputrec *ir)
+bool inputrecPbcXY2Walls(const t_inputrecir)
 {
-    return (ir->ePBC == epbcXY && ir->nwall == 2);
+    return (ir->pbcType == PbcType::XY && ir->nwall == 2);
 }
 
-bool integratorHasConservedEnergyQuantity(const t_inputrec *ir)
+bool inputrecFrozenAtoms(const t_inputrec* ir)
+{
+    return ((ir->opts.nFreeze != nullptr)
+            && (ir->opts.ngfrz > 1 || ir->opts.nFreeze[0][XX] != 0 || ir->opts.nFreeze[0][YY] != 0
+                || ir->opts.nFreeze[0][ZZ] != 0));
+}
+
+bool integratorHasConservedEnergyQuantity(const t_inputrec* ir)
 {
     if (!EI_MD(ir->eI))
-    {
+    { // NOLINT bugprone-branch-clone
         // Energy minimization or stochastic integrator: no conservation
         return false;
     }
-    else if (ir->etc == etcNO && ir->epc == epcNO)
+    else if (ir->etc == TemperatureCoupling::No && ir->epc == PressureCoupling::No)
     {
         // The total energy is conserved, no additional conserved quanitity
         return false;
@@ -1480,14 +1654,20 @@ bool integratorHasConservedEnergyQuantity(const t_inputrec *ir)
     {
         // Shear stress with Parrinello-Rahman is not supported (tedious)
         bool shearWithPR =
-            ((ir->epc == epcPARRINELLORAHMAN || ir->epc == epcMTTK) &&
-             (ir->ref_p[YY][XX] != 0 || ir->ref_p[ZZ][XX] != 0 || ir->ref_p[ZZ][YY] != 0));
+                ((ir->epc == PressureCoupling::ParrinelloRahman || ir->epc == PressureCoupling::Mttk)
+                 && (ir->ref_p[YY][XX] != 0 || ir->ref_p[ZZ][XX] != 0 || ir->ref_p[ZZ][YY] != 0));
 
         return !ETC_ANDERSEN(ir->etc) && !shearWithPR;
     }
 }
 
-int inputrec2nboundeddim(const t_inputrec *ir)
+bool integratorHasReferenceTemperature(const t_inputrec* ir)
+{
+    return ((ir->etc != TemperatureCoupling::No) || EI_SD(ir->eI)
+            || (ir->eI == IntegrationAlgorithm::BD) || EI_TPI(ir->eI));
+}
+
+int inputrec2nboundeddim(const t_inputrec* ir)
 {
     if (inputrecPbcXY2Walls(ir))
     {
@@ -1495,29 +1675,68 @@ int inputrec2nboundeddim(const t_inputrec *ir)
     }
     else
     {
-        return ePBC2npbcdim(ir->ePBC);
+        return numPbcDimensions(ir->pbcType);
     }
 }
 
-int ndof_com(const t_inputrec *ir)
+int ndof_com(const t_inputrecir)
 {
     int n = 0;
 
-    switch (ir->ePBC)
+    switch (ir->pbcType)
     {
-        case epbcXYZ:
-        case epbcNONE:
-            n = 3;
-            break;
-        case epbcXY:
-            n = (ir->nwall == 0 ? 3 : 2);
-            break;
-        case epbcSCREW:
-            n = 1;
-            break;
-        default:
-            gmx_incons("Unknown pbc in calc_nrdf");
+        case PbcType::Xyz:
+        case PbcType::No: n = 3; break;
+        case PbcType::XY: n = (ir->nwall == 0 ? 3 : 2); break;
+        case PbcType::Screw: n = 1; break;
+        default: gmx_incons("Unknown pbc in calc_nrdf");
     }
 
     return n;
 }
+
+real maxReferenceTemperature(const t_inputrec& ir)
+{
+    if (EI_ENERGY_MINIMIZATION(ir.eI) || ir.eI == IntegrationAlgorithm::NM)
+    {
+        return 0;
+    }
+
+    if (EI_MD(ir.eI) && ir.etc == TemperatureCoupling::No)
+    {
+        return -1;
+    }
+
+    /* SD and BD also use ref_t and tau_t for setting the reference temperature.
+     * TPI can be treated as MD, since it needs an ensemble temperature.
+     */
+
+    real maxTemperature = 0;
+    for (int i = 0; i < ir.opts.ngtc; i++)
+    {
+        if (ir.opts.tau_t[i] >= 0)
+        {
+            maxTemperature = std::max(maxTemperature, ir.opts.ref_t[i]);
+        }
+    }
+
+    return maxTemperature;
+}
+
+bool haveEwaldSurfaceContribution(const t_inputrec& ir)
+{
+    return EEL_PME_EWALD(ir.coulombtype)
+           && (ir.ewald_geometry == EwaldGeometry::ThreeDC || ir.epsilon_surface != 0);
+}
+
+bool haveFreeEnergyType(const t_inputrec& ir, const int fepType)
+{
+    for (int i = 0; i < ir.fepvals->n_lambda; i++)
+    {
+        if (ir.fepvals->all_lambda[fepType][i] > 0)
+        {
+            return true;
+        }
+    }
+    return false;
+}