Refactor md_enums
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / readir.cpp
index 6595bff142eefec8798a74cefad4c2d81c450c9b..2276b84de2ece6a26ab28878eefa816d00017318 100644 (file)
@@ -37,6 +37,7 @@
  */
 #include "gmxpre.h"
 
+#include "gromacs/utility/enumerationhelpers.h"
 #include "readir.h"
 
 #include <cctype>
@@ -107,10 +108,10 @@ struct gmx_inputrec_strings
             energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
             couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
             wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN], imd_grp[STRLEN];
-    char                     fep_lambda[efptNR][STRLEN];
-    char                     lambda_weights[STRLEN];
-    std::vector<std::string> pullGroupNames;
-    std::vector<std::string> rotateGroupNames;
+    gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, std::string> fep_lambda;
+    char                                                                   lambda_weights[STRLEN];
+    std::vector<std::string>                                               pullGroupNames;
+    std::vector<std::string>                                               rotateGroupNames;
     char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
 };
 
@@ -150,7 +151,7 @@ static const char* constraints[eshNR + 1] = { "none",     "h-bonds",    "all-bon
 
 static const char* couple_lam[ecouplamNR + 1] = { "vdw-q", "vdw", "q", "none", nullptr };
 
-static void GetSimTemps(int ntemps, t_simtemp* simtemp, double* temperature_lambdas)
+static void getSimTemps(int ntemps, t_simtemp* simtemp, gmx::ArrayRef<double> temperature_lambdas)
 {
 
     int i;
@@ -158,20 +159,20 @@ static void GetSimTemps(int ntemps, t_simtemp* simtemp, double* temperature_lamb
     for (i = 0; i < ntemps; i++)
     {
         /* simple linear scaling -- allows more control */
-        if (simtemp->eSimTempScale == esimtempLINEAR)
+        if (simtemp->eSimTempScale == SimulatedTempering::Linear)
         {
             simtemp->temperatures[i] =
                     simtemp->simtemp_low
                     + (simtemp->simtemp_high - simtemp->simtemp_low) * temperature_lambdas[i];
         }
         else if (simtemp->eSimTempScale
-                 == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
+                 == SimulatedTempering::Geometric) /* should give roughly equal acceptance for constant heat capacity . . . */
         {
             simtemp->temperatures[i] = simtemp->simtemp_low
                                        * std::pow(simtemp->simtemp_high / simtemp->simtemp_low,
                                                   static_cast<real>((1.0 * i) / (ntemps - 1)));
         }
-        else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
+        else if (simtemp->eSimTempScale == SimulatedTempering::Exponential)
         {
             simtemp->temperatures[i] = simtemp->simtemp_low
                                        + (simtemp->simtemp_high - simtemp->simtemp_low)
@@ -180,7 +181,7 @@ static void GetSimTemps(int ntemps, t_simtemp* simtemp, double* temperature_lamb
         else
         {
             char errorstr[128];
-            sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
+            sprintf(errorstr, "eSimTempScale=%s not defined", enumValueToString(simtemp->eSimTempScale));
             gmx_fatal(FARGS, "%s", errorstr);
         }
     }
@@ -209,11 +210,11 @@ static void check_nst(const char* desc_nst, int nst, const char* desc_p, int* p,
 }
 
 //! Convert legacy mdp entries to modern ones.
-static void process_interaction_modifier(int* eintmod)
+static void process_interaction_modifier(InteractionModifiers* eintmod)
 {
-    if (*eintmod == eintmodPOTSHIFT_VERLET_UNSUPPORTED)
+    if (*eintmod == InteractionModifiers::PotShiftVerletUnsupported)
     {
-        *eintmod = eintmodPOTSHIFT;
+        *eintmod = InteractionModifiers::PotShift;
     }
 }
 
@@ -235,8 +236,8 @@ void check_ir(const char*                   mdparin,
     char        err_buf[256], warn_buf[STRLEN];
     int         i, j;
     real        dt_pcoupl;
-    t_lambda*   fep    = ir->fepvals;
-    t_expanded* expand = ir->expandedvals;
+    t_lambda*   fep    = ir->fepvals.get();
+    t_expanded* expand = ir->expandedvals.get();
 
     set_warning_line(wi, mdparin, -1);
 
@@ -253,10 +254,12 @@ void check_ir(const char*                   mdparin,
         }
     }
 
-    if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
+    if (ir->coulombtype == CoulombInteractionType::RFNecUnsupported)
     {
-        sprintf(warn_buf, "%s electrostatics is no longer supported", eel_names[eelRF_NEC_UNSUPPORTED]);
-        warning_error(wi, warn_buf);
+        std::string message =
+                gmx::formatString("%s electrostatics is no longer supported",
+                                  enumValueToString(CoulombInteractionType::RFNecUnsupported));
+        warning_error(wi, message);
     }
 
     /* BASIC CUT-OFF STUFF */
@@ -268,7 +271,7 @@ void check_ir(const char*                   mdparin,
     {
         warning_error(wi, "rvdw should be >= 0");
     }
-    if (ir->rlist < 0 && !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
+    if (ir->rlist < 0 && !(ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol > 0))
     {
         warning_error(wi, "rlist should be >= 0");
     }
@@ -281,13 +284,13 @@ void check_ir(const char*                   mdparin,
     process_interaction_modifier(&ir->coulomb_modifier);
     process_interaction_modifier(&ir->vdw_modifier);
 
-    if (ir->cutoff_scheme == ecutsGROUP)
+    if (ir->cutoff_scheme == CutoffScheme::Group)
     {
         gmx_fatal(FARGS,
                   "The group cutoff scheme has been removed since GROMACS 2020. "
                   "Please use the Verlet cutoff scheme.");
     }
-    if (ir->cutoff_scheme == ecutsVERLET)
+    if (ir->cutoff_scheme == CutoffScheme::Verlet)
     {
         real rc_max;
 
@@ -302,8 +305,9 @@ void check_ir(const char*                   mdparin,
         {
             // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
             // for PME load balancing, we can support this exception.
-            bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == evdwCUT
-                                            && ir->rcoulomb > ir->rvdw);
+            bool bUsesPmeTwinRangeKernel =
+                    (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == VanDerWaalsType::Cut
+                     && ir->rcoulomb > ir->rvdw);
             if (!bUsesPmeTwinRangeKernel)
             {
                 warning_error(wi,
@@ -312,47 +316,51 @@ void check_ir(const char*                   mdparin,
             }
         }
 
-        if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
+        if (ir->vdwtype == VanDerWaalsType::Shift || ir->vdwtype == VanDerWaalsType::Switch)
         {
-            if (ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT)
+            if (ir->vdw_modifier == InteractionModifiers::None
+                || ir->vdw_modifier == InteractionModifiers::PotShift)
             {
-                ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
+                ir->vdw_modifier =
+                        (ir->vdwtype == VanDerWaalsType::Shift ? InteractionModifiers::ForceSwitch
+                                                               : InteractionModifiers::PotSwitch);
 
                 sprintf(warn_buf,
                         "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and "
                         "vdw_modifier=%s",
-                        evdw_names[ir->vdwtype],
-                        evdw_names[evdwCUT],
-                        eintmod_names[ir->vdw_modifier]);
+                        enumValueToString(ir->vdwtype),
+                        enumValueToString(VanDerWaalsType::Cut),
+                        enumValueToString(ir->vdw_modifier));
                 warning_note(wi, warn_buf);
 
-                ir->vdwtype = evdwCUT;
+                ir->vdwtype = VanDerWaalsType::Cut;
             }
             else
             {
                 sprintf(warn_buf,
                         "Unsupported combination of vdwtype=%s and vdw_modifier=%s",
-                        evdw_names[ir->vdwtype],
-                        eintmod_names[ir->vdw_modifier]);
+                        enumValueToString(ir->vdwtype),
+                        enumValueToString(ir->vdw_modifier));
                 warning_error(wi, warn_buf);
             }
         }
 
-        if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
+        if (!(ir->vdwtype == VanDerWaalsType::Cut || ir->vdwtype == VanDerWaalsType::Pme))
         {
             warning_error(wi,
                           "With Verlet lists only cut-off and PME LJ interactions are supported");
         }
-        if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) || EEL_PME(ir->coulombtype)
-              || ir->coulombtype == eelEWALD))
+        if (!(ir->coulombtype == CoulombInteractionType::Cut || EEL_RF(ir->coulombtype)
+              || EEL_PME(ir->coulombtype) || ir->coulombtype == CoulombInteractionType::Ewald))
         {
             warning_error(wi,
                           "With Verlet lists only cut-off, reaction-field, PME and Ewald "
                           "electrostatics are supported");
         }
-        if (!(ir->coulomb_modifier == eintmodNONE || ir->coulomb_modifier == eintmodPOTSHIFT))
+        if (!(ir->coulomb_modifier == InteractionModifiers::None
+              || ir->coulomb_modifier == InteractionModifiers::PotShift))
         {
-            sprintf(warn_buf, "coulomb_modifier=%s is not supported", eintmod_names[ir->coulomb_modifier]);
+            sprintf(warn_buf, "coulomb_modifier=%s is not supported", enumValueToString(ir->coulomb_modifier));
             warning_error(wi, warn_buf);
         }
 
@@ -360,7 +368,7 @@ void check_ir(const char*                   mdparin,
         {
             sprintf(warn_buf,
                     "Coulomb type %s is not supported with the verlet scheme",
-                    eel_names[ir->coulombtype]);
+                    enumValueToString(ir->coulombtype));
             warning_error(wi, warn_buf);
         }
 
@@ -459,8 +467,8 @@ void check_ir(const char*                   mdparin,
                         "implicitly. See the documentation for more information on which "
                         "parameters affect temperature for %s.",
                         enumValueToString(ir->etc),
-                        ei_names[ir->eI],
-                        ei_names[ir->eI]);
+                        enumValueToString(ir->eI),
+                        enumValueToString(ir->eI));
             }
             else
             {
@@ -468,20 +476,20 @@ void check_ir(const char*                   mdparin,
                         "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to "
                         "%s.",
                         enumValueToString(ir->etc),
-                        ei_names[ir->eI]);
+                        enumValueToString(ir->eI));
             }
             warning_note(wi, warn_buf);
         }
         ir->etc = TemperatureCoupling::No;
     }
-    if (ir->eI == eiVVAK)
+    if (ir->eI == IntegrationAlgorithm::VVAK)
     {
         sprintf(warn_buf,
                 "Integrator method %s is implemented primarily for validation purposes; for "
                 "molecular dynamics, you should probably be using %s or %s",
-                ei_names[eiVVAK],
-                ei_names[eiMD],
-                ei_names[eiVV]);
+                enumValueToString(IntegrationAlgorithm::VVAK),
+                enumValueToString(IntegrationAlgorithm::MD),
+                enumValueToString(IntegrationAlgorithm::VV));
         warning_note(wi, warn_buf);
     }
     if (!EI_DYNAMICS(ir->eI))
