Fix grompp position restraint reference check

The check for position restraints and absolute reference always
produced a warning, also only with position restraints.

Fixes #3996

diff --git a/docs/release-notes/2021/2021.2.rst b/docs/release-notes/2021/2021.2.rst
index b8993d0aa15ceced21d20d45ab1f2f9a37dd6389..cc202689ed77f8ba31f2956220de67ed83b4f2cf 100644 (file)
--- a/docs/release-notes/2021/2021.2.rst
+++ b/docs/release-notes/2021/2021.2.rst
@@ -19,6 +19,14 @@ Fixes where mdrun could behave incorrectly
 Fixes for ``gmx`` tools
 ^^^^^^^^^^^^^^^^^^^^^^^
 
+Fix grompp check for position restraints with absolute reference
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+Fixed that grompp with position restraints would always issue a warning about
+using an absolute reference, even when an absolute reference was not used.
+
+issue:`3996`
+
 Fix error when using VMD plugin
 """""""""""""""""""""""""""""""
 

diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp
index 9f4d2fa678f25f47a043ff983b4740fc5e6012b8..2ca9db05869697bb49dd33d1649bc8b9e1bcb3f2 100644 (file)
--- a/src/gromacs/gmxpreprocess/readir.cpp
+++ b/src/gromacs/gmxpreprocess/readir.cpp
@@ -93,6 +93,8 @@
 #define MAXPTR 254
 #define NOGID 255
 
+using gmx::BasicVector;
+
 /* Resource parameters
  * Do not change any of these until you read the instruction
  * in readinp.h. Some cpp's do not take spaces after the backslash
@@ -3995,84 +3997,84 @@ static void check_disre(const gmx_mtop_t* mtop)
     }
 }
 
-static bool absolute_reference(const t_inputrec* ir, const gmx_mtop_t* sys, const bool posres_only, ivec AbsRef)
+//! Returns whether dimensions have an absolute reference due to walls, pbc or freezing
+static BasicVector<bool> haveAbsoluteReference(const t_inputrec& ir)
 {
-    int                  d, g, i;
-    gmx_mtop_ilistloop_t iloop;
-    int                  nmol;
-    const t_iparams*     pr;
-
-    clear_ivec(AbsRef);
+    BasicVector<bool> absRef = { false, false, false };
 
-    if (!posres_only)
+    /* Check the degrees of freedom of the COM (not taking COMM removal into account) */
+    for (int d = 0; d < DIM; d++)
     {
-        /* Check the COM */
-        for (d = 0; d < DIM; d++)
-        {
-            AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
-        }
-        /* Check for freeze groups */
-        for (g = 0; g < ir->opts.ngfrz; g++)
+        absRef[d] = (d >= ndof_com(&ir));
+    }
+    /* Check for freeze groups */
+    for (int g = 0; g < ir.opts.ngfrz; g++)
+    {
+        for (int d = 0; d < DIM; d++)
         {
-            for (d = 0; d < DIM; d++)
+            if (ir.opts.nFreeze[g][d] != 0)
             {
-                if (ir->opts.nFreeze[g][d] != 0)
-                {
-                    AbsRef[d] = 1;
-                }
+                absRef[d] = true;
             }
         }
     }
 
