Fix center of mass motion removal with frozen atoms

When frozen atoms were part of center of mass motion removal groups,
they would still contribute to the mass of those groups. This meant
that the COM velocity correction was (slightly) too small. Now
completely frozen atoms are removed from COM removal groups by grompp.
When atoms are only frozen along one or two dimensions and part of
a COM removal group, grompp now issues a warning.

Also fixed an nullptr or incorrect string buffer passed to warning()
with invalid freeze group dimension user input.

Fixes #2553

Change-Id: I20a03fea511e75a131cb27880acc1f4ee4a2bfb8

diff --git a/docs/release-notes/2020/2020.1.rst b/docs/release-notes/2020/2020.1.rst
index 6467742572fa17d6f79bd61be9028654ec13f613..b75be1badf49ef85594bddd6d26b5ca5ab871956 100644 (file)
--- a/docs/release-notes/2020/2020.1.rst
+++ b/docs/release-notes/2020/2020.1.rst
@@ -147,6 +147,18 @@ segmentation fault.
 
 :issue:`3389`
 
+Fix center of mass motion removal with frozen atoms
+"""""""""""""""""""""""""""""""""""""""""""""""""""
+
+When frozen atoms were part of center of mass motion removal groups,
+they would still contribute to the mass of those groups. This meant that
+the COM velocity correction was (slightly) too small. Now completely
+frozen atoms are removed from COM removal groups by grompp.
+When atoms are only frozen along one or two dimensions and part of
+a COM removal group, grompp now issues a warning.
+
+:issue:`2553`
+
 Fix temperature calculation when center of mass motion is removed for part of the system
 """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 

diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp
index 46d8fecb80baa82b6861c8af328be4ebb87d2438..5b662c6e31e8689c1e8d9f2a005283fa5961ce99 100644 (file)
--- a/src/gromacs/gmxpreprocess/readir.cpp
+++ b/src/gromacs/gmxpreprocess/readir.cpp
@@ -2706,7 +2706,7 @@ int search_string(const char* s, int ng, char* gn[])
               s);
 }
 