@@ -491,7 +499,7 @@ void check_ir(const char*                   mdparin,
             sprintf(warn_buf,
                     "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.",
                     enumValueToString(ir->epc),
-                    ei_names[ir->eI]);
+                    enumValueToString(ir->eI));
             warning_note(wi, warn_buf);
         }
         ir->epc = PressureCoupling::No;
@@ -517,7 +525,7 @@ void check_ir(const char*                   mdparin,
             }
         }
         else if ((ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy)
-                 || (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
+                 || (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl > 0
                      && (ir->nstcalcenergy > ir->fepvals->nstdhdl)))
 
         {
@@ -527,7 +535,7 @@ void check_ir(const char*                   mdparin,
             int         min_nst  = ir->nstenergy;
 
             /* find the smallest of ( nstenergy, nstdhdl ) */
-            if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
+            if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl > 0
                 && (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
             {
                 min_nst  = ir->fepvals->nstdhdl;
@@ -555,7 +563,7 @@ void check_ir(const char*                   mdparin,
 
         if (ir->nstcalcenergy > 0)
         {
-            if (ir->efep != efepNO)
+            if (ir->efep != FreeEnergyPerturbationType::No)
             {
                 /* nstdhdl should be a multiple of nstcalcenergy */
                 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstdhdl", &ir->fepvals->nstdhdl, wi);
@@ -592,7 +600,7 @@ void check_ir(const char*                   mdparin,
     }
 
     /* LD STUFF */
-    if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
+    if ((EI_SD(ir->eI) || ir->eI == IntegrationAlgorithm::BD) && ir->bContinuation && ir->ld_seed != -1)
     {
         warning_note(wi,
                      "You are doing a continuation with SD or BD, make sure that ld_seed is "
@@ -617,7 +625,7 @@ void check_ir(const char*                   mdparin,
         warning(wi, warn_buf);
     }
 
-    if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
+    if ((EI_SD(ir->eI) || ir->eI == IntegrationAlgorithm::BD) && ir->bContinuation && ir->ld_seed != -1)
     {
         warning_note(wi,
                      "You are doing a continuation with SD or BD, make sure that ld_seed is "
@@ -633,10 +641,12 @@ void check_ir(const char*                   mdparin,
             sprintf(err_buf,
                     "Entry %d for %s must be between 0 and 1, instead is %g",
                     i,
-                    efpt_names[efptTEMPERATURE],
-                    fep->all_lambda[efptTEMPERATURE][i]);
-            CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
-            if (fep->all_lambda[efptTEMPERATURE][i] > 0)
+                    enumValueToString(FreeEnergyPerturbationCouplingType::Temperature),
+                    fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i]);
+            CHECK((fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i] < 0)
+                  || (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i]
+                      > 1));
+            if (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i] > 0)
             {
                 bAllTempZero = FALSE;
             }
@@ -645,7 +655,7 @@ void check_ir(const char*                   mdparin,
         CHECK(bAllTempZero == TRUE);
 
         sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
-        CHECK(ir->eI != eiVV);
+        CHECK(ir->eI != IntegrationAlgorithm::VV);
 
         /* check compatability of the temperature coupling with simulated tempering */
 
@@ -680,9 +690,9 @@ void check_ir(const char*                   mdparin,
 
     /* verify free energy options */
 
-    if (ir->efep != efepNO)
+    if (ir->efep != FreeEnergyPerturbationType::No)
     {
-        fep = ir->fepvals;
+        fep = ir->fepvals.get();
         sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power);
         CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
 
@@ -701,13 +711,13 @@ void check_ir(const char*                   mdparin,
         sprintf(err_buf,
                 "Can't use positive delta-lambda (%g) with expanded ensemble simulations",
                 fep->delta_lambda);
-        CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
+        CHECK(fep->delta_lambda > 0 && (ir->efep == FreeEnergyPerturbationType::Expanded));
 
         sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
-        CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
+        CHECK(!(EI_VV(ir->eI)) && (ir->efep == FreeEnergyPerturbationType::Expanded));
 
         sprintf(err_buf, "Free-energy not implemented for Ewald");
-        CHECK(ir->coulombtype == eelEWALD);
+        CHECK(ir->coulombtype == CoulombInteractionType::Ewald);
 
         /* check validty of lambda inputs */
         if (fep->n_lambda == 0)
@@ -741,7 +751,7 @@ void check_ir(const char*                   mdparin,
         {
             int n_lambda_terms;
             n_lambda_terms = 0;
-            for (i = 0; i < efptNR; i++)
+            for (i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
             {
                 if (fep->separate_dvdl[i])
                 {
@@ -765,14 +775,15 @@ void check_ir(const char*                   mdparin,
             }
         }
 
-        for (j = 0; j < efptNR; j++)
+        for (j = 0; j < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); j++)
         {
             for (i = 0; i < fep->n_lambda; i++)
             {
+                auto enumValue = static_cast<FreeEnergyPerturbationCouplingType>(j);
                 sprintf(err_buf,
                         "Entry %d for %s must be between 0 and 1, instead is %g",
                         i,
-                        efpt_names[j],
+                        enumValueToString(enumValue),
                         fep->all_lambda[j][i]);
                 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
             }
@@ -787,11 +798,14 @@ void check_ir(const char*                   mdparin,
                         "coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to "
                         "crashes, and is not supported.",
                         i,
-                        fep->all_lambda[efptVDW][i],
-                        fep->all_lambda[efptCOUL][i]);
+                        fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)][i],
+                        fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)][i]);
                 CHECK((fep->sc_alpha > 0)
-                      && (((fep->all_lambda[efptCOUL][i] > 0.0) && (fep->all_lambda[efptCOUL][i] < 1.0))
-                          && ((fep->all_lambda[efptVDW][i] > 0.0) && (fep->all_lambda[efptVDW][i] < 1.0))));
+                      && (((fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)][i] > 0.0)
+                           && (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)][i] < 1.0))
+                          && ((fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)][i] > 0.0)
+                              && (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)][i]
+                                  < 1.0))));
             }
         }
 
@@ -821,19 +835,20 @@ void check_ir(const char*                   mdparin,
         /*  Free Energy Checks -- In an ideal world, slow growth and FEP would
             be treated differently, but that's the next step */
 
-        for (i = 0; i < efptNR; i++)
+        for (i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
         {
+            auto enumValue = static_cast<FreeEnergyPerturbationCouplingType>(i);
             for (j = 0; j < fep->n_lambda; j++)
             {
-                sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
+                sprintf(err_buf, "%s[%d] must be between 0 and 1", enumValueToString(enumValue), j);
                 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
             }
         }
     }
 
-    if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
+    if ((ir->bSimTemp) || (ir->efep == FreeEnergyPerturbationType::Expanded))
     {
-        fep = ir->fepvals;
+        fep = ir->fepvals.get();
 
         /* checking equilibration of weights inputs for validity */
 
@@ -841,72 +856,75 @@ void check_ir(const char*                   mdparin,
                 "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal "
                 "to %s",
                 expand->equil_n_at_lam,
-                elmceq_names[elmceqNUMATLAM]);
-        CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
+                enumValueToString(LambdaWeightWillReachEquilibrium::NumAtLambda));
+        CHECK((expand->equil_n_at_lam > 0)
+              && (expand->elmceq != LambdaWeightWillReachEquilibrium::NumAtLambda));
 
         sprintf(err_buf,
                 "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to "
                 "%s",
                 expand->equil_samples,
-                elmceq_names[elmceqSAMPLES]);
-        CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
+                enumValueToString(LambdaWeightWillReachEquilibrium::Samples));
+        CHECK((expand->equil_samples > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::Samples));
 
         sprintf(err_buf,
                 "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
                 expand->equil_steps,
-                elmceq_names[elmceqSTEPS]);
-        CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
+                enumValueToString(LambdaWeightWillReachEquilibrium::Steps));
+        CHECK((expand->equil_steps > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::Steps));
 
         sprintf(err_buf,
                 "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
                 expand->equil_samples,
-                elmceq_names[elmceqWLDELTA]);
-        CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
+                enumValueToString(LambdaWeightWillReachEquilibrium::WLDelta));
+        CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::WLDelta));
 
         sprintf(err_buf,
                 "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
                 expand->equil_ratio,
-                elmceq_names[elmceqRATIO]);
-        CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
+                enumValueToString(LambdaWeightWillReachEquilibrium::Ratio));
+        CHECK((expand->equil_ratio > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::Ratio));
 
         sprintf(err_buf,
                 "weight-equil-number-all-lambda (%d) must be a positive integer if "
                 "lmc-weights-equil=%s",
                 expand->equil_n_at_lam,
-                elmceq_names[elmceqNUMATLAM]);
-        CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
+                enumValueToString(LambdaWeightWillReachEquilibrium::NumAtLambda));
+        CHECK((expand->equil_n_at_lam <= 0)
+              && (expand->elmceq == LambdaWeightWillReachEquilibrium::NumAtLambda));
 
         sprintf(err_buf,
                 "weight-equil-number-samples (%d) must be a positive integer if "
                 "lmc-weights-equil=%s",
                 expand->equil_samples,
-                elmceq_names[elmceqSAMPLES]);
-        CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
+                enumValueToString(LambdaWeightWillReachEquilibrium::Samples));
+        CHECK((expand->equil_samples <= 0) && (expand->elmceq == LambdaWeightWillReachEquilibrium::Samples));
 
         sprintf(err_buf,
                 "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
                 expand->equil_steps,
-                elmceq_names[elmceqSTEPS]);
-        CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
+                enumValueToString(LambdaWeightWillReachEquilibrium::Steps));
+        CHECK((expand->equil_steps <= 0) && (expand->elmceq == LambdaWeightWillReachEquilibrium::Steps));
 
         sprintf(err_buf,
                 "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
                 expand->equil_wl_delta,
-                elmceq_names[elmceqWLDELTA]);
-        CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
+                enumValueToString(LambdaWeightWillReachEquilibrium::WLDelta));
+        CHECK((expand->equil_wl_delta <= 0)
+              && (expand->elmceq == LambdaWeightWillReachEquilibrium::WLDelta));
 
         sprintf(err_buf,
                 "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
                 expand->equil_ratio,
-                elmceq_names[elmceqRATIO]);
-        CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
+                enumValueToString(LambdaWeightWillReachEquilibrium::Ratio));
+        CHECK((expand->equil_ratio <= 0) && (expand->elmceq == LambdaWeightWillReachEquilibrium::Ratio));
 
         sprintf(err_buf,
                 "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
-                elmceq_names[elmceqWLDELTA],
-                elamstats_names[elamstatsWL],
-                elamstats_names[elamstatsWWL]);
-        CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
+                enumValueToString(LambdaWeightWillReachEquilibrium::WLDelta),
+                enumValueToString(LambdaWeightCalculation::WL),
+                enumValueToString(LambdaWeightCalculation::WWL));
+        CHECK((expand->elmceq == LambdaWeightWillReachEquilibrium::WLDelta) && (!EWL(expand->elamstats)));
 
         sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
         CHECK((expand->lmc_repeats <= 0));
