Fix random typos
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / grompp.cpp
index cb052c723a681d73c0c4226888d9ac1c680928a5..899b5008c2becadbe5008ba31f7652865f5c7fa8 100644 (file)
@@ -4,7 +4,7 @@
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -39,6 +39,7 @@
 
 #include "grompp.h"
 
+#include <array>
 #include <cerrno>
 #include <climits>
 #include <cmath>
@@ -79,6 +80,7 @@
 #include "gromacs/mdlib/calc_verletbuf.h"
 #include "gromacs/mdlib/compute_io.h"
 #include "gromacs/mdlib/constr.h"
+#include "gromacs/mdlib/md_support.h"
 #include "gromacs/mdlib/perf_est.h"
 #include "gromacs/mdlib/vsite.h"
 #include "gromacs/mdrun/mdmodules.h"
 #include "gromacs/utility/listoflists.h"
 #include "gromacs/utility/logger.h"
 #include "gromacs/utility/loggerbuilder.h"
-#include "gromacs/utility/mdmodulenotification.h"
+#include "gromacs/utility/mdmodulesnotifiers.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/snprintf.h"
 
 InteractionOfType::InteractionOfType(gmx::ArrayRef<const int>  atoms,
                                      gmx::ArrayRef<const real> params,
                                      const std::string&        name) :