-static bool do_numbering(int                        natoms,
+static void do_numbering(int                        natoms,
                          SimulationGroups*          groups,
                          gmx::ArrayRef<std::string> groupsFromMdpFile,
                          t_blocka*                  block,
@@ -2721,7 +2721,6 @@ static bool do_numbering(int                        natoms,
     AtomGroupIndices* grps = &(groups->groups[gtype]);
     int               j, gid, aj, ognr, ntot = 0;
     const char*       title;
-    bool              bRest;
     char              warn_buf[STRLEN];
 
     title = shortName(gtype);
@@ -2777,7 +2776,6 @@ static bool do_numbering(int                        natoms,
     }
 
     /* Now check whether we have done all atoms */
-    bRest = FALSE;
     if (ntot != natoms)
     {
         if (grptp == egrptpALL)
@@ -2795,7 +2793,6 @@ static bool do_numbering(int                        natoms,
             if (cbuf[j] == NOGID)
             {
                 cbuf[j] = grps->size();
-                bRest   = TRUE;
             }
         }
         if (grptp != egrptpPART)
@@ -2834,8 +2831,6 @@ static bool do_numbering(int                        natoms,
     }
 
     sfree(cbuf);
-
-    return (bRest && grptp == egrptpPART);
 }
 
 static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
@@ -3238,6 +3233,88 @@ static void make_IMD_group(t_IMD* IMDgroup, char* IMDgname, t_blocka* grps, char
     }
 }
 
+/* Checks whether atoms are both part of a COM removal group and frozen.
+ * If a fully frozen atom is part of a COM removal group, it is removed
+ * from the COM removal group. A note is issued if such atoms are present.
+ * A warning is issued for atom with one or two dimensions frozen that
+ * are part of a COM removal group (mdrun would need to compute COM mass
+ * per dimension to handle this correctly).
+ * Also issues a warning when non-frozen atoms are not part of a COM
+ * removal group while COM removal is active.
+ */
+static void checkAndUpdateVcmFreezeGroupConsistency(SimulationGroups* groups,
+                                                    const int         numAtoms,
+                                                    const t_grpopts&  opts,
+                                                    warninp_t         wi)
+{
+    const int vcmRestGroup =
+            std::max(int(groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()), 1);
+
+    int numFullyFrozenVcmAtoms     = 0;
+    int numPartiallyFrozenVcmAtoms = 0;
+    int numNonVcmAtoms             = 0;
+    for (int a = 0; a < numAtoms; a++)
+    {
+        const int freezeGroup   = getGroupType(*groups, SimulationAtomGroupType::Freeze, a);
+        int       numFrozenDims = 0;
+        for (int d = 0; d < DIM; d++)
+        {
+            numFrozenDims += opts.nFreeze[freezeGroup][d];
+        }
+
+        const int vcmGroup = getGroupType(*groups, SimulationAtomGroupType::MassCenterVelocityRemoval, a);
+        if (vcmGroup < vcmRestGroup)
+        {
+            if (numFrozenDims == DIM)
+            {
+                /* Do not remove COM motion for this fully frozen atom */
+                if (groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].empty())
+                {
+                    groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].resize(numAtoms, 0);
+                }
+                groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval][a] = vcmRestGroup;
+                numFullyFrozenVcmAtoms++;
+            }
+            else if (numFrozenDims > 0)
+            {
+                numPartiallyFrozenVcmAtoms++;
+            }
+        }
+        else if (numFrozenDims < DIM)
+        {
+            numNonVcmAtoms++;
+        }
+    }
+
+    if (numFullyFrozenVcmAtoms > 0)
+    {
+        std::string warningText = gmx::formatString(
+                "There are %d atoms that are fully frozen and part of COMM removal group(s), "
+                "removing these atoms from the COMM removal group(s)",
+                numFullyFrozenVcmAtoms);
+        warning_note(wi, warningText.c_str());
+    }
+    if (numPartiallyFrozenVcmAtoms > 0 && numPartiallyFrozenVcmAtoms < numAtoms)
+    {
+        std::string warningText = gmx::formatString(
+                "There are %d atoms that are frozen along less then %d dimensions and part of COMM "
+                "removal group(s), due to limitations in the code these still contribute to the "
+                "mass of the COM along frozen dimensions and therefore the COMM correction will be "
+                "too small.",
+                numPartiallyFrozenVcmAtoms, DIM);
+        warning(wi, warningText.c_str());
+    }
+    if (numNonVcmAtoms > 0)
+    {
+        std::string warningText = gmx::formatString(
+                "%d atoms are not part of any center of mass motion removal group.\n"
+                "This may lead to artifacts.\n"
+                "In most cases one should use one group for the whole system.",
+                numNonVcmAtoms);
+        warning(wi, warningText.c_str());
+    }
+}
+
 void do_index(const char*                   mdparin,
               const char*                   ndx,
               gmx_mtop_t*                   mtop,
@@ -3250,12 +3327,12 @@ void do_index(const char*                   mdparin,
     int       natoms;
     t_symtab* symtab;
     t_atoms   atoms_all;
-    char      warnbuf[STRLEN], **gnames;
+    char**    gnames;
     int       nr;
     real      tau_min;
     int       nstcmin;
     int       i, j, k, restnm;
-    bool      bExcl, bTable, bAnneal, bRest;
+    bool      bExcl, bTable, bAnneal;
     char      warn_buf[STRLEN];
 
     if (bVerbose)
@@ -3655,7 +3732,7 @@ void do_index(const char*                   mdparin,
             {
                 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
                 {
-                    sprintf(warnbuf,
+                    sprintf(warn_buf,
                             "Please use Y(ES) or N(O) for freezedim only "
                             "(not %s)",
                             freezeDims[k].c_str());
@@ -3678,15 +3755,13 @@ void do_index(const char*                   mdparin,
     add_wall_energrps(groups, ir->nwall, symtab);
     ir->opts.ngener    = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
     auto vcmGroupNames = gmx::splitString(is->vcm);
-    bRest              = do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
-                         SimulationAtomGroupType::MassCenterVelocityRemoval, restnm,
-                         vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
-    if (bRest)
+    do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
+                 SimulationAtomGroupType::MassCenterVelocityRemoval, restnm,
+                 vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
+
+    if (ir->comm_mode != ecmNO)
     {
-        warning(wi,
-                "Some atoms are not part of any center of mass motion removal group.\n"
-                "This may lead to artifacts.\n"
-                "In most cases one should use one group for the whole system.");
+        checkAndUpdateVcmFreezeGroupConsistency(groups, natoms, ir->opts, wi);
     }
 
     /* Now we have filled the freeze struct, so we can calculate NRDF */