-    /* Check for position restraints */
-    iloop = gmx_mtop_ilistloop_init(sys);
+    return absRef;
+}
+
+//! Returns whether position restraints are used for dimensions
+static BasicVector<bool> havePositionRestraints(const gmx_mtop_t& sys)
+{
+    BasicVector<bool> havePosres = { false, false, false };
+
+    gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(&sys);
+    int                  nmol;
     while (const InteractionLists* ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
     {
-        if (nmol > 0 && (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
+        if (nmol > 0 && (!havePosres[XX] || !havePosres[YY] || !havePosres[ZZ]))
         {
-            for (i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
+            for (int i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
             {
-                pr = &sys->ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
-                for (d = 0; d < DIM; d++)
+                const t_iparams& pr = sys.ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
+                for (int d = 0; d < DIM; d++)
                 {
-                    if (pr->posres.fcA[d] != 0)
+                    if (pr.posres.fcA[d] != 0)
                     {
-                        AbsRef[d] = 1;
+                        havePosres[d] = true;
                     }
                 }
             }
-            for (i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
+            for (int i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
             {
                 /* Check for flat-bottom posres */
-                pr = &sys->ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
-                if (pr->fbposres.k != 0)
+                const t_iparams& pr = sys.ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
+                if (pr.fbposres.k != 0)
                 {
-                    switch (pr->fbposres.geom)
+                    switch (pr.fbposres.geom)
                     {
-                        case efbposresSPHERE: AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1; break;
-                        case efbposresCYLINDERX: AbsRef[YY] = AbsRef[ZZ] = 1; break;
-                        case efbposresCYLINDERY: AbsRef[XX] = AbsRef[ZZ] = 1; break;
+                        case efbposresSPHERE: havePosres = { true, true, true }; break;
+                        case efbposresCYLINDERX: havePosres[YY] = havePosres[ZZ] = true; break;
+                        case efbposresCYLINDERY: havePosres[XX] = havePosres[ZZ] = true; break;
                         case efbposresCYLINDER:
                         /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
-                        case efbposresCYLINDERZ: AbsRef[XX] = AbsRef[YY] = 1; break;
+                        case efbposresCYLINDERZ: havePosres[XX] = havePosres[YY] = true; break;
                         case efbposresX: /* d=XX */
                         case efbposresY: /* d=YY */
                         case efbposresZ: /* d=ZZ */
-                            d         = pr->fbposres.geom - efbposresX;
-                            AbsRef[d] = 1;
+                            havePosres[pr.fbposres.geom - efbposresX] = true;
                             break;
                         default:
                             gmx_fatal(FARGS,
-                                      " Invalid geometry for flat-bottom position restraint.\n"
+                                      "Invalid geometry for flat-bottom position restraint.\n"
                                       "Expected nr between 1 and %d. Found %d\n",
-                                      efbposresNR - 1, pr->fbposres.geom);
+                                      efbposresNR - 1, pr.fbposres.geom);
                     }
                 }
             }
         }
     }
 
-    return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
+    return havePosres;
 }
 
 static void check_combination_rule_differences(const gmx_mtop_t* mtop,
@@ -4205,6 +4207,11 @@ static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t* mtop
     }
 }
 
+static bool allTrue(const BasicVector<bool>& boolVector)
+{
+    return boolVector[0] && boolVector[1] && boolVector[2];
+}
+
 void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
 {
     // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
@@ -4216,12 +4223,11 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_
     real *                    mgrp, mt;
     rvec                      acc;
     gmx_mtop_atomloop_block_t aloopb;
-    ivec                      AbsRef;
     char                      warn_buf[STRLEN];
 
     set_warning_line(wi, mdparin, -1);
 
-    if (absolute_reference(ir, sys, false, AbsRef))
+    if (allTrue(haveAbsoluteReference(*ir)) && allTrue(havePositionRestraints(*sys)))
     {
         warning_note(wi,
                      "Removing center of mass motion in the presence of position restraints might "
@@ -4307,7 +4313,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
-        && !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) && !ETC_ANDERSEN(ir->etc))
+        && !(allTrue(haveAbsoluteReference(*ir)) || allTrue(havePositionRestraints(*sys)) || ir->nsteps <= 10)
+        && !ETC_ANDERSEN(ir->etc))
     {
         warning(wi,
                 "You are not using center of mass motion removal (mdp option comm-mode), numerical "
@@ -4335,11 +4342,11 @@ 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 != epcNO && ir->refcoord_scaling == erscNO)
     {
-        absolute_reference(ir, sys, TRUE, AbsRef);
+        const BasicVector<bool> havePosres = havePositionRestraints(*sys);
         {
             for (m = 0; m < DIM; m++)
             {
-                if (AbsRef[m] && norm2(ir->compress[m]) > 0)
+                if (havePosres[m] && norm2(ir->compress[m]) > 0)
                 {
                     warning(wi,
                             "You are using pressure coupling with absolute position restraints, "
@@ -4466,10 +4473,11 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_
         {
             if (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0)
             {
-                absolute_reference(ir, sys, FALSE, AbsRef);
+                const auto absRef     = haveAbsoluteReference(*ir);
+                const auto havePosres = havePositionRestraints(*sys);
                 for (m = 0; m < DIM; m++)
                 {
-                    if (ir->pull->coord[i].dim[m] && !AbsRef[m])
+                    if (ir->pull->coord[i].dim[m] && !(absRef[m] || havePosres[m]))
                     {
                         warning(wi,
                                 "You are using an absolute reference for pulling, but the rest of "