@@ -920,7 +938,7 @@ void check_ir(const char*                   mdparin,
                 fep->init_fep_state,
                 expand->lmc_forced_nstart);
         CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0)
-              && (expand->elmcmove != elmcmoveNO));
+              && (expand->elmcmove != LambdaMoveCalculation::No));
         sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
         CHECK((expand->lmc_forced_nstart < 0));
         sprintf(err_buf,
@@ -936,8 +954,8 @@ void check_ir(const char*                   mdparin,
         CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
 
         /* if there is no temperature control, we need to specify an MC temperature */
-        if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO)
-            && (expand->mc_temp <= 0.0))
+        if (!integratorHasReferenceTemperature(ir)
+            && (expand->elmcmove != LambdaMoveCalculation::No) && (expand->mc_temp <= 0.0))
         {
             sprintf(err_buf,
                     "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set "
@@ -988,7 +1006,7 @@ void check_ir(const char*                   mdparin,
         sprintf(err_buf,
                 "Can not have dispersion correction with pbc=%s",
                 c_pbcTypeNames[ir->pbcType].c_str());
-        CHECK(ir->eDispCorr != edispcNO);
+        CHECK(ir->eDispCorr != DispersionCorrectionType::No);
     }
 
     if (ir->rlist == 0.0)
@@ -998,10 +1016,11 @@ void check_ir(const char*                   mdparin,
                 "with coulombtype = %s or coulombtype = %s\n"
                 "without periodic boundary conditions (pbc = %s) and\n"
                 "rcoulomb and rvdw set to zero",
-                eel_names[eelCUT],
-                eel_names[eelUSER],
+                enumValueToString(CoulombInteractionType::Cut),
+                enumValueToString(CoulombInteractionType::User),
                 c_pbcTypeNames[PbcType::No].c_str());
-        CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER))
+        CHECK(((ir->coulombtype != CoulombInteractionType::Cut)
+               && (ir->coulombtype != CoulombInteractionType::User))
               || (ir->pbcType != PbcType::No) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
 
         if (ir->nstlist > 0)
@@ -1017,9 +1036,9 @@ void check_ir(const char*                   mdparin,
     {
         // TODO Change this behaviour. There should be exactly one way
         // to turn off an algorithm.
-        ir->comm_mode = ecmNO;
+        ir->comm_mode = ComRemovalAlgorithm::No;
     }
-    if (ir->comm_mode != ecmNO)
+    if (ir->comm_mode != ComRemovalAlgorithm::No)
     {
         if (ir->nstcomm < 0)
         {
@@ -1042,7 +1061,7 @@ void check_ir(const char*                   mdparin,
             ir->nstcomm = ir->nstcalcenergy;
         }
 
-        if (ir->comm_mode == ecmANGULAR)
+        if (ir->comm_mode == ComRemovalAlgorithm::Angular)
         {
             sprintf(err_buf,
                     "Can not remove the rotation around the center of mass with periodic "
@@ -1058,13 +1077,14 @@ void check_ir(const char*                   mdparin,
         }
     }
 
-    if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->pbcType == PbcType::No && ir->comm_mode != ecmANGULAR)
+    if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->pbcType == PbcType::No
+        && ir->comm_mode != ComRemovalAlgorithm::Angular)
     {
         sprintf(warn_buf,
                 "Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
                 "in a non-periodic system. You should probably set comm_mode = ANGULAR or use "
                 "integrator = %s.",
-                ei_names[eiSD1]);
+                enumValueToString(IntegrationAlgorithm::SD1));
         warning_note(wi, warn_buf);
     }
 
@@ -1102,12 +1122,12 @@ void check_ir(const char*                   mdparin,
         ir->opts.nhchainlength = 0;
     }
 
-    if (ir->eI == eiVVAK)
+    if (ir->eI == IntegrationAlgorithm::VVAK)
     {
         sprintf(err_buf,
                 "%s implemented primarily for validation, and requires nsttcouple = 1 and "
                 "nstpcouple = 1.",
-                ei_names[eiVVAK]);
+                enumValueToString(IntegrationAlgorithm::VVAK));
         CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
     }
 
@@ -1116,7 +1136,7 @@ void check_ir(const char*                   mdparin,
         sprintf(err_buf,
                 "%s temperature control not supported for integrator %s.",
                 enumValueToString(ir->etc),
-                ei_names[ir->eI]);
+                enumValueToString(ir->eI));
         CHECK(!(EI_VV(ir->eI)));
 
         if (ir->nstcomm > 0 && (ir->etc == TemperatureCoupling::Andersen))
@@ -1218,13 +1238,13 @@ void check_ir(const char*                   mdparin,
     /* ELECTROSTATICS */
     /* More checks are in triple check (grompp.c) */
 
-    if (ir->coulombtype == eelSWITCH)
+    if (ir->coulombtype == CoulombInteractionType::Switch)
     {
         sprintf(warn_buf,
                 "coulombtype = %s is only for testing purposes and can lead to serious "
                 "artifacts, advice: use coulombtype = %s",
-                eel_names[ir->coulombtype],
-                eel_names[eelRF_ZERO]);
+                enumValueToString(ir->coulombtype),
+                enumValueToString(CoulombInteractionType::RFZero));
         warning(wi, warn_buf);
     }
 
@@ -1259,11 +1279,11 @@ void check_ir(const char*                   mdparin,
     {
         /* reaction field (at the cut-off) */
 
-        if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
+        if (ir->coulombtype == CoulombInteractionType::RFZero && ir->epsilon_rf != 0)
         {
             sprintf(warn_buf,
                     "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
-                    eel_names[ir->coulombtype]);
+                    enumValueToString(ir->coulombtype));
             warning(wi, warn_buf);
             ir->epsilon_rf = 0.0;
         }
@@ -1274,7 +1294,7 @@ void check_ir(const char*                   mdparin,
         {
             sprintf(warn_buf,
                     "Using epsilon-rf = epsilon-r with %s does not make sense",
-                    eel_names[ir->coulombtype]);
+                    enumValueToString(ir->coulombtype));
             warning(wi, warn_buf);
         }
     }
@@ -1289,28 +1309,28 @@ void check_ir(const char*                   mdparin,
             sprintf(err_buf,
                     "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the "
                     "potential modifier options!",
-                    eel_names[ir->coulombtype]);
+                    enumValueToString(ir->coulombtype));
             CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
         }
     }
 
-    if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
+    if (ir->coulombtype == CoulombInteractionType::Switch || ir->coulombtype == CoulombInteractionType::Shift)
     {
         sprintf(err_buf,
                 "Explicit switch/shift coulomb interactions cannot be used in combination with a "
                 "secondary coulomb-modifier.");
-        CHECK(ir->coulomb_modifier != eintmodNONE);
+        CHECK(ir->coulomb_modifier != InteractionModifiers::None);
     }
-    if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
+    if (ir->vdwtype == VanDerWaalsType::Switch || ir->vdwtype == VanDerWaalsType::Shift)
     {
         sprintf(err_buf,
                 "Explicit switch/shift vdw interactions cannot be used in combination with a "
                 "secondary vdw-modifier.");
-        CHECK(ir->vdw_modifier != eintmodNONE);
+        CHECK(ir->vdw_modifier != InteractionModifiers::None);
     }
 
