Reimplement constant acceleration groups

[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / readir.cpp

diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp
index ff47ab5a35ea0e2ee71578db6db51c298f087985..557a361c8f2f472ffc94c0b719a7a6633d9b2cde 100644 (file)
--- a/src/gromacs/gmxpreprocess/readir.cpp
+++ b/src/gromacs/gmxpreprocess/readir.cpp
@@ -107,10 +107,11 @@ using gmx::BasicVector;
 
 struct gmx_inputrec_strings
 {
-    char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], freeze[STRLEN], frdim[STRLEN],
-            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 tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], accelerationGroups[STRLEN],
+            acceleration[STRLEN], freeze[STRLEN], frdim[STRLEN], 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];
     gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, std::string> fep_lambda;
     char                                                                   lambda_weights[STRLEN];
     std::vector<std::string>                                               pullGroupNames;
@@ -1771,6 +1772,33 @@ static void convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs,
     }
 }
 
+static void convertRvecs(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, rvec* outputs)
+{
+    int i = 0, d = 0;
+    for (const auto& input : inputs)
+    {
+        try
+        {
+            outputs[i][d] = gmx::fromString<real>(input);
+        }
+        catch (gmx::GromacsException&)
+        {
+            auto message = gmx::formatString(
+                    "Invalid value for mdp option %s. %s should only consist of real numbers "
+                    "separated by spaces.",
+                    name,
+                    name);
+            warning_error(wi, message);
+        }
+        ++d;
+        if (d == DIM)
+        {
+            d = 0;
+            ++i;
+        }
+    }
+}
+
 static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_density, t_gromppopts* opts, warninp_t wi)
 {
     opts->wall_atomtype[0] = nullptr;
@@ -2340,6 +2368,8 @@ void get_ir(const char*     mdparin,
 
     /* Non-equilibrium MD stuff */
     printStringNewline(&inp, "Non-equilibrium MD stuff");
+    setStringEntry(&inp, "acc-grps", inputrecStrings->accelerationGroups, nullptr);
+    setStringEntry(&inp, "accelerate", inputrecStrings->acceleration, nullptr);
     setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
     setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
     ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
@@ -3907,6 +3937,31 @@ void do_index(const char*                    mdparin,
             *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
     mdModulesNotifiers.preProcessingNotifier_.notify(defaultIndexGroupsAndNames);
 
+    auto accelerations          = gmx::splitString(inputrecStrings->acceleration);
+    auto accelerationGroupNames = gmx::splitString(inputrecStrings->accelerationGroups);
+    if (accelerationGroupNames.size() * DIM != accelerations.size())
+    {
+        gmx_fatal(FARGS,
+                  "Invalid Acceleration input: %zu groups and %zu acc. values",
+                  accelerationGroupNames.size(),
+                  accelerations.size());
+    }
+    do_numbering(natoms,
+                 groups,
+                 accelerationGroupNames,
+                 defaultIndexGroups,
+                 gnames,
+                 SimulationAtomGroupType::Acceleration,
+                 restnm,
+                 GroupCoverage::AllGenerateRest,
+                 bVerbose,
+                 wi);
+    nr = groups->groups[SimulationAtomGroupType::Acceleration].size();
+    snew(ir->opts.acceleration, nr);
+    ir->opts.ngacc = nr;
+
+    convertRvecs(wi, accelerations, "accelerations", ir->opts.acceleration);
+
     auto freezeDims       = gmx::splitString(inputrecStrings->frdim);
     auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
     if (freezeDims.size() != DIM * freezeGroupNames.size())
@@ -4370,6 +4425,7 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_
     char                      err_buf[STRLEN];
     int                       i, m, c, nmol;
     bool                      bCharge;
+    real *                    mgrp, mt;
     gmx_mtop_atomloop_block_t aloopb;
     char                      warn_buf[STRLEN];
 
@@ -4571,6 +4627,57 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_
                       "constant by hand.");
     }
 
+    ir->useConstantAcceleration = false;
+    for (int i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
+    {
+        if (norm2(ir->opts.acceleration[i]) != 0)
+        {
+            ir->useConstantAcceleration = true;
+        }
+    }
+    if (ir->useConstantAcceleration)
+    {
+        gmx::RVec acceleration = { 0.0_real, 0.0_real, 0.0_real };
+        snew(mgrp, sys->groups.groups[SimulationAtomGroupType::Acceleration].size());
+        for (const AtomProxy atomP : AtomRange(*sys))
+        {
+            const t_atom& local = atomP.atom();
+            int           i     = atomP.globalAtomNumber();
+            mgrp[getGroupType(sys->groups, SimulationAtomGroupType::Acceleration, i)] += local.m;
+        }
+        mt = 0.0;
+        for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
+        {
+            for (m = 0; (m < DIM); m++)
+            {
+                acceleration[m] += ir->opts.acceleration[i][m] * mgrp[i];
+            }
+            mt += mgrp[i];
+        }
+        for (m = 0; (m < DIM); m++)
+        {
+            if (fabs(acceleration[m]) > 1e-6)
+            {
+                const char* dim[DIM] = { "X", "Y", "Z" };
+                fprintf(stderr,
+                        "Net Acceleration in %s direction, will %s be corrected\n",
+                        dim[m],
+                        ir->nstcomm != 0 ? "" : "not");
+                if (ir->nstcomm != 0 && m < ndof_com(ir))
+                {
+                    acceleration[m] /= mt;
+                    for (i = 0;
+                         (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration]));
+                         i++)
+                    {
+                        ir->opts.acceleration[i][m] -= acceleration[m];
+                    }
+                }
+            }
+        }
+        sfree(mgrp);
+    }
+
     if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->sc_alpha != 0
         && !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
     {