Fix random typos
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / grompp.cpp
index 97d6d6ffcfe31e598e52d0a6f7d9bf4ad068b8d0..899b5008c2becadbe5008ba31f7652865f5c7fa8 100644 (file)
@@ -39,6 +39,7 @@
 
 #include "grompp.h"
 
+#include <array>
 #include <cerrno>
 #include <climits>
 #include <cmath>
 #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;
@@ -468,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]);
         }
@@ -482,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++;
         }
@@ -679,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,
@@ -1149,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;
+        }
     }
 }
 
@@ -1167,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++;
         }
@@ -1378,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)
@@ -1450,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)
             {
@@ -1571,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)
@@ -1852,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 */
@@ -1920,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)
     {
@@ -1931,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)
     {
@@ -2170,6 +2189,13 @@ 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 != IntegrationAlgorithm::BD)
@@ -2195,7 +2221,10 @@ 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 == CutoffScheme::Verlet && ir->verletbuf_tol > 0)
     {
@@ -2258,7 +2287,7 @@ int gmx_grompp(int argc, char* argv[])
                      * 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)
                     {
@@ -2330,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)
@@ -2385,12 +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[FreeEnergyPerturbationCouplingType::Mass],
-                             wi);
+        pull = set_pull_init(
+                ir, sys, state.x, state.box, state.lambda[FreeEnergyPerturbationCouplingType::Mass], wi);
     }
 
     /* Modules that supply external potential for pull coordinates
@@ -2406,17 +2431,16 @@ int gmx_grompp(int argc, char* argv[])
             copy_mat(ir->compress, compressibility);
         }
         setStateDependentAwhParams(
-                ir->awhParams,
+                ir->awhParams.get(),
                 *ir->pull,
                 pull,
                 state.box,
                 ir->pbcType,
                 compressibility,
                 &ir->opts,
-                ir->efep != FreeEnergyPerturbationType::No
-                        ? ir->fepvals->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Fep)]
-                                                 [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);
     }
@@ -2481,12 +2505,21 @@ 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());
     }
@@ -2516,15 +2549,17 @@ int gmx_grompp(int argc, char* argv[])
     }
 
     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