-    if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT || ir->vdwtype == evdwSWITCH
-        || ir->vdwtype == evdwSHIFT)
+    if (ir->coulombtype == CoulombInteractionType::Switch || ir->coulombtype == CoulombInteractionType::Shift
+        || ir->vdwtype == VanDerWaalsType::Switch || ir->vdwtype == VanDerWaalsType::Shift)
     {
         sprintf(warn_buf,
                 "The switch/shift interaction settings are just for compatibility; you will get "
@@ -1319,7 +1339,8 @@ void check_ir(const char*                   mdparin,
         warning_note(wi, warn_buf);
     }
 
-    if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
+    if (ir->coulombtype == CoulombInteractionType::PmeSwitch
+        || ir->coulomb_modifier == InteractionModifiers::PotSwitch)
     {
         if (ir->rcoulomb_switch / ir->rcoulomb < 0.9499)
         {
@@ -1336,7 +1357,7 @@ void check_ir(const char*                   mdparin,
         }
     }
 
-    if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
+    if (ir->vdwtype == VanDerWaalsType::Switch || ir->vdw_modifier == InteractionModifiers::PotSwitch)
     {
         if (ir->rvdw_switch == 0)
         {
@@ -1351,10 +1372,13 @@ void check_ir(const char*                   mdparin,
 
     if (EEL_FULL(ir->coulombtype))
     {
-        if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER
-            || ir->coulombtype == eelPMEUSERSWITCH)
+        if (ir->coulombtype == CoulombInteractionType::PmeSwitch
+            || ir->coulombtype == CoulombInteractionType::PmeUser
+            || ir->coulombtype == CoulombInteractionType::PmeUserSwitch)
         {
-            sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist", eel_names[ir->coulombtype]);
+            sprintf(err_buf,
+                    "With coulombtype = %s, rcoulomb must be <= rlist",
+                    enumValueToString(ir->coulombtype));
             CHECK(ir->rcoulomb > ir->rlist);
         }
     }
@@ -1363,13 +1387,13 @@ void check_ir(const char*                   mdparin,
     {
         // TODO: Move these checks into the ewald module with the options class
         int orderMin = 3;
-        int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
+        int orderMax = (ir->coulombtype == CoulombInteractionType::P3mAD ? 8 : 12);
 
         if (ir->pme_order < orderMin || ir->pme_order > orderMax)
         {
             sprintf(warn_buf,
                     "With coulombtype = %s, you should have %d <= pme-order <= %d",
-                    eel_names[ir->coulombtype],
+                    enumValueToString(ir->coulombtype),
                     orderMin,
                     orderMax);
             warning_error(wi, warn_buf);
@@ -1378,24 +1402,25 @@ void check_ir(const char*                   mdparin,
 
     if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
     {
-        if (ir->ewald_geometry == eewg3D)
+        if (ir->ewald_geometry == EwaldGeometry::ThreeD)
         {
             sprintf(warn_buf,
                     "With pbc=%s you should use ewald-geometry=%s",
                     c_pbcTypeNames[ir->pbcType].c_str(),
-                    eewg_names[eewg3DC]);
+                    enumValueToString(EwaldGeometry::ThreeDC));
             warning(wi, warn_buf);
         }
         /* This check avoids extra pbc coding for exclusion corrections */
         sprintf(err_buf, "wall-ewald-zfac should be >= 2");
         CHECK(ir->wall_ewald_zfac < 2);
     }
-    if ((ir->ewald_geometry == eewg3DC) && (ir->pbcType != PbcType::XY) && EEL_FULL(ir->coulombtype))
+    if ((ir->ewald_geometry == EwaldGeometry::ThreeDC) && (ir->pbcType != PbcType::XY)
+        && EEL_FULL(ir->coulombtype))
     {
         sprintf(warn_buf,
                 "With %s and ewald_geometry = %s you should use pbc = %s",
-                eel_names[ir->coulombtype],
-                eewg_names[eewg3DC],
+                enumValueToString(ir->coulombtype),
+                enumValueToString(EwaldGeometry::ThreeDC),
                 c_pbcTypeNames[PbcType::XY].c_str());
         warning(wi, warn_buf);
     }
@@ -1429,20 +1454,21 @@ void check_ir(const char*                   mdparin,
         }
     }
 
-    if (ir->vdwtype == evdwPME)
+    if (ir->vdwtype == VanDerWaalsType::Pme)
     {
-        if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
+        if (!(ir->vdw_modifier == InteractionModifiers::None
+              || ir->vdw_modifier == InteractionModifiers::PotShift))
         {
             sprintf(err_buf,
                     "With vdwtype = %s, the only supported modifiers are %s and %s",
-                    evdw_names[ir->vdwtype],
-                    eintmod_names[eintmodPOTSHIFT],
-                    eintmod_names[eintmodNONE]);
+                    enumValueToString(ir->vdwtype),
+                    enumValueToString(InteractionModifiers::PotShift),
+                    enumValueToString(InteractionModifiers::None));
             warning_error(wi, err_buf);
         }
     }
 
-    if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
+    if (ir->vdwtype == VanDerWaalsType::User && ir->eDispCorr != DispersionCorrectionType::No)
     {
         warning_note(wi,
                      "You have selected user tables with dispersion correction, the dispersion "
@@ -1451,20 +1477,22 @@ void check_ir(const char*                   mdparin,
                      "really want dispersion correction to -C6/r^6.");
     }
 
-    if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT) && ir->rvdw != 0)
+    if (ir->eI == IntegrationAlgorithm::LBFGS
+        && (ir->coulombtype == CoulombInteractionType::Cut || ir->vdwtype == VanDerWaalsType::Cut)
+        && ir->rvdw != 0)
     {
         warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
     }
 
-    if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
+    if (ir->eI == IntegrationAlgorithm::LBFGS && ir->nbfgscorr <= 0)
     {
         warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
     }
 
     /* IMPLICIT SOLVENT */
-    if (ir->coulombtype == eelGB_NOTUSED)
+    if (ir->coulombtype == CoulombInteractionType::GBNotused)
     {
-        sprintf(warn_buf, "Invalid option %s for coulombtype", eel_names[ir->coulombtype]);
+        sprintf(warn_buf, "Invalid option %s for coulombtype", enumValueToString(ir->coulombtype));
         warning_error(wi, warn_buf);
     }
 
@@ -1479,7 +1507,7 @@ void check_ir(const char*                   mdparin,
     }
 
     // cosine acceleration is only supported in leap-frog
-    if (ir->cos_accel != 0.0 && ir->eI != eiMD)
+    if (ir->cos_accel != 0.0 && ir->eI != IntegrationAlgorithm::MD)
     {
         warning_error(wi, "cos-acceleration is only supported by integrator = md");
     }
@@ -1490,17 +1518,17 @@ void check_ir(const char*                   mdparin,
    str = the input string
    n = the (pre-allocated) number of doubles read
    r = the output array of doubles. */
-static void parse_n_real(char* str, int* n, real** r, warninp_t wi)
+static std::vector<real> parse_n_real(const std::string& str, int* n, warninp_t wi)
 {
     auto values = gmx::splitString(str);
     *n          = values.size();
 
-    snew(*r, *n);
+    std::vector<real> r;
     for (int i = 0; i < *n; i++)
     {
         try
         {
-            (*r)[i] = gmx::fromString<real>(values[i]);
+            r.emplace_back(gmx::fromString<real>(values[i]));
         }
         catch (gmx::GromacsException&)
         {
@@ -1509,33 +1537,33 @@ static void parse_n_real(char* str, int* n, real** r, warninp_t wi)
                                   + " in string in mdp file. Expected a real number.");
         }
     }
+    return r;
 }
 
 
-static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
+static void do_fep_params(t_inputrec* ir, gmx::ArrayRef<std::string> fep_lambda, char weights[STRLEN], warninp_t wi)
 {
 
-    int         i, j, max_n_lambda, nweights, nfep[efptNR];
-    t_lambda*   fep    = ir->fepvals;
-    t_expanded* expand = ir->expandedvals;
-    real**      count_fep_lambdas;
-    bool        bOneLambda = TRUE;
-
-    snew(count_fep_lambdas, efptNR);
+    int         i, j, max_n_lambda, nweights;
+    t_lambda*   fep    = ir->fepvals.get();
+    t_expanded* expand = ir->expandedvals.get();
+    gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, std::vector<real>> count_fep_lambdas;
+    bool                                                                         bOneLambda = TRUE;
+    gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, int>               nfep;
 
     /* FEP input processing */
     /* first, identify the number of lambda values for each type.
        All that are nonzero must have the same number */
 
-    for (i = 0; i < efptNR; i++)
+    for (auto i : keysOf(nfep))
     {
-        parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
+        count_fep_lambdas[i] = parse_n_real(fep_lambda[static_cast<int>(i)], &(nfep[i]), wi);
     }
 
     /* now, determine the number of components.  All must be either zero, or equal. */
 
     max_n_lambda = 0;
-    for (i = 0; i < efptNR; i++)
+    for (auto i : keysOf(nfep))
     {
         if (nfep[i] > max_n_lambda)
         {
@@ -1545,7 +1573,7 @@ static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weight
         }
     }
 
-    for (i = 0; i < efptNR; i++)
+    for (auto i : keysOf(nfep))
     {
         if (nfep[i] == 0)
         {
@@ -1553,7 +1581,7 @@ static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weight
         }
         else if (nfep[i] == max_n_lambda)
         {
-            if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute
+            if (i != FreeEnergyPerturbationCouplingType::Temperature) /* we treat this differently -- not really a reason to compute
                                          the derivative with respect to the temperature currently */
             {
                 ir->fepvals->separate_dvdl[i] = TRUE;
@@ -1565,28 +1593,26 @@ static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weight
                       "Number of lambdas (%d) for FEP type %s not equal to number of other types "
                       "(%d)",
                       nfep[i],
-                      efpt_names[i],
+                      enumValueToString(i),
                       max_n_lambda);
         }
     }
     /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
-    ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
+    ir->fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Temperature] = FALSE;
 
     /* the number of lambdas is the number we've read in, which is either zero
        or the same for all */
     fep->n_lambda = max_n_lambda;
 
-    /* allocate space for the array of lambda values */
-    snew(fep->all_lambda, efptNR);
     /* if init_lambda is defined, we need to set lambda */
     if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
     {
-        ir->fepvals->separate_dvdl[efptFEP] = TRUE;
+        ir->fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Fep] = TRUE;
     }
     /* otherwise allocate the space for all of the lambdas, and transfer the data */
-    for (i = 0; i < efptNR; i++)
+    for (auto i : keysOf(nfep))
     {
-        snew(fep->all_lambda[i], fep->n_lambda);
+        fep->all_lambda[i].resize(fep->n_lambda);
         if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
                             are zero */
         {
@@ -1594,19 +1620,17 @@ static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weight
             {
                 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
             }
-            sfree(count_fep_lambdas[i]);
         }
     }
-    sfree(count_fep_lambdas);
 
     /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for
        internal bookkeeping -- for now, init_lambda */
 
-    if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
+    if ((nfep[FreeEnergyPerturbationCouplingType::Fep] == 0) && (fep->init_lambda >= 0))
     {
         for (i = 0; i < fep->n_lambda; i++)
         {
-            fep->all_lambda[efptFEP][i] = fep->init_lambda;
+            fep->all_lambda[FreeEnergyPerturbationCouplingType::Fep][i] = fep->init_lambda;
         }
     }
 
@@ -1619,9 +1643,9 @@ static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weight
     }
     else
     {
-        for (i = 0; i < efptNR; i++)
+        for (auto i : keysOf(nfep))
         {
-            if ((nfep[i] != 0) && (i != efptFEP))
+            if ((nfep[i] != 0) && (i != FreeEnergyPerturbationCouplingType::Fep))
             {
                 bOneLambda = FALSE;
             }
@@ -1636,23 +1660,23 @@ static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weight
        specified (i.e. nfep[i] == 0).  This means if fep is not defined,
        they are all zero. */
 
-    for (i = 0; i < efptNR; i++)
+    for (auto i : keysOf(nfep))
     {
-        if ((nfep[i] == 0) && (i != efptFEP))
+        if ((nfep[i] == 0) && (i != FreeEnergyPerturbationCouplingType::Fep))
         {
             for (j = 0; j < fep->n_lambda; j++)
             {
-                fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
+                fep->all_lambda[i][j] = fep->all_lambda[FreeEnergyPerturbationCouplingType::Fep][j];
             }
         }
     }
 
 
     /* now read in the weights */
-    parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
+    expand->init_lambda_weights = parse_n_real(weights, &nweights, wi);
     if (nweights == 0)
     {
-        snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
+        expand->init_lambda_weights.resize(fep->n_lambda); /* initialize to zero */
     }
     else if (nweights != fep->n_lambda)
     {
@@ -1661,7 +1685,7 @@ static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weight
                   nweights,
                   fep->n_lambda);
     }
-    if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
+    if ((expand->nstexpanded < 0) && (ir->efep != FreeEnergyPerturbationType::No))
     {
         expand->nstexpanded = fep->nstdhdl;
         /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
@@ -1673,7 +1697,9 @@ static void do_simtemp_params(t_inputrec* ir)
 {
 
     snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
-    GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
+    getSimTemps(ir->fepvals->n_lambda,
+                ir->simtempvals.get(),
+                ir->fepvals->all_lambda[FreeEnergyPerturbationCouplingType::Temperature]);
 }
 
 template<typename T>
@@ -1747,7 +1773,7 @@ static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_densi
             opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
         }
 
-        if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
+        if (ir->wall_type == WallType::NineThree || ir->wall_type == WallType::TenFour)
         {
             auto wallDensity = gmx::splitString(wall_density);
             if (wallDensity.size() != size_t(ir->nwall))
@@ -1786,10 +1812,10 @@ static void read_expandedparams(std::vector<t_inpfile>* inp, t_expanded* expand,
 {
     /* read expanded ensemble parameters */
     printStringNewline(inp, "expanded ensemble variables");
-    expand->nstexpanded    = get_eint(inp, "nstexpanded", -1, wi);
-    expand->elamstats      = get_eeenum(inp, "lmc-stats", elamstats_names, wi);
-    expand->elmcmove       = get_eeenum(inp, "lmc-move", elmcmove_names, wi);
-    expand->elmceq         = get_eeenum(inp, "lmc-weights-equil", elmceq_names, wi);
+    expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
+    expand->elamstats   = getEnum<LambdaWeightCalculation>(inp, "lmc-stats", wi);
+    expand->elmcmove    = getEnum<LambdaMoveCalculation>(inp, "lmc-move", wi);
+    expand->elmceq      = getEnum<LambdaWeightWillReachEquilibrium>(inp, "lmc-weights-equil", wi);
     expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
     expand->equil_samples  = get_eint(inp, "weight-equil-number-samples", -1, wi);
     expand->equil_steps    = get_eint(inp, "weight-equil-number-steps", -1, wi);
@@ -1802,61 +1828,14 @@ static void read_expandedparams(std::vector<t_inpfile>* inp, t_expanded* expand,
     expand->gibbsdeltalam     = get_eint(inp, "lmc-gibbsdelta", -1, wi);
     expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
     expand->bSymmetrizedTMatrix =
-            (get_eeenum(inp, "symmetrized-transition-matrix", yesno_names, wi) != 0);
+            (getEnum<Boolean>(inp, "symmetrized-transition-matrix", wi) != Boolean::No);
     expand->nstTij        = get_eint(inp, "nst-transition-matrix", -1, wi);
     expand->minvarmin     = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
     expand->c_range       = get_eint(inp, "weight-c-range", 0, wi);    /* default is just C=0 */
     expand->wl_scale      = get_ereal(inp, "wl-scale", 0.8, wi);
     expand->wl_ratio      = get_ereal(inp, "wl-ratio", 0.8, wi);
     expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
-    expand->bWLoneovert   = (get_eeenum(inp, "wl-oneovert", yesno_names, wi) != 0);
-}
-
-template<typename EnumType>
-EnumType getEnum(std::vector<t_inpfile>* inp, const char* name, warninp_t wi)
-{
-    // If we there's no valid option, we'll use the first enum entry as default.
-    // Note, this assumes the enum is zero based, which is also assumed by
-    // EnumerationWrapper and EnumerationArray.
-    const auto  defaultEnumValue = EnumType::Default;
-    const auto& defaultName      = enumValueToString(defaultEnumValue);
-    // Get index of option in input
-    const auto ii = get_einp(inp, name);
-    if (ii == -1)
-    {
-        // If the option wasn't set, we use the first enum entry as default
-        inp->back().value_.assign(defaultName);
-        return defaultEnumValue;
-    }
-
-    // Check if option string can be mapped to a valid enum value
-    const auto* optionString = (*inp)[ii].value_.c_str();
-    for (auto enumValue : gmx::EnumerationWrapper<EnumType>{})
-    {
-        if (gmx_strcasecmp_min(enumValueToString(enumValue), optionString) == 0)
-        {
-            return enumValue;
-        }
-    }
-
-    // If we get here, the option set is invalid. Print error.
-    std::string errorMessage = gmx::formatString(
-            "Invalid enum '%s' for variable %s, using '%s'\n", optionString, name, defaultName);
-    errorMessage += gmx::formatString("Next time, use one of:");
-    for (auto enumValue : gmx::EnumerationWrapper<EnumType>{})
-    {
-        errorMessage += gmx::formatString(" '%s'", enumValueToString(enumValue));
-    }
-    if (wi != nullptr)
-    {
-        warning_error(wi, errorMessage);
-    }
-    else
-    {
-        fprintf(stderr, "%s\n", errorMessage.c_str());
-    }
-    (*inp)[ii].value_.assign(defaultName);
-    return defaultEnumValue;
+    expand->bWLoneovert   = (getEnum<Boolean>(inp, "wl-oneovert", wi) != Boolean::No);
 }
 
 /*! \brief Return whether an end state with the given coupling-lambda
@@ -1920,8 +1899,8 @@ void get_ir(const char*     mdparin,
     double      dumdub[2][6];
     int         i, j, m;
     char        warn_buf[STRLEN];
-    t_lambda*   fep    = ir->fepvals;
-    t_expanded* expand = ir->expandedvals;
+    t_lambda*   fep    = ir->fepvals.get();
+    t_expanded* expand = ir->expandedvals.get();
 
     const char* no_names[] = { "no", nullptr };
 
@@ -1988,7 +1967,7 @@ void get_ir(const char*     mdparin,
     setStringEntry(&inp, "define", opts->define, nullptr);
 
     printStringNewline(&inp, "RUN CONTROL PARAMETERS");
-    ir->eI = get_eeenum(&inp, "integrator", ei_names, wi);
+    ir->eI = getEnum<IntegrationAlgorithm>(&inp, "integrator", wi);
     printStringNoNewline(&inp, "Start time and timestep in ps");
     ir->init_t  = get_ereal(&inp, "tinit", 0.0, wi);
     ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
@@ -1999,7 +1978,7 @@ void get_ir(const char*     mdparin,
             &inp, "Part index is updated automatically on checkpointing (keeps files separate)");
     ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
     printStringNoNewline(&inp, "Multiple time-stepping");
-    ir->useMts = (get_eeenum(&inp, "mts", yesno_names, wi) != 0);
+    ir->useMts = (getEnum<Boolean>(&inp, "mts", wi) != Boolean::No);
     if (ir->useMts)
     {
         gmx::GromppMtsOpts& mtsOpts = opts->mtsOpts;
@@ -2014,7 +1993,7 @@ void get_ir(const char*     mdparin,
         }
     }
     printStringNoNewline(&inp, "mode for center of mass motion removal");
-    ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
+    ir->comm_mode = getEnum<ComRemovalAlgorithm>(&inp, "comm-mode", wi);
     printStringNoNewline(&inp, "number of steps for center of mass motion removal");
     ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
     printStringNoNewline(&inp, "group(s) for center of mass motion removal");
@@ -2064,7 +2043,7 @@ void get_ir(const char*     mdparin,
     /* Neighbor searching */
     printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
     printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs)");
-    ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
+    ir->cutoff_scheme = getEnum<CutoffScheme>(&inp, "cutoff-scheme", wi);
     printStringNoNewline(&inp, "nblist update frequency");
     ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
     printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
@@ -2075,7 +2054,7 @@ void get_ir(const char*     mdparin,
         pbcTypesNamesChar.push_back(pbcTypeName.c_str());
     }
     ir->pbcType       = static_cast<PbcType>(get_eeenum(&inp, "pbc", pbcTypesNamesChar.data(), wi));
-    ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
+    ir->bPeriodicMols = getEnum<Boolean>(&inp, "periodic-molecules", wi) != Boolean::No;
     printStringNoNewline(&inp,
                          "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
     printStringNoNewline(&inp, "a value of -1 means: use rlist");
@@ -2087,8 +2066,8 @@ void get_ir(const char*     mdparin,
     /* Electrostatics */
     printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
     printStringNoNewline(&inp, "Method for doing electrostatics");
-    ir->coulombtype      = get_eeenum(&inp, "coulombtype", eel_names, wi);
-    ir->coulomb_modifier = get_eeenum(&inp, "coulomb-modifier", eintmod_names, wi);
+    ir->coulombtype      = getEnum<CoulombInteractionType>(&inp, "coulombtype", wi);
+    ir->coulomb_modifier = getEnum<InteractionModifiers>(&inp, "coulomb-modifier", wi);
     printStringNoNewline(&inp, "cut-off lengths");
     ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
     ir->rcoulomb        = get_ereal(&inp, "rcoulomb", 1.0, wi);
@@ -2096,13 +2075,13 @@ void get_ir(const char*     mdparin,
     ir->epsilon_r  = get_ereal(&inp, "epsilon-r", 1.0, wi);
     ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
     printStringNoNewline(&inp, "Method for doing Van der Waals");
-    ir->vdwtype      = get_eeenum(&inp, "vdw-type", evdw_names, wi);
-    ir->vdw_modifier = get_eeenum(&inp, "vdw-modifier", eintmod_names, wi);
+    ir->vdwtype      = getEnum<VanDerWaalsType>(&inp, "vdw-type", wi);
+    ir->vdw_modifier = getEnum<InteractionModifiers>(&inp, "vdw-modifier", wi);
     printStringNoNewline(&inp, "cut-off lengths");
     ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
     ir->rvdw        = get_ereal(&inp, "rvdw", 1.0, wi);
     printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
-    ir->eDispCorr = get_eeenum(&inp, "DispCorr", edispc_names, wi);
+    ir->eDispCorr = getEnum<DispersionCorrectionType>(&inp, "DispCorr", wi);
     printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
     ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
     printStringNoNewline(&inp, "Separate tables between energy group pairs");
@@ -2117,8 +2096,8 @@ void get_ir(const char*     mdparin,
     ir->pme_order              = get_eint(&inp, "pme-order", 4, wi);
     ir->ewald_rtol             = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
     ir->ewald_rtol_lj          = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
-    ir->ljpme_combination_rule = get_eeenum(&inp, "lj-pme-comb-rule", eljpme_names, wi);
-    ir->ewald_geometry         = get_eeenum(&inp, "ewald-geometry", eewg_names, wi);
+    ir->ljpme_combination_rule = getEnum<LongRangeVdW>(&inp, "lj-pme-comb-rule", wi);
+    ir->ewald_geometry         = getEnum<EwaldGeometry>(&inp, "ewald-geometry", wi);
     ir->epsilon_surface        = get_ereal(&inp, "epsilon-surface", 0.0, wi);
 
     /* Implicit solvation is no longer supported, but we need grompp
@@ -2132,7 +2111,7 @@ void get_ir(const char*     mdparin,
     ir->etc                = getEnum<TemperatureCoupling>(&inp, "tcoupl", wi);
     ir->nsttcouple         = get_eint(&inp, "nsttcouple", -1, wi);
     ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
-    ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
+    ir->bPrintNHChains = (getEnum<Boolean>(&inp, "print-nose-hoover-chain-variables", wi) != Boolean::No);
     printStringNoNewline(&inp, "Groups to couple separately");
     setStringEntry(&inp, "tc-grps", inputrecStrings->tcgrps, nullptr);
     printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
@@ -2140,18 +2119,18 @@ void get_ir(const char*     mdparin,
     setStringEntry(&inp, "ref-t", inputrecStrings->ref_t, nullptr);
     printStringNoNewline(&inp, "pressure coupling");
     ir->epc        = getEnum<PressureCoupling>(&inp, "pcoupl", wi);
-    ir->epct       = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
+    ir->epct       = getEnum<PressureCouplingType>(&inp, "pcoupltype", wi);
     ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
     printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
     ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
     setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
     setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
     printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
-    ir->refcoord_scaling = get_eeenum(&inp, "refcoord-scaling", erefscaling_names, wi);
+    ir->refcoord_scaling = getEnum<RefCoordScaling>(&inp, "refcoord-scaling", wi);
 
     /* QMMM */
     printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
-    ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
+    ir->bQMMM = (getEnum<Boolean>(&inp, "QMMM", wi) != Boolean::No);
     printStringNoNewline(&inp, "Groups treated with MiMiC");
     setStringEntry(&inp, "QMMM-grps", inputrecStrings->QMMM, nullptr);
 
@@ -2169,7 +2148,7 @@ void get_ir(const char*     mdparin,
 
     /* Startup run */
     printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
-    opts->bGenVel = (get_eeenum(&inp, "gen-vel", yesno_names, wi) != 0);
+    opts->bGenVel = (getEnum<Boolean>(&inp, "gen-vel", wi) != Boolean::No);
     opts->tempi   = get_ereal(&inp, "gen-temp", 300.0, wi);
     opts->seed    = get_eint(&inp, "gen-seed", -1, wi);
 
@@ -2177,12 +2156,12 @@ void get_ir(const char*     mdparin,
     printStringNewline(&inp, "OPTIONS FOR BONDS");
     opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
     printStringNoNewline(&inp, "Type of constraint algorithm");
-    ir->eConstrAlg = get_eeenum(&inp, "constraint-algorithm", econstr_names, wi);
+    ir->eConstrAlg = getEnum<ConstraintAlgorithm>(&inp, "constraint-algorithm", wi);
     printStringNoNewline(&inp, "Do not constrain the start configuration");
-    ir->bContinuation = (get_eeenum(&inp, "continuation", yesno_names, wi) != 0);
+    ir->bContinuation = (getEnum<Boolean>(&inp, "continuation", wi) != Boolean::No);
     printStringNoNewline(&inp,
                          "Use successive overrelaxation to reduce the number of shake iterations");
-    ir->bShakeSOR = (get_eeenum(&inp, "Shake-SOR", yesno_names, wi) != 0);
+    ir->bShakeSOR = (getEnum<Boolean>(&inp, "Shake-SOR", wi) != Boolean::No);
     printStringNoNewline(&inp, "Relative tolerance of shake");
     ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
     printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
@@ -2195,7 +2174,7 @@ void get_ir(const char*     mdparin,
     printStringNoNewline(&inp, "rotates over more degrees than");
     ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
     printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
-    opts->bMorse = (get_eeenum(&inp, "morse", yesno_names, wi) != 0);
+    opts->bMorse = (getEnum<Boolean>(&inp, "morse", wi) != Boolean::No);
 
     /* Energy group exclusions */
     printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
@@ -2208,7 +2187,7 @@ void get_ir(const char*     mdparin,
     printStringNoNewline(
             &inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
     ir->nwall         = get_eint(&inp, "nwall", 0, wi);
-    ir->wall_type     = get_eeenum(&inp, "wall-type", ewt_names, wi);
+    ir->wall_type     = getEnum<WallType>(&inp, "wall-type", wi);
     ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
     setStringEntry(&inp, "wall-atomtype", inputrecStrings->wall_atomtype, nullptr);
     setStringEntry(&inp, "wall-density", inputrecStrings->wall_density, nullptr);
@@ -2216,7 +2195,7 @@ void get_ir(const char*     mdparin,
 
     /* COM pulling */
     printStringNewline(&inp, "COM PULLING");
-    ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
+    ir->bPull = (getEnum<Boolean>(&inp, "pull", wi) != Boolean::No);
     if (ir->bPull)
     {
         ir->pull                        = std::make_unique<pull_params_t>();
@@ -2226,7 +2205,7 @@ void get_ir(const char*     mdparin,
         {
             for (int c = 0; c < ir->pull->ncoord; c++)
             {
-                if (ir->pull->coord[c].eType == epullCONSTRAINT)
+                if (ir->pull->coord[c].eType == PullingAlgorithm::Constraint)
                 {
                     warning_error(wi,
                                   "Constraint COM pulling is not supported in combination with "
@@ -2240,7 +2219,7 @@ void get_ir(const char*     mdparin,
     /* AWH biasing
        NOTE: needs COM pulling or free energy input */
     printStringNewline(&inp, "AWH biasing");
-    ir->bDoAwh = (get_eeenum(&inp, "awh", yesno_names, wi) != 0);
+    ir->bDoAwh = (getEnum<Boolean>(&inp, "awh", wi) != Boolean::No);
     if (ir->bDoAwh)
     {
         ir->awhParams = gmx::readAwhParams(&inp, wi);
@@ -2249,7 +2228,7 @@ void get_ir(const char*     mdparin,
     /* Enforced rotation */
     printStringNewline(&inp, "ENFORCED ROTATION");
     printStringNoNewline(&inp, "Enforced rotation: No or Yes");
-    ir->bRot = (get_eeenum(&inp, "rotation", yesno_names, wi) != 0);
+    ir->bRot = (getEnum<Boolean>(&inp, "rotation", wi) != Boolean::No);
     if (ir->bRot)
     {
         snew(ir->rot, 1);
@@ -2269,18 +2248,18 @@ void get_ir(const char*     mdparin,
     /* Refinement */
     printStringNewline(&inp, "NMR refinement stuff");
     printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
-    ir->eDisre = get_eeenum(&inp, "disre", edisre_names, wi);
+    ir->eDisre = getEnum<DistanceRestraintRefinement>(&inp, "disre", wi);
     printStringNoNewline(
             &inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
-    ir->eDisreWeighting = get_eeenum(&inp, "disre-weighting", edisreweighting_names, wi);
+    ir->eDisreWeighting = getEnum<DistanceRestraintWeighting>(&inp, "disre-weighting", wi);
     printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
-    ir->bDisreMixed = (get_eeenum(&inp, "disre-mixed", yesno_names, wi) != 0);
+    ir->bDisreMixed = (getEnum<Boolean>(&inp, "disre-mixed", wi) != Boolean::No);
     ir->dr_fc       = get_ereal(&inp, "disre-fc", 1000.0, wi);
     ir->dr_tau      = get_ereal(&inp, "disre-tau", 0.0, wi);
     printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
     ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
     printStringNoNewline(&inp, "Orientation restraints: No or Yes");
-    opts->bOrire = (get_eeenum(&inp, "orire", yesno_names, wi) != 0);
+    opts->bOrire = (getEnum<Boolean>(&inp, "orire", wi) != Boolean::No);
     printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
     ir->orires_fc  = get_ereal(&inp, "orire-fc", 0.0, wi);
     ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
@@ -2290,11 +2269,11 @@ void get_ir(const char*     mdparin,
 
     /* free energy variables */
     printStringNewline(&inp, "Free energy variables");
-    ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
+    ir->efep = getEnum<FreeEnergyPerturbationType>(&inp, "free-energy", wi);
     setStringEntry(&inp, "couple-moltype", inputrecStrings->couple_moltype, nullptr);
     opts->couple_lam0  = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
     opts->couple_lam1  = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
-    opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
+    opts->bCoupleIntra = (getEnum<Boolean>(&inp, "couple-intramol", wi) != Boolean::No);
 
     fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
                                                                          we can recognize if
@@ -2302,25 +2281,32 @@ void get_ir(const char*     mdparin,
     fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
     fep->delta_lambda   = get_ereal(&inp, "delta-lambda", 0.0, wi);
     fep->nstdhdl        = get_eint(&inp, "nstdhdl", 50, wi);
-    setStringEntry(&inp, "fep-lambdas", inputrecStrings->fep_lambda[efptFEP], nullptr);
-    setStringEntry(&inp, "mass-lambdas", inputrecStrings->fep_lambda[efptMASS], nullptr);
-    setStringEntry(&inp, "coul-lambdas", inputrecStrings->fep_lambda[efptCOUL], nullptr);
-    setStringEntry(&inp, "vdw-lambdas", inputrecStrings->fep_lambda[efptVDW], nullptr);
-    setStringEntry(&inp, "bonded-lambdas", inputrecStrings->fep_lambda[efptBONDED], nullptr);
-    setStringEntry(&inp, "restraint-lambdas", inputrecStrings->fep_lambda[efptRESTRAINT], nullptr);
-    setStringEntry(&inp, "temperature-lambdas", inputrecStrings->fep_lambda[efptTEMPERATURE], nullptr);
+    inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Fep] =
+            setStringEntry(&inp, "fep-lambdas", "");
+    inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Mass] =
+            setStringEntry(&inp, "mass-lambdas", "");
+    inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Coul] =
+            setStringEntry(&inp, "coul-lambdas", "");
+    inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Vdw] =
+            setStringEntry(&inp, "vdw-lambdas", "");
+    inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Bonded] =
+            setStringEntry(&inp, "bonded-lambdas", "");
+    inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Restraint] =
+            setStringEntry(&inp, "restraint-lambdas", "");
+    inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Temperature] =
+            setStringEntry(&inp, "temperature-lambdas", "");
     fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
     setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr);
-    fep->edHdLPrintEnergy   = get_eeenum(&inp, "dhdl-print-energy", edHdLPrintEnergy_names, wi);
+    fep->edHdLPrintEnergy   = getEnum<FreeEnergyPrintEnergy>(&inp, "dhdl-print-energy", wi);
     fep->sc_alpha           = get_ereal(&inp, "sc-alpha", 0.0, wi);
     fep->sc_power           = get_eint(&inp, "sc-power", 1, wi);
     fep->sc_r_power         = get_ereal(&inp, "sc-r-power", 6.0, wi);
     fep->sc_sigma           = get_ereal(&inp, "sc-sigma", 0.3, wi);
-    fep->bScCoul            = (get_eeenum(&inp, "sc-coul", yesno_names, wi) != 0);
+    fep->bScCoul            = (getEnum<Boolean>(&inp, "sc-coul", wi) != Boolean::No);
     fep->dh_hist_size       = get_eint(&inp, "dh_hist_size", 0, wi);
     fep->dh_hist_spacing    = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
-    fep->separate_dhdl_file = get_eeenum(&inp, "separate-dhdl-file", separate_dhdl_file_names, wi);
-    fep->dhdl_derivatives   = get_eeenum(&inp, "dhdl-derivatives", dhdl_derivatives_names, wi);
+    fep->separate_dhdl_file = getEnum<SeparateDhdlFile>(&inp, "separate-dhdl-file", wi);
+    fep->dhdl_derivatives   = getEnum<DhDlDerivativeCalculation>(&inp, "dhdl-derivatives", wi);
     fep->dh_hist_size       = get_eint(&inp, "dh_hist_size", 0, wi);
     fep->dh_hist_spacing    = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
 
@@ -2333,13 +2319,13 @@ void get_ir(const char*     mdparin,
 
     /* simulated tempering variables */
     printStringNewline(&inp, "simulated tempering variables");
-    ir->bSimTemp = (get_eeenum(&inp, "simulated-tempering", yesno_names, wi) != 0);
-    ir->simtempvals->eSimTempScale = get_eeenum(&inp, "simulated-tempering-scaling", esimtemp_names, wi);
+    ir->bSimTemp = (getEnum<Boolean>(&inp, "simulated-tempering", wi) != Boolean::No);
+    ir->simtempvals->eSimTempScale = getEnum<SimulatedTempering>(&inp, "simulated-tempering-scaling", wi);
     ir->simtempvals->simtemp_low  = get_ereal(&inp, "sim-temp-low", 300.0, wi);
     ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
 
     /* expanded ensemble variables */
-    if (ir->efep == efepEXPANDED || ir->bSimTemp)
+    if (ir->efep == FreeEnergyPerturbationType::Expanded || ir->bSimTemp)
     {
         read_expandedparams(&inp, expand, wi);
     }
@@ -2367,8 +2353,8 @@ void get_ir(const char*     mdparin,
     printStringNewline(&inp,
                        "Ion/water position swapping for computational electrophysiology setups");
     printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
-    ir->eSwapCoords = get_eeenum(&inp, "swapcoords", eSwapTypes_names, wi);
-    if (ir->eSwapCoords != eswapNO)
+    ir->eSwapCoords = getEnum<SwapType>(&inp, "swapcoords", wi);
+    if (ir->eSwapCoords != SwapType::No)
     {
         char buf[STRLEN];
         int  nIonTypes;
@@ -2383,7 +2369,7 @@ void get_ir(const char*     mdparin,
         {
             warning_error(wi, "You need to provide at least one ion type for position exchanges.");
         }
-        ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
+        ir->swap->ngrp = nIonTypes + static_cast<int>(SwapGroupSplittingType::Count);
         snew(ir->swap->grp, ir->swap->ngrp);
         for (i = 0; i < ir->swap->ngrp; i++)
         {
@@ -2391,16 +2377,25 @@ void get_ir(const char*     mdparin,
         }
         printStringNoNewline(&inp,
                              "Two index groups that contain the compartment-partitioning atoms");
-        setStringEntry(&inp, "split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
-        setStringEntry(&inp, "split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
+        setStringEntry(&inp,
+                       "split-group0",
+                       ir->swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname,
+                       nullptr);
+        setStringEntry(&inp,
+                       "split-group1",
+                       ir->swap->grp[static_cast<int>(SwapGroupSplittingType::Split1)].molname,
+                       nullptr);
         printStringNoNewline(&inp,
                              "Use center of mass of split groups (yes/no), otherwise center of "
                              "geometry is used");
-        ir->swap->massw_split[0] = (get_eeenum(&inp, "massw-split0", yesno_names, wi) != 0);
-        ir->swap->massw_split[1] = (get_eeenum(&inp, "massw-split1", yesno_names, wi) != 0);
+        ir->swap->massw_split[0] = (getEnum<Boolean>(&inp, "massw-split0", wi) != Boolean::No);
+        ir->swap->massw_split[1] = (getEnum<Boolean>(&inp, "massw-split1", wi) != Boolean::No);
 
         printStringNoNewline(&inp, "Name of solvent molecules");
-        setStringEntry(&inp, "solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
+        setStringEntry(&inp,
+                       "solvent-group",
+                       ir->swap->grp[static_cast<int>(SwapGroupSplittingType::Solvent)].molname,
+                       nullptr);
 
         printStringNoNewline(&inp,
                              "Split cylinder: radius, upper and lower extension (nm) (this will "
@@ -2430,7 +2425,7 @@ void get_ir(const char*     mdparin,
         printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
         for (i = 0; i < nIonTypes; i++)
         {
-            int ig = eSwapFixedGrpNR + i;
+            int ig = static_cast<int>(SwapGroupSplittingType::Count) + i;
 
             sprintf(buf, "iontype%d-name", i);
             setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
@@ -2508,7 +2503,7 @@ void get_ir(const char*     mdparin,
         {
             switch (ir->epct)
             {
-                case epctISOTROPIC:
+                case PressureCouplingType::Isotropic:
                     if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
                     {
                         warning_error(
@@ -2517,8 +2512,8 @@ void get_ir(const char*     mdparin,
                     }
                     dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
                     break;
-                case epctSEMIISOTROPIC:
-                case epctSURFACETENSION:
+                case PressureCouplingType::SemiIsotropic:
+                case PressureCouplingType::SurfaceTension:
                     if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
                     {
                         warning_error(
@@ -2527,7 +2522,7 @@ void get_ir(const char*     mdparin,
                     }
                     dumdub[m][YY] = dumdub[m][XX];
                     break;
-                case epctANISOTROPIC:
+                case PressureCouplingType::Anisotropic:
                     if (sscanf(dumstr[m],
                                "%lf%lf%lf%lf%lf%lf",
                                &(dumdub[m][XX]),
@@ -2546,7 +2541,7 @@ void get_ir(const char*     mdparin,
                 default:
                     gmx_fatal(FARGS,
                               "Pressure coupling type %s not implemented yet",
-                              epcoupltype_names[ir->epct]);
+                              enumValueToString(ir->epct));
             }
         }
     }
@@ -2557,7 +2552,7 @@ void get_ir(const char*     mdparin,
         ir->ref_p[i][i]    = dumdub[1][i];
         ir->compress[i][i] = dumdub[0][i];
     }
-    if (ir->epct == epctANISOTROPIC)
+    if (ir->epct == PressureCouplingType::Anisotropic)
     {
         ir->ref_p[XX][YY] = dumdub[1][3];
         ir->ref_p[XX][ZZ] = dumdub[1][4];
@@ -2581,7 +2576,7 @@ void get_ir(const char*     mdparin,
         }
     }
 
-    if (ir->comm_mode == ecmNO)
+    if (ir->comm_mode == ComRemovalAlgorithm::No)
     {
         ir->nstcomm = 0;
     }
@@ -2589,14 +2584,15 @@ void get_ir(const char*     mdparin,
     opts->couple_moltype = nullptr;
     if (strlen(inputrecStrings->couple_moltype) > 0)
     {
-        if (ir->efep != efepNO)
+        if (ir->efep != FreeEnergyPerturbationType::No)
         {
             opts->couple_moltype = gmx_strdup(inputrecStrings->couple_moltype);
             if (opts->couple_lam0 == opts->couple_lam1)
             {
                 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
             }
-            if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam1 == ecouplamNONE))
+            if (ir->eI == IntegrationAlgorithm::MD
+                && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam1 == ecouplamNONE))
             {
                 warning_note(
                         wi,
@@ -2612,23 +2608,23 @@ void get_ir(const char*     mdparin,
         }
     }
     /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
-    if (ir->efep != efepNO)
+    if (ir->efep != FreeEnergyPerturbationType::No)
     {
         if (fep->delta_lambda != 0)
         {
-            ir->efep = efepSLOWGROWTH;
+            ir->efep = FreeEnergyPerturbationType::SlowGrowth;
         }
     }
 
-    if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
+    if (fep->edHdLPrintEnergy == FreeEnergyPrintEnergy::Yes)
     {
-        fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
+        fep->edHdLPrintEnergy = FreeEnergyPrintEnergy::Total;
         warning_note(wi,
                      "Old option for dhdl-print-energy given: "
                      "changing \"yes\" to \"total\"\n");
     }
 
-    if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
+    if (ir->bSimTemp && (fep->edHdLPrintEnergy == FreeEnergyPrintEnergy::No))
     {
         /* always print out the energy to dhdl if we are doing
            expanded ensemble, since we need the total energy for
@@ -2637,13 +2633,13 @@ void get_ir(const char*     mdparin,
            we will allow that if the appropriate mdp setting has
            been enabled. Otherwise, total it is:
          */
-        fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
+        fep->edHdLPrintEnergy = FreeEnergyPrintEnergy::Total;
     }
 
-    if ((ir->efep != efepNO) || ir->bSimTemp)
+    if ((ir->efep != FreeEnergyPerturbationType::No) || ir->bSimTemp)
     {
         ir->bExpanded = FALSE;
-        if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
+        if ((ir->efep == FreeEnergyPerturbationType::Expanded) || ir->bSimTemp)
         {
             ir->bExpanded = TRUE;
         }
@@ -2661,7 +2657,8 @@ void get_ir(const char*     mdparin,
          * If the (advanced) user does FEP through manual topology changes,
          * this check will not be triggered.
          */
-        if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 && ir->fepvals->sc_alpha != 0
+        if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->n_lambda == 0
+            && ir->fepvals->sc_alpha != 0
             && (couple_lambda_has_vdw_on(opts->couple_lam0) && couple_lambda_has_vdw_on(opts->couple_lam1)))
         {
             warning(wi,
@@ -2756,7 +2753,7 @@ void get_ir(const char*     mdparin,
     }
 
     /* Ion/water position swapping checks */
-    if (ir->eSwapCoords != eswapNO)
+    if (ir->eSwapCoords != SwapType::No)
     {
         if (ir->swap->nstswap < 1)
         {
@@ -3111,7 +3108,7 @@ static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
 
         for (int i = 0; i < pull->ncoord; i++)
         {
-            if (pull->coord[i].eType != epullCONSTRAINT)
+            if (pull->coord[i].eType != PullingAlgorithm::Constraint)
             {
                 continue;
             }
@@ -3169,8 +3166,8 @@ static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
         {
             switch (ir->comm_mode)
             {
-                case ecmLINEAR:
-                case ecmLINEAR_ACCELERATION_CORRECTION:
+                case ComRemovalAlgorithm::Linear:
+                case ComRemovalAlgorithm::LinearAccelerationCorrection:
                     nrdf_vcm_sub[j] = 0;
                     for (int d = 0; d < ndim_rm_vcm; d++)
                     {
@@ -3180,7 +3177,7 @@ static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
                         }
                     }
                     break;
-                case ecmANGULAR: nrdf_vcm_sub[j] = 6; break;
+                case ComRemovalAlgorithm::Angular: nrdf_vcm_sub[j] = 6; break;
                 default: gmx_incons("Checking comm_mode");
             }
         }
@@ -3306,9 +3303,13 @@ static void make_swap_groups(t_swapcoords* swap, t_blocka* grps, char** gnames)
 
 
     /* Just a quick check here, more thorough checks are in mdrun */
-    if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
+    if (strcmp(swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname,
+               swap->grp[static_cast<int>(SwapGroupSplittingType::Split1)].molname)
+        == 0)
     {
-        gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
+        gmx_fatal(FARGS,
+                  "The split groups can not both be '%s'.",
+                  swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname);
     }
 
     /* Get the index atoms of the split0, split1, solvent, and swap groups */
@@ -3322,7 +3323,7 @@ static void make_swap_groups(t_swapcoords* swap, t_blocka* grps, char** gnames)
         {
             fprintf(stderr,
                     "%s group '%s' contains %d atoms.\n",
-                    ig < 3 ? eSwapFixedGrp_names[ig] : "Swap",
+                    ig < 3 ? enumValueToString(static_cast<SwapGroupSplittingType>(ig)) : "Swap",
                     swap->grp[ig].molname,
                     swapg->nat);
             snew(swapg->ind, swapg->nat);
@@ -3530,7 +3531,7 @@ void do_index(const char*                   mdparin,
     snew(ir->opts.nrdf, nr);
     snew(ir->opts.tau_t, nr);
     snew(ir->opts.ref_t, nr);
-    if (ir->eI == eiBD && ir->bd_fric == 0)
+    if (ir->eI == IntegrationAlgorithm::BD && ir->bd_fric == 0)
     {
         fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
     }
@@ -3546,9 +3547,11 @@ void do_index(const char*                   mdparin,
         convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
         for (i = 0; (i < nr); i++)
         {
-            if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
+            if ((ir->eI == IntegrationAlgorithm::BD) && ir->opts.tau_t[i] <= 0)
             {
-                sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
+                sprintf(warn_buf,
+                        "With integrator %s tau-t should be larger than 0",
+                        enumValueToString(ir->eI));
                 warning_error(wi, warn_buf);
             }
 
@@ -3671,7 +3674,7 @@ void do_index(const char*                   mdparin,
         snew(ir->opts.anneal_temp, nr);
         for (i = 0; i < nr; i++)
         {
-            ir->opts.annealing[i]      = eannNO;
+            ir->opts.annealing[i]      = SimulatedAnnealing::No;
             ir->opts.anneal_npoints[i] = 0;
             ir->opts.anneal_time[i]    = nullptr;
             ir->opts.anneal_temp[i]    = nullptr;
@@ -3683,16 +3686,16 @@ void do_index(const char*                   mdparin,
             {
                 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
                 {
-                    ir->opts.annealing[i] = eannNO;
+                    ir->opts.annealing[i] = SimulatedAnnealing::No;
                 }
                 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
                 {
-                    ir->opts.annealing[i] = eannSINGLE;
+                    ir->opts.annealing[i] = SimulatedAnnealing::Single;
                     bAnneal               = TRUE;
                 }
                 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
                 {
-                    ir->opts.annealing[i] = eannPERIODIC;
+                    ir->opts.annealing[i] = SimulatedAnnealing::Periodic;
                     bAnneal               = TRUE;
                 }
             }
@@ -3784,13 +3787,13 @@ void do_index(const char*                   mdparin,
                 /* Print out some summary information, to make sure we got it right */
                 for (i = 0; i < nr; i++)
                 {
-                    if (ir->opts.annealing[i] != eannNO)
+                    if (ir->opts.annealing[i] != SimulatedAnnealing::No)
                     {
                         j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
                         fprintf(stderr,
                                 "Simulated annealing for group %s: %s, %d timepoints\n",
                                 *(groups->groupNames[j]),
-                                eann_names[ir->opts.annealing[i]],
+                                enumValueToString(ir->opts.annealing[i]),
                                 ir->opts.anneal_npoints[i]);
                         fprintf(stderr, "Time (ps)   Temperature (K)\n");
                         /* All terms except the last one */
@@ -3804,7 +3807,7 @@ void do_index(const char*                   mdparin,
 
                         /* Finally the last one */
                         j = ir->opts.anneal_npoints[i] - 1;
-                        if (ir->opts.annealing[i] == eannSINGLE)
+                        if (ir->opts.annealing[i] == SimulatedAnnealing::Single)
                         {
                             fprintf(stderr,
                                     "%9.1f-     %5.1f\n",
@@ -3850,7 +3853,7 @@ void do_index(const char*                   mdparin,
         make_rotation_groups(ir->rot, inputrecStrings->rotateGroupNames, defaultIndexGroups, gnames);
     }
 
-    if (ir->eSwapCoords != eswapNO)
+    if (ir->eSwapCoords != SwapType::No)
     {
         make_swap_groups(ir->swap, defaultIndexGroups, gnames);
     }
@@ -3938,7 +3941,7 @@ void do_index(const char*                   mdparin,
                  bVerbose,
                  wi);
 
-    if (ir->comm_mode != ecmNO)
+    if (ir->comm_mode != ComRemovalAlgorithm::No)
     {
         checkAndUpdateVcmFreezeGroupConsistency(groups, natoms, ir->opts, wi);
     }
@@ -4029,7 +4032,7 @@ void do_index(const char*                   mdparin,
     snew(ir->opts.egp_flags, nr * nr);
 
     bExcl = do_egp_flag(ir, groups, "energygrp-excl", inputrecStrings->egpexcl, EGP_EXCL);
-    if (bExcl && ir->cutoff_scheme == ecutsVERLET)
+    if (bExcl && ir->cutoff_scheme == CutoffScheme::Verlet)
     {
         warning_error(wi, "Energy group exclusions are currently not supported");
     }
@@ -4039,8 +4042,10 @@ void do_index(const char*                   mdparin,
     }
 
     bTable = do_egp_flag(ir, groups, "energygrp-table", inputrecStrings->egptable, EGP_TABLE);
-    if (bTable && !(ir->vdwtype == evdwUSER) && !(ir->coulombtype == eelUSER)
-        && !(ir->coulombtype == eelPMEUSER) && !(ir->coulombtype == eelPMEUSERSWITCH))
+    if (bTable && !(ir->vdwtype == VanDerWaalsType::User)
+        && !(ir->coulombtype == CoulombInteractionType::User)
+        && !(ir->coulombtype == CoulombInteractionType::PmeUser)
+        && !(ir->coulombtype == CoulombInteractionType::PmeUserSwitch))
     {
         gmx_fatal(FARGS,
                   "Can only have energy group pair tables in combination with user tables for VdW "
@@ -4273,7 +4278,7 @@ static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t* mtop
 
     check_combination_rule_differences(
             mtop, 0, &bC6ParametersWorkWithGeometricRules, &bC6ParametersWorkWithLBRules, &bLBRulesPossible);
-    if (ir->ljpme_combination_rule == eljpmeLB)
+    if (ir->ljpme_combination_rule == LongRangeVdW::LB)
     {
         if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
         {
@@ -4287,7 +4292,7 @@ static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t* mtop
     {
         if (!bC6ParametersWorkWithGeometricRules)
         {
-            if (ir->eDispCorr != edispcNO)
+            if (ir->eDispCorr != DispersionCorrectionType::No)
             {
                 warning_note(wi,
                              "You are using geometric combination rules in "
@@ -4335,7 +4340,7 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_
                      "macro-molecule, the artifacts are usually negligible.");
     }
 
-    if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0 && ir->nstlist > 1
+    if (ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol > 0 && ir->nstlist > 1
         && ((EI_MD(ir->eI) || EI_SD(ir->eI))
             && (ir->etc == TemperatureCoupling::VRescale || ir->etc == TemperatureCoupling::Berendsen)))
     {
@@ -4402,7 +4407,7 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_
             CHECK(ir->opts.tau_t[i] < 0);
         }
 
-        if (ir->etc == TemperatureCoupling::AndersenMassive && ir->comm_mode != ecmNO)
+        if (ir->etc == TemperatureCoupling::AndersenMassive && ir->comm_mode != ComRemovalAlgorithm::No)
         {
             for (i = 0; i < ir->opts.ngtc; i++)
             {
@@ -4422,7 +4427,8 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_
         }
     }
 
-    if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD && ir->comm_mode == ecmNO
+    if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != IntegrationAlgorithm::BD
+        && ir->comm_mode == ComRemovalAlgorithm::No
         && !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) && !ETC_ANDERSEN(ir->etc))
     {
         warning(wi,
@@ -4453,7 +4459,7 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_
     }
 
     /* Check for pressure coupling with absolute position restraints */
-    if (ir->epc != PressureCoupling::No && ir->refcoord_scaling == erscNO)
+    if (ir->epc != PressureCoupling::No && ir->refcoord_scaling == RefCoordScaling::No)
     {
         absolute_reference(ir, sys, TRUE, AbsRef);
         {
@@ -4489,19 +4495,19 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_
                     "You are using full electrostatics treatment %s for a system without charges.\n"
                     "This costs a lot of performance for just processing zeros, consider using %s "
                     "instead.\n",
-                    EELTYPE(ir->coulombtype),
-                    EELTYPE(eelCUT));
+                    enumValueToString(ir->coulombtype),
+                    enumValueToString(CoulombInteractionType::Cut));
             warning(wi, err_buf);
         }
     }
     else
     {
-        if (ir->coulombtype == eelCUT && ir->rcoulomb > 0)
+        if (ir->coulombtype == CoulombInteractionType::Cut && ir->rcoulomb > 0)
         {
             sprintf(err_buf,
                     "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
                     "You might want to consider using %s electrostatics.\n",
-                    EELTYPE(eelPME));
+                    enumValueToString(CoulombInteractionType::Pme));
             warning_note(wi, err_buf);
         }
     }
@@ -4513,7 +4519,7 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_
     }
 
     /* Generalized reaction field */
-    if (ir->coulombtype == eelGRF_NOTUSED)
+    if (ir->coulombtype == CoulombInteractionType::GRFNotused)
     {
         warning_error(wi,
                       "Generalized reaction-field electrostatics is no longer supported. "
@@ -4521,7 +4527,7 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_
                       "constant by hand.");
     }
 
-    if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0
+    if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->sc_alpha != 0
         && !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
     {
         gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
@@ -4560,12 +4566,13 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_
                 {
                     for (c = 0; c < ir->pull->ncoord; c++)
                     {
-                        if (ir->pull->coord[c].eGeom == epullgDIRPBC && ir->pull->coord[c].vec[m] != 0)
+                        if (ir->pull->coord[c].eGeom == PullGroupGeometry::DirectionPBC
+                            && ir->pull->coord[c].vec[m] != 0)
                         {
                             gmx_fatal(FARGS,
                                       "Can not have dynamic box while using pull geometry '%s' "
                                       "(dim %c)",
-                                      EPULLGEOM(ir->pull->coord[c].eGeom),
+                                      enumValueToString(ir->pull->coord[c].eGeom),
                                       'x' + m);
                         }
                     }
@@ -4588,7 +4595,7 @@ void double_check(t_inputrec* ir, matrix box, bool bHasNormalConstraints, bool b
         warning_error(wi, ptr);
     }
 
-    if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
+    if (bHasNormalConstraints && ir->eConstrAlg == ConstraintAlgorithm::Shake)
     {
         if (ir->shake_tol <= 0.0)
         {
@@ -4597,22 +4604,23 @@ void double_check(t_inputrec* ir, matrix box, bool bHasNormalConstraints, bool b
         }
     }
 
-    if ((ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
+    if ((ir->eConstrAlg == ConstraintAlgorithm::Lincs) && bHasNormalConstraints)
     {
         /* If we have Lincs constraints: */
-        if (ir->eI == eiMD && ir->etc == TemperatureCoupling::No && ir->eConstrAlg == econtLINCS
-            && ir->nLincsIter == 1)
+        if (ir->eI == IntegrationAlgorithm::MD && ir->etc == TemperatureCoupling::No
+            && ir->eConstrAlg == ConstraintAlgorithm::Lincs && ir->nLincsIter == 1)
         {
             sprintf(warn_buf,
                     "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
             warning_note(wi, warn_buf);
         }
 
-        if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
+        if ((ir->eI == IntegrationAlgorithm::CG || ir->eI == IntegrationAlgorithm::LBFGS)
+            && (ir->nProjOrder < 8))
         {
             sprintf(warn_buf,
                     "For accurate %s with LINCS constraints, lincs-order should be 8 or more.",
-                    ei_names[ir->eI]);
+                    enumValueToString(ir->eI));
             warning_note(wi, warn_buf);
         }
         if (ir->epc == PressureCoupling::Mttk)