Merge branch release-2021 into merge-2021-into-master
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / readir.cpp
index 55756cec59aceafdb68a2bf6d3fc13c5ed0dc791..54008d57b91a7e56c63621fdbc9e502c6c4e874b 100644 (file)
@@ -96,6 +96,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
@@ -3268,8 +3270,8 @@ static bool do_egp_flag(t_inputrec* ir, SimulationGroups* groups, const char* op
         j = 0;
         while ((j < nr)
                && gmx_strcasecmp(
-                          names[2 * i].c_str(),
-                          *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
+                       names[2 * i].c_str(),
+                       *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
         {
             j++;
         }
@@ -3280,8 +3282,8 @@ static bool do_egp_flag(t_inputrec* ir, SimulationGroups* groups, const char* op
         k = 0;
         while ((k < nr)
                && gmx_strcasecmp(
-                          names[2 * i + 1].c_str(),
-                          *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
+                       names[2 * i + 1].c_str(),
+                       *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
         {
             k++;
         }
@@ -4112,82 +4114,85 @@ 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;
-    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 */
-    for (const auto ilist : IListRange(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 };
+
+    for (const auto ilists : IListRange(sys))
     {
-        if (ilist.nmol() > 0 && (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
+        const auto& posResList   = ilists.list()[F_POSRES];
+        const auto& fbPosResList = ilists.list()[F_FBPOSRES];
+        if (ilists.nmol() > 0 && (!havePosres[XX] || !havePosres[YY] || !havePosres[ZZ]))
         {
-            for (int i = 0; i < ilist.list()[F_POSRES].size(); i += 2)
+            for (int i = 0; i < posResList.size(); i += 2)
             {
-                pr = &sys.ffparams.iparams[ilist.list()[F_POSRES].iatoms[i]];
-                for (d = 0; d < DIM; d++)
+                const t_iparams& pr = sys.ffparams.iparams[posResList.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.list()[F_FBPOSRES].size(); i += 2)
+            for (int i = 0; i < fbPosResList.size(); i += 2)
             {
                 /* Check for flat-bottom posres */
-                pr = &sys.ffparams.iparams[ilist.list()[F_FBPOSRES].iatoms[i]];
-                if (pr->fbposres.k != 0)
+                const t_iparams& pr = sys.ffparams.iparams[fbPosResList.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);
+                                      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,
@@ -4212,7 +4217,7 @@ static void check_combination_rule_differences(const gmx_mtop_t& mtop,
     ptr = getenv("GMX_LJCOMB_TOL");
     if (ptr != nullptr)
     {
-        double dbl;
+        double            dbl;
         double gmx_unused canary;
 
         if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
@@ -4320,6 +4325,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
@@ -4329,12 +4339,11 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_
     int                       i, m, c, nmol;
     bool                      bCharge;
     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 "
@@ -4431,7 +4440,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 != IntegrationAlgorithm::BD
         && ir->comm_mode == ComRemovalAlgorithm::No
-        && !(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 "
@@ -4463,11 +4473,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 != PressureCoupling::No && ir->refcoord_scaling == RefCoordScaling::No)
     {
-        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, "
@@ -4544,10 +4554,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 "