-    atoms_(atoms.begin(), atoms.end()),
-    interactionTypeName_(name)
+    atoms_(atoms.begin(), atoms.end()), interactionTypeName_(name)
 {
     GMX_RELEASE_ASSERT(
             params.size() <= forceParam_.size(),
             gmx::formatString("Cannot have more parameters than the maximum number possible (%d)", MAXFORCEPARAM)
                     .c_str());
-    auto forceParamIt = forceParam_.begin();
+    std::array<real, MAXFORCEPARAM>::iterator forceParamIt = forceParam_.begin();
     for (const auto param : params)
     {
         *forceParamIt++ = param;
@@ -467,7 +468,8 @@ static void check_vel(gmx_mtop_t* mtop, rvec v[])
     {
         const t_atom& local = atomP.atom();
         int           i     = atomP.globalAtomNumber();
-        if (local.ptype == eptShell || local.ptype == eptBond || local.ptype == eptVSite)
+        if (local.ptype == ParticleType::Shell || local.ptype == ParticleType::Bond
+            || local.ptype == ParticleType::VSite)
         {
             clear_rvec(v[i]);
         }
@@ -481,7 +483,7 @@ static void check_shells_inputrec(gmx_mtop_t* mtop, t_inputrec* ir, warninp* wi)
     for (const AtomProxy atomP : AtomRange(*mtop))
     {
         const t_atom& local = atomP.atom();
-        if (local.ptype == eptShell || local.ptype == eptBond)
+        if (local.ptype == ParticleType::Shell || local.ptype == ParticleType::Bond)
         {
             nshells++;
         }
@@ -585,7 +587,7 @@ static void new_status(const char*                           topfile,
                        std::vector<MoleculeInformation>*     mi,
                        std::unique_ptr<MoleculeInformation>* intermolecular_interactions,
                        gmx::ArrayRef<InteractionsOfType>     interactions,
-                       int*                                  comb,
+                       CombinationRule*                      comb,
                        double*                               reppow,
                        real*                                 fudgeQQ,
                        gmx_bool                              bMorse,
@@ -645,7 +647,7 @@ static void new_status(const char*                           topfile,
         convert_harmonics(*mi, atypes);
     }
 
-    if (ir->eDisre == edrNone)
+    if (ir->eDisre == DistanceRestraintRefinement::None)
     {
         i = rm_interactions(F_DISRES, *mi);
         if (i > 0)
@@ -678,7 +680,7 @@ static void new_status(const char*                           topfile,
      * constraints only. Do not print note with large timesteps or vsites.
      */
     if (opts->nshake == eshALLBONDS && ffParametrizedWithHBondConstraints && ir->delta_t < 0.0026
-        && gmx_mtop_ftype_count(sys, F_VSITE3FD) == 0)
+        && gmx_mtop_ftype_count(*sys, F_VSITE3FD) == 0)
     {
         set_warning_line(wi, "unknown", -1);
         warning_note(wi,
@@ -713,10 +715,10 @@ static void new_status(const char*                           topfile,
     /* It would be nice to get rid of the copies below, but we don't know
      * a priori if the number of atoms in confin matches what we expect.
      */
-    state->flags |= (1 << estX);
+    state->flags |= enumValueToBitMask(StateEntry::X);
     if (v != nullptr)
     {
-        state->flags |= (1 << estV);
+        state->flags |= enumValueToBitMask(StateEntry::V);
     }
     state_change_natoms(state, state->natoms);
     std::copy(x, x + state->natoms, state->x.data());
@@ -777,7 +779,7 @@ static void new_status(const char*                           topfile,
             opts->seed = static_cast<int>(gmx::makeRandomSeed());
             GMX_LOG(logger.info).asParagraph().appendTextFormatted("Setting gen_seed to %d", opts->seed);
         }
-        state->flags |= (1 << estV);
+        state->flags |= enumValueToBitMask(StateEntry::V);
         maxwell_speed(opts->tempi, opts->seed, sys, state->v.rvec_array(), logger);
 
         stop_cm(logger, state->natoms, mass, state->x.rvec_array(), state->v.rvec_array());
@@ -906,7 +908,7 @@ static void cont_status(const char*             slog,
     GMX_LOG(logger.info).asParagraph().appendTextFormatted("Using frame at t = %g ps", use_time);
     GMX_LOG(logger.info).asParagraph().appendTextFormatted("Starting time for run is %g ps", ir->init_t);
 
-    if ((ir->epc != epcNO || ir->etc == etcNOSEHOOVER) && ener)
+    if ((ir->epc != PressureCoupling::No || ir->etc == TemperatureCoupling::NoseHoover) && ener)
     {
         get_enx_state(ener, use_time, sys->groups, ir, state);
         preserve_box_shape(ir, state->box_rel, state->boxv);
@@ -917,7 +919,7 @@ static void read_posres(gmx_mtop_t*                              mtop,
                         gmx::ArrayRef<const MoleculeInformation> molinfo,
                         gmx_bool                                 bTopB,
                         const char*                              fn,
-                        int                                      rc_scaling,
+                        RefCoordScaling                          rc_scaling,
                         PbcType                                  pbcType,
                         rvec                                     com,
                         warninp*                                 wi,
@@ -954,7 +956,7 @@ static void read_posres(gmx_mtop_t*                              mtop,
     npbcdim = numPbcDimensions(pbcType);
     GMX_RELEASE_ASSERT(npbcdim <= DIM, "Invalid npbcdim");
     clear_rvec(com);
-    if (rc_scaling != erscNO)
+    if (rc_scaling != RefCoordScaling::No)
     {
         copy_mat(box, invbox);
         for (int j = npbcdim; j < DIM; j++)
@@ -992,7 +994,7 @@ static void read_posres(gmx_mtop_t*                              mtop,
                               natoms);
                 }
                 hadAtom[ai] = TRUE;
-                if (rc_scaling == erscCOM)
+                if (rc_scaling == RefCoordScaling::Com)
                 {
                     /* Determine the center of mass of the posres reference coordinates */
                     for (int j = 0; j < npbcdim; j++)
@@ -1016,7 +1018,7 @@ static void read_posres(gmx_mtop_t*                              mtop,
                               fn,
                               natoms);
                 }
-                if (rc_scaling == erscCOM && !hadAtom[ai])
+                if (rc_scaling == RefCoordScaling::Com && !hadAtom[ai])
                 {
                     /* Determine the center of mass of the posres reference coordinates */
                     for (int j = 0; j < npbcdim; j++)
@@ -1045,7 +1047,7 @@ static void read_posres(gmx_mtop_t*                              mtop,
         }
         a += nat_molb;
     }
-    if (rc_scaling == erscCOM)
+    if (rc_scaling == RefCoordScaling::Com)
     {
         if (totmass == 0)
         {
@@ -1064,7 +1066,7 @@ static void read_posres(gmx_mtop_t*                              mtop,
                         com[ZZ]);
     }
 
-    if (rc_scaling != erscNO)
+    if (rc_scaling != RefCoordScaling::No)
     {
         GMX_ASSERT(npbcdim <= DIM, "Only DIM dimensions can have PBC");
 
@@ -1078,7 +1080,7 @@ static void read_posres(gmx_mtop_t*                              mtop,
                 {
                     for (int j = 0; j < npbcdim; j++)
                     {
-                        if (rc_scaling == erscALL)
+                        if (rc_scaling == RefCoordScaling::All)
                         {
                             /* Convert from Cartesian to crystal coordinates */
                             xp[i][j] *= invbox[j][j];
@@ -1087,7 +1089,7 @@ static void read_posres(gmx_mtop_t*                              mtop,
                                 xp[i][j] += invbox[k][j] * xp[i][k];
                             }
                         }
-                        else if (rc_scaling == erscCOM)
+                        else if (rc_scaling == RefCoordScaling::Com)
                         {
                             /* Subtract the center of mass */
                             xp[i][j] -= com[j];
@@ -1097,7 +1099,7 @@ static void read_posres(gmx_mtop_t*                              mtop,
             }
         }
 
-        if (rc_scaling == erscCOM)
+        if (rc_scaling == RefCoordScaling::Com)
         {
             /* Convert the COM from Cartesian to crystal coordinates */
             for (int j = 0; j < npbcdim; j++)
@@ -1120,7 +1122,7 @@ static void gen_posres(gmx_mtop_t*                              mtop,
                        gmx::ArrayRef<const MoleculeInformation> mi,
                        const char*                              fnA,
                        const char*                              fnB,
-                       int                                      rc_scaling,
+                       RefCoordScaling                          rc_scaling,
                        PbcType                                  pbcType,
                        rvec                                     com,
                        rvec                                     comB,
@@ -1148,13 +1150,17 @@ static void set_wall_atomtype(PreprocessingAtomTypes* at,
     }
     for (i = 0; i < ir->nwall; i++)
     {
-        ir->wall_atomtype[i] = at->atomTypeFromName(opts->wall_atomtype[i]);
-        if (ir->wall_atomtype[i] == NOTSET)
+        auto atomType = at->atomTypeFromName(opts->wall_atomtype[i]);
+        if (!atomType.has_value())
         {
             std::string warningMessage = gmx::formatString(
                     "Specified wall atom type %s is not defined", opts->wall_atomtype[i]);
             warning_error(wi, warningMessage.c_str());
         }
+        else
+        {
+            ir->wall_atomtype[i] = *atomType;
+        }
     }
 }
 
@@ -1166,7 +1172,8 @@ static int nrdf_internal(const t_atoms* atoms)
     for (i = 0; i < atoms->nr; i++)
     {
         /* Vsite ptype might not be set here yet, so also check the mass */
-        if ((atoms->atom[i].ptype == eptAtom || atoms->atom[i].ptype == eptNucleus) && atoms->atom[i].m > 0)
+        if ((atoms->atom[i].ptype == ParticleType::Atom || atoms->atom[i].ptype == ParticleType::Nucleus)
+            && atoms->atom[i].m > 0)
         {
             nmass++;
         }
@@ -1377,7 +1384,7 @@ static real calc_temp(const gmx_mtop_t* mtop, const t_inputrec* ir, rvec* v)
         nrdf += ir->opts.nrdf[g];
     }
 
-    return sum_mv2 / (nrdf * BOLTZ);
+    return sum_mv2 / (nrdf * gmx::c_boltz);
 }
 
 static real get_max_reference_temp(const t_inputrec* ir, warninp* wi)
@@ -1449,7 +1456,7 @@ static void checkForUnboundAtoms(const gmx_moltype_t* molt, gmx_bool bVerbose, w
     int numDanglingAtoms = 0;
     for (int a = 0; a < atoms->nr; a++)
     {
-        if (atoms->atom[a].ptype != eptVSite && count[a] == 0)
+        if (atoms->atom[a].ptype != ParticleType::VSite && count[a] == 0)
         {
             if (bVerbose)
             {
@@ -1570,7 +1577,7 @@ static bool haveDecoupledModeInMol(const gmx_moltype_t&           molt,
 
 /*! \brief Checks if the Verlet buffer and constraint accuracy is sufficient for decoupled dynamic modes.
  *
- * When decoupled modes are present and the accuray in insufficient,
+ * When decoupled modes are present and the accuracy in insufficient,
  * this routine issues a warning if the accuracy is insufficient.
  */
 static void checkDecoupledModeAccuracy(const gmx_mtop_t* mtop, const t_inputrec* ir, warninp* wi)
@@ -1604,11 +1611,12 @@ static void checkDecoupledModeAccuracy(const gmx_mtop_t* mtop, const t_inputrec*
     const int  lincsOrderThreshold      = 4;
     const real shakeToleranceThreshold  = 0.005 * ir->delta_t;
 
-    bool lincsWithSufficientTolerance = (ir->eConstrAlg == econtLINCS && ir->nLincsIter >= lincsIterationThreshold
-                                         && ir->nProjOrder >= lincsOrderThreshold);
-    bool shakeWithSufficientTolerance =
-            (ir->eConstrAlg == econtSHAKE && ir->shake_tol <= 1.1 * shakeToleranceThreshold);
-    if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol <= 1.1 * bufferToleranceThreshold
+    bool lincsWithSufficientTolerance =
+            (ir->eConstrAlg == ConstraintAlgorithm::Lincs
+             && ir->nLincsIter >= lincsIterationThreshold && ir->nProjOrder >= lincsOrderThreshold);
+    bool shakeWithSufficientTolerance = (ir->eConstrAlg == ConstraintAlgorithm::Shake
+                                         && ir->shake_tol <= 1.1 * shakeToleranceThreshold);
+    if (ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol <= 1.1 * bufferToleranceThreshold
         && (lincsWithSufficientTolerance || shakeWithSufficientTolerance))
     {
         return;
@@ -1630,7 +1638,7 @@ static void checkDecoupledModeAccuracy(const gmx_mtop_t* mtop, const t_inputrec*
                 "and with masses that differ by more than a factor of %g. This means "
                 "that there are likely dynamic modes that are only very weakly coupled.",
                 massFactorThreshold);
-        if (ir->cutoff_scheme == ecutsVERLET)
+        if (ir->cutoff_scheme == CutoffScheme::Verlet)
         {
             message += gmx::formatString(
                     " To ensure good equipartitioning, you need to either not use "
@@ -1638,7 +1646,7 @@ static void checkDecoupledModeAccuracy(const gmx_mtop_t* mtop, const t_inputrec*
                     "hydrogens) or use integrator = %s or decrease one or more tolerances: "
                     "verlet-buffer-tolerance <= %g, LINCS iterations >= %d, LINCS order "
                     ">= %d or SHAKE tolerance <= %g",
-                    ei_names[eiSD1],
+                    enumValueToString(IntegrationAlgorithm::SD1),
                     bufferToleranceThreshold,
                     lincsIterationThreshold,
                     lincsOrderThreshold,
@@ -1650,7 +1658,7 @@ static void checkDecoupledModeAccuracy(const gmx_mtop_t* mtop, const t_inputrec*
                     " To ensure good equipartitioning, we suggest to switch to the %s "
                     "cutoff-scheme, since that allows for better control over the Verlet "
                     "buffer size and thus over the energy drift.",
-                    ecutscheme_names[ecutsVERLET]);
+                    enumValueToString(CutoffScheme::Verlet));
         }
         warning(wi, message);
     }
@@ -1742,9 +1750,7 @@ int gmx_grompp(int argc, char* argv[])
         "Then a coordinate file is read and velocities can be generated",
         "from a Maxwellian distribution if requested.",
         "[THISMODULE] also reads parameters for [gmx-mdrun] ",
-        "(eg. number of MD steps, time step, cut-off), and others such as",
-        "NEMD parameters, which are corrected so that the net acceleration",
-        "is zero.",
+        "(eg. number of MD steps, time step, cut-off).",
         "Eventually a binary file is produced that can serve as the sole input",
         "file for the MD program.[PAR]",
 
@@ -1829,7 +1835,8 @@ int gmx_grompp(int argc, char* argv[])
     };
     std::vector<MoleculeInformation>     mi;
     std::unique_ptr<MoleculeInformation> intermolecular_interactions;
-    int                                  nvsite, comb;
+    int                                  nvsite;
+    CombinationRule                      comb;
     real                                 fudgeQQ;
     double                               reppow;
     const char*                          mdparin;
@@ -1851,7 +1858,9 @@ int gmx_grompp(int argc, char* argv[])
                        { efEDR, "-e", nullptr, ffOPTRD },
                        /* This group is needed by the VMD viewer as the start configuration for IMD sessions: */
                        { efGRO, "-imd", "imdgroup", ffOPTWR },
-                       { efTRN, "-ref", "rotref", ffOPTRW | ffALLOW_MISSING } };
+                       { efTRN, "-ref", "rotref", ffOPTRW | ffALLOW_MISSING },
+                       /* This group is needed by the QMMM MDModule: */
+                       { efQMI, "-qmi", nullptr, ffOPTRD } };
 #define NFILE asize(fnm)
 
     /* Command line options */
@@ -1919,10 +1928,21 @@ int gmx_grompp(int argc, char* argv[])
     }
     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 
-    // Now that the MdModules have their options assigned from get_ir, subscribe
+    // Now that the MDModules have their options assigned from get_ir, subscribe
     // to eventual notifications during pre-processing their data
     mdModules.subscribeToPreProcessingNotifications();
 
+    // Notify MDModules of existing logger
+    mdModules.notifiers().preProcessingNotifier_.notify(logger);
+
+    // Notify MDModules of existing warninp
+    mdModules.notifiers().preProcessingNotifier_.notify(wi);
+
+    // Notify QMMM MDModule of external QM input file command-line option
+    {
+        gmx::QMInputFileName qmInputFileName = { ftp2bSet(efQMI, NFILE, fnm), ftp2fn(efQMI, NFILE, fnm) };
+        mdModules.notifiers().preProcessingNotifier_.notify(qmInputFileName);
+    }
 
     if (bVerbose)
     {
@@ -1930,7 +1950,7 @@ int gmx_grompp(int argc, char* argv[])
                 .asParagraph()
                 .appendTextFormatted("checking input for internal consistency...");
     }
-    check_ir(mdparin, mdModules.notifier(), ir, opts, wi);
+    check_ir(mdparin, mdModules.notifiers(), ir, opts, wi);
 
     if (ir->ld_seed == -1)
     {
@@ -2016,27 +2036,28 @@ int gmx_grompp(int argc, char* argv[])
         }
     }
 
-    if ((count_constraints(&sys, mi, wi) != 0) && (ir->eConstrAlg == econtSHAKE))
+    if ((count_constraints(&sys, mi, wi) != 0) && (ir->eConstrAlg == ConstraintAlgorithm::Shake))
     {
-        if (ir->eI == eiCG || ir->eI == eiLBFGS)
+        if (ir->eI == IntegrationAlgorithm::CG || ir->eI == IntegrationAlgorithm::LBFGS)
         {
-            std::string warningMessage = gmx::formatString("Can not do %s with %s, use %s",
-                                                           EI(ir->eI),
-                                                           econstr_names[econtSHAKE],
-                                                           econstr_names[econtLINCS]);
+            std::string warningMessage =
+                    gmx::formatString("Can not do %s with %s, use %s",
+                                      enumValueToString(ir->eI),
+                                      enumValueToString(ConstraintAlgorithm::Shake),
+                                      enumValueToString(ConstraintAlgorithm::Lincs));
             warning_error(wi, warningMessage);
         }
         if (ir->bPeriodicMols)
         {
             std::string warningMessage =
                     gmx::formatString("Can not do periodic molecules with %s, use %s",
-                                      econstr_names[econtSHAKE],
-                                      econstr_names[econtLINCS]);
+                                      enumValueToString(ConstraintAlgorithm::Shake),
+                                      enumValueToString(ConstraintAlgorithm::Lincs));
             warning_error(wi, warningMessage);
         }
     }
 
-    if (EI_SD(ir->eI) && ir->etc != etcNO)
+    if (EI_SD(ir->eI) && ir->etc != TemperatureCoupling::No)
     {
         warning_note(wi, "Temperature coupling is ignored with SD integrators.");
     }
@@ -2048,14 +2069,14 @@ int gmx_grompp(int argc, char* argv[])
 
     if (nint_ftype(&sys, mi, F_POSRES) > 0 || nint_ftype(&sys, mi, F_FBPOSRES) > 0)
     {
-        if (ir->epc == epcPARRINELLORAHMAN || ir->epc == epcMTTK)
+        if (ir->epc == PressureCoupling::ParrinelloRahman || ir->epc == PressureCoupling::Mttk)
         {
             std::string warningMessage = gmx::formatString(
                     "You are combining position restraints with %s pressure coupling, which can "
                     "lead to instabilities. If you really want to combine position restraints with "
                     "pressure coupling, we suggest to use %s pressure coupling instead.",
-                    EPCOUPLTYPE(ir->epc),
-                    EPCOUPLTYPE(epcBERENDSEN));
+                    enumValueToString(ir->epc),
+                    enumValueToString(PressureCoupling::Berendsen));
             warning_note(wi, warningMessage);
         }
 
@@ -2168,9 +2189,16 @@ int gmx_grompp(int argc, char* argv[])
     /* check masses */
     check_mol(&sys, wi);
 
+    if (haveFepPerturbedMassesInSettles(sys))
+    {
+        warning_error(wi,
+                      "SETTLE is not implemented for atoms whose mass is perturbed. "
+                      "You might instead use normal constraints.");
+    }
+
     checkForUnboundAtoms(&sys, bVerbose, wi, logger);
 
-    if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD)
+    if (EI_DYNAMICS(ir->eI) && ir->eI != IntegrationAlgorithm::BD)
     {
         check_bonds_timestep(&sys, ir->delta_t, wi);
     }
@@ -2193,15 +2221,18 @@ int gmx_grompp(int argc, char* argv[])
     {
         GMX_LOG(logger.info).asParagraph().appendTextFormatted("initialising group options...");
     }
-    do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm), &sys, bVerbose, mdModules.notifier(), ir, wi);
+    do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm), &sys, bVerbose, mdModules.notifiers(), ir, wi);
+
+    // Notify topology to MdModules for pre-processing after all indexes were built
+    mdModules.notifiers().preProcessingNotifier_.notify(&sys);
 
-    if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0)
+    if (ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol > 0)
     {
         if (EI_DYNAMICS(ir->eI) && inputrec2nboundeddim(ir) == 3)
         {
             real buffer_temp;
 
-            if (EI_MD(ir->eI) && ir->etc == etcNO)
+            if (EI_MD(ir->eI) && ir->etc == TemperatureCoupling::No)
             {
                 if (bGenVel)
                 {
@@ -2233,7 +2264,7 @@ int gmx_grompp(int argc, char* argv[])
                 buffer_temp = get_max_reference_temp(ir, wi);
             }
 
-            if (EI_MD(ir->eI) && ir->etc == etcNO && buffer_temp == 0)
+            if (EI_MD(ir->eI) && ir->etc == TemperatureCoupling::No && buffer_temp == 0)
             {
                 /* NVE with initial T=0: we add a fixed ratio to rlist.
                  * Since we don't actually use verletbuf_tol, we set it to -1
@@ -2249,14 +2280,14 @@ int gmx_grompp(int argc, char* argv[])
                  * Note that we can't warn when nsteps=0, since we don't
                  * know how many steps the user intends to run.
                  */
-                if (EI_MD(ir->eI) && ir->etc == etcNO && ir->nstlist > 1 && ir->nsteps > 0)
+                if (EI_MD(ir->eI) && ir->etc == TemperatureCoupling::No && ir->nstlist > 1 && ir->nsteps > 0)
                 {
                     const real driftTolerance = 0.01;
                     /* We use 2 DOF per atom = 2kT pot+kin energy,
                      * to be on the safe side with constraints.
                      */
                     const real totalEnergyDriftPerAtomPerPicosecond =
-                            2 * BOLTZ * buffer_temp / (ir->nsteps * ir->delta_t);
+                            2 * gmx::c_boltz * buffer_temp / (ir->nsteps * ir->delta_t);
 
                     if (ir->verletbuf_tol > 1.1 * driftTolerance * totalEnergyDriftPerAtomPerPicosecond)
                     {
@@ -2294,7 +2325,7 @@ int gmx_grompp(int argc, char* argv[])
         pr_symtab(debug, 0, "After close", &sys.symtab);
     }
 
-    if (ir->eI == eiMimic)
+    if (ir->eI == IntegrationAlgorithm::Mimic)
     {
         generate_qmexcl(&sys, ir, logger);
     }
@@ -2328,7 +2359,7 @@ int gmx_grompp(int argc, char* argv[])
     {
         /* Calculate the optimal grid dimensions */
         matrix          scaledBox;
-        EwaldBoxZScaler boxScaler(*ir);
+        EwaldBoxZScaler boxScaler(inputrecPbcXY2Walls(ir), ir->wall_ewald_zfac);
         boxScaler.scaleBox(state.box, scaledBox);
 
         if (ir->nkx > 0 && ir->nky > 0 && ir->nkz > 0)
@@ -2355,10 +2386,10 @@ int gmx_grompp(int argc, char* argv[])
     /* MRS: eventually figure out better logic for initializing the fep
        values that makes declaring the lambda and declaring the state not
        potentially conflict if not handled correctly. */
-    if (ir->efep != efepNO)
+    if (ir->efep != FreeEnergyPerturbationType::No)
     {
         state.fep_state = ir->fepvals->init_fep_state;
-        for (i = 0; i < efptNR; i++)
+        for (i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
         {
             /* init_lambda trumps state definitions*/
             if (ir->fepvals->init_lambda >= 0)
@@ -2367,7 +2398,7 @@ int gmx_grompp(int argc, char* argv[])
             }
             else
             {
-                if (ir->fepvals->all_lambda[i] == nullptr)
+                if (ir->fepvals->all_lambda[i].empty())
                 {
                     gmx_fatal(FARGS, "Values of lambda not set for a free energy calculation!");
                 }
@@ -2383,7 +2414,8 @@ int gmx_grompp(int argc, char* argv[])
 
     if (ir->bPull)
     {
-        pull = set_pull_init(ir, &sys, state.x.rvec_array(), state.box, state.lambda[efptMASS], wi);
+        pull = set_pull_init(
+                ir, sys, state.x, state.box, state.lambda[FreeEnergyPerturbationCouplingType::Mass], wi);
     }
 
     /* Modules that supply external potential for pull coordinates
@@ -2394,19 +2426,21 @@ int gmx_grompp(int argc, char* argv[])
     if (ir->bDoAwh)
     {
         tensor compressibility = { { 0 } };
-        if (ir->epc != epcNO)
+        if (ir->epc != PressureCoupling::No)
         {
             copy_mat(ir->compress, compressibility);
         }
         setStateDependentAwhParams(
-                ir->awhParams,
+                ir->awhParams.get(),
                 *ir->pull,
                 pull,
                 state.box,
                 ir->pbcType,
                 compressibility,
                 &ir->opts,
-                ir->efep != efepNO ? ir->fepvals->all_lambda[efptFEP][ir->fepvals->init_fep_state] : 0,
+                ir->efep != FreeEnergyPerturbationType::No ? ir->fepvals->all_lambda[static_cast<int>(
+                        FreeEnergyPerturbationCouplingType::Fep)][ir->fepvals->init_fep_state]
+                                                           : 0,
                 sys,
                 wi);
     }
@@ -2439,14 +2473,15 @@ int gmx_grompp(int argc, char* argv[])
          * charges. This will double the cost, but the optimal performance will
          * then probably be at a slightly larger cut-off and grid spacing.
          */
-        if ((ir->efep == efepNO && ratio > 1.0 / 2.0) || (ir->efep != efepNO && ratio > 2.0 / 3.0))
+        if ((ir->efep == FreeEnergyPerturbationType::No && ratio > 1.0 / 2.0)
+            || (ir->efep != FreeEnergyPerturbationType::No && ratio > 2.0 / 3.0))
         {
             warning_note(
                     wi,
                     "The optimal PME mesh load for parallel simulations is below 0.5\n"
                     "and for highly parallel simulations between 0.25 and 0.33,\n"
                     "for higher performance, increase the cut-off and the PME grid spacing.\n");
-            if (ir->efep != efepNO)
+            if (ir->efep != FreeEnergyPerturbationType::No)
             {
                 warning_note(wi,
                              "For free energy simulations, the optimal load limit increases from "
@@ -2470,31 +2505,61 @@ int gmx_grompp(int argc, char* argv[])
         }
     }
 
+    // Hand over box and coordiantes to MdModules before they evaluate their final parameters
+    {
+        gmx::CoordinatesAndBoxPreprocessed coordinatesAndBoxPreprocessed;
+        coordinatesAndBoxPreprocessed.coordinates_ = state.x.arrayRefWithPadding();
+        copy_mat(state.box, coordinatesAndBoxPreprocessed.box_);
+        coordinatesAndBoxPreprocessed.pbc_ = ir->pbcType;
+        mdModules.notifiers().preProcessingNotifier_.notify(coordinatesAndBoxPreprocessed);
+    }
+
     // Add the md modules internal parameters that are not mdp options
     // e.g., atom indices
 
     {
         gmx::KeyValueTreeBuilder internalParameterBuilder;
-        mdModules.notifier().preProcessingNotifications_.notify(internalParameterBuilder.rootObject());
+        mdModules.notifiers().preProcessingNotifier_.notify(internalParameterBuilder.rootObject());
         ir->internalParameters =
                 std::make_unique<gmx::KeyValueTreeObject>(internalParameterBuilder.build());
     }
 
+    if (ir->comm_mode != ComRemovalAlgorithm::No)
+    {
+        const int nstglobalcomm = computeGlobalCommunicationPeriod(ir);
+        if (ir->nstcomm % nstglobalcomm != 0)
+        {
+            warning_note(
+                    wi,
+                    gmx::formatString(
+                            "COM removal frequency is set to (%d).\n"
+                            "Other settings require a global communication frequency of %d.\n"
+                            "Note that this will require additional global communication steps,\n"
+                            "which will reduce performance when using multiple ranks.\n"
+                            "Consider setting nstcomm to a multiple of %d.",
+                            ir->nstcomm,
+                            nstglobalcomm,
+                            nstglobalcomm));
+        }
+    }
+
     if (bVerbose)
     {
         GMX_LOG(logger.info).asParagraph().appendTextFormatted("writing run input file...");
     }
 
     done_warning(wi, FARGS);
-    write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, &sys);
+    write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, sys);
 
     /* Output IMD group, if bIMD is TRUE */
-    gmx::write_IMDgroup_to_file(ir->bIMD, ir, &state, &sys, NFILE, fnm);
+    gmx::write_IMDgroup_to_file(ir->bIMD, ir, &state, sys, NFILE, fnm);
 
     sfree(opts->define);
     sfree(opts->wall_atomtype[0]);
     sfree(opts->wall_atomtype[1]);
     sfree(opts->include);
+    sfree(opts->couple_moltype);
+
     for (auto& mol : mi)
     {
         // Some of the contents of molinfo have been stolen, so