Separate flexible constraint counting

Functions should do one thing, particularly when the thing is only
needed some of the time and complicates returning multiple things.

Minimized scope of some variables, avoiding risky re-use of variables
with the same name.

Refs #2423

Change-Id: I1860447b4780e76c20a53964a2a946512689b872

diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp
index a4fe12a05b9dce491adca27e66a52c3e8c7bddf2..cab849e7c094b55615e4dfb0fc7d3b539e57aff9 100644 (file)
--- a/src/gromacs/gmxpreprocess/grompp.cpp
+++ b/src/gromacs/gmxpreprocess/grompp.cpp
@@ -1385,10 +1385,9 @@ static bool haveDecoupledModeInMol(const gmx_moltype_t *molt,
 
     const t_atom * atom = molt->atoms.atom;
 
-    int            numFlexibleConstraints;
     t_blocka       atomToConstraints = gmx::make_at2con(0, molt->atoms.nr,
                                                         molt->ilist, iparams,
-                                                        FALSE, &numFlexibleConstraints);
+                                                        FALSE);
 
     bool           haveDecoupledMode = false;
     for (int ftype = 0; ftype < F_NRE; ftype++)

diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp
index 7e7c956f582193f4992d6b8428c3103f319767c4..a17ed8ac31c5aae0487a614e0690bc05feb7b656 100644 (file)
--- a/src/gromacs/mdlib/constr.cpp
+++ b/src/gromacs/mdlib/constr.cpp
@@ -710,45 +710,39 @@ real Constraints::rmsd() const
 
 t_blocka make_at2con(int start, int natoms,
                      const t_ilist *ilist, const t_iparams *iparams,
-                     bool bDynamics, int *nflexiblecons)
+                     bool bDynamics)
 {
-    int      *count, ncon, con, con_tot, nflexcon, ftype, i, a;
-    t_iatom  *ia;
+    int      *count;
     t_blocka  at2con;
-    bool      bFlexCon;
 
     snew(count, natoms);
-    nflexcon = 0;
-    for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
+    for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
     {
-        ncon = ilist[ftype].nr/3;
-        ia   = ilist[ftype].iatoms;
-        for (con = 0; con < ncon; con++)
+        int      ncon    = ilist[ftype].nr/3;
+        t_iatom *ia      = ilist[ftype].iatoms;
+        for (int con = 0; con < ncon; con++)
         {
-            bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
-                        iparams[ia[0]].constr.dB == 0);
-            if (bFlexCon)
-            {
-                nflexcon++;
-            }
+            bool bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
+                             iparams[ia[0]].constr.dB == 0);
+            // Flexible constraints are only applied with dynamical
+            // integrators, not e.g. energy minimization.
             if (bDynamics || !bFlexCon)
             {
-                for (i = 1; i < 3; i++)
+                for (int i = 1; i < 3; i++)
                 {
-                    a = ia[i] - start;
+                    int a = ia[i] - start;
                     count[a]++;
                 }
             }
             ia += 3;
         }
     }
-    *nflexiblecons = nflexcon;
 
     at2con.nr           = natoms;
     at2con.nalloc_index = at2con.nr+1;
     snew(at2con.index, at2con.nalloc_index);
     at2con.index[0] = 0;
-    for (a = 0; a < natoms; a++)
+    for (int a = 0; a < natoms; a++)
     {
         at2con.index[a+1] = at2con.index[a] + count[a];
         count[a]          = 0;
@@ -760,20 +754,20 @@ t_blocka make_at2con(int start, int natoms,
     /* The F_CONSTRNC constraints have constraint numbers
      * that continue after the last F_CONSTR constraint.
      */
-    con_tot = 0;
-    for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
+    int con_tot = 0;
+    for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
     {
-        ncon = ilist[ftype].nr/3;
-        ia   = ilist[ftype].iatoms;
-        for (con = 0; con < ncon; con++)
+        int      ncon = ilist[ftype].nr/3;
+        t_iatom *ia   = ilist[ftype].iatoms;
+        for (int con = 0; con < ncon; con++)
         {
-            bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
-                        iparams[ia[0]].constr.dB == 0);
+            bool bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
+                             iparams[ia[0]].constr.dB == 0);
             if (bDynamics || !bFlexCon)
             {
-                for (i = 1; i < 3; i++)
+                for (int i = 1; i < 3; i++)
                 {
-                    a = ia[i] - start;
+                    int a = ia[i] - start;
                     at2con.a[at2con.index[a]+count[a]++] = con_tot;
                 }
             }
@@ -787,6 +781,29 @@ t_blocka make_at2con(int start, int natoms,
     return at2con;
 }
 
+int countFlexibleConstraints(const t_ilist   *ilist,
+                             const t_iparams *iparams)
+{
+    int nflexcon = 0;
+    for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
+    {
+        const int numIatomsPerConstraint = 3;
+        int       ncon                   = ilist[ftype].nr /  numIatomsPerConstraint;
+        t_iatom  *ia                     = ilist[ftype].iatoms;
+        for (int con = 0; con < ncon; con++)
+        {
+            if (iparams[ia[0]].constr.dA == 0 &&
+                iparams[ia[0]].constr.dB == 0)
+            {
+                nflexcon++;
+            }
+            ia += numIatomsPerConstraint;
+        }
+    }
+
+    return nflexcon;
+}
+
 //! Returns the index of the settle to which each atom belongs.
 static int *make_at2settle(int natoms, const t_ilist *ilist)
 {
@@ -921,18 +938,17 @@ Constraints::Impl::Impl(const gmx_mtop_t     &mtop_p,
         snew(at2con_mt, n_at2con_mt);
         for (int mt = 0; mt < static_cast<int>(mtop.moltype.size()); mt++)
         {
-            int nflexcon;
             at2con_mt[mt] = make_at2con(0, mtop.moltype[mt].atoms.nr,
                                         mtop.moltype[mt].ilist,
                                         mtop.ffparams.iparams,
-                                        EI_DYNAMICS(ir.eI), &nflexcon);
-            for (const gmx_molblock_t &molblock : mtop.molblock)
-            {
-                if (molblock.type == mt)
-                {
-                    nflexcon += molblock.nmol*nflexcon;
-                }
-            }
+                                        EI_DYNAMICS(ir.eI));
+        }
+
+        for (const gmx_molblock_t &molblock : mtop.molblock)
+        {
+            int count = countFlexibleConstraints(mtop.moltype[molblock.type].ilist,
+                                                 mtop.ffparams.iparams);
+            nflexcon += molblock.nmol*count;
         }
 
         if (nflexcon > 0)

diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h
index 7abbe8ce1706acaf74c443d111f69572375dba19..7b331bad08a69e56aad36c35969bc3d9f39e1cf3 100644 (file)
--- a/src/gromacs/mdlib/constr.h
+++ b/src/gromacs/mdlib/constr.h
@@ -204,7 +204,11 @@ void too_many_constraint_warnings(int eConstrAlg, int warncount);
 /*! \brief Returns a block struct to go from atoms to constraints */
 t_blocka make_at2con(int start, int natoms,
                      const t_ilist *ilist, const t_iparams *iparams,
-                     bool bDynamics, int *nflexiblecons);
+                     bool bDynamics);
+
+//! Return the number of flexible constraints in the \c ilist and \c iparams.
+int countFlexibleConstraints(const t_ilist   *ilist,
+                             const t_iparams *iparams);
 
 /*! \brief Macro for getting the constraint iatoms for a constraint number con
  * which comes from a list where F_CONSTR and F_CONSTRNC constraints

diff --git a/src/gromacs/mdlib/constraintrange.cpp b/src/gromacs/mdlib/constraintrange.cpp
index 2e573bbe02a21ea7eebb90e07f32ef197ca77df6..381c5a449fb726d8bf4da68a61bd4c482a84db8f 100644 (file)
--- a/src/gromacs/mdlib/constraintrange.cpp
+++ b/src/gromacs/mdlib/constraintrange.cpp
@@ -157,7 +157,7 @@ static real constr_r_max_moltype(const gmx_moltype_t *molt,
                                  const t_iparams     *iparams,
                                  const t_inputrec    *ir)
 {
-    int      natoms, nflexcon, *path, at, count;
+    int      natoms, *path, at, count;
 
     t_blocka at2con;
     real     r0, r1, r2maxA, r2maxB, rmax, lam0, lam1;
@@ -171,7 +171,7 @@ static real constr_r_max_moltype(const gmx_moltype_t *molt,
     natoms = molt->atoms.nr;
 
     at2con = make_at2con(0, natoms, molt->ilist, iparams,
-                         EI_DYNAMICS(ir->eI), &nflexcon);
+                         EI_DYNAMICS(ir->eI));
     snew(path, 1+ir->nProjOrder);
     for (at = 0; at < 1+ir->nProjOrder; at++)
     {

diff --git a/src/gromacs/mdlib/lincs.cpp b/src/gromacs/mdlib/lincs.cpp
index e471a98292844d26359fad0fc2699f6a3d330e2e..ae096002478c97abb131840021ccab10e4c472fe 100644 (file)
--- a/src/gromacs/mdlib/lincs.cpp
+++ b/src/gromacs/mdlib/lincs.cpp
@@ -1971,7 +1971,7 @@ void set_lincs(const t_idef         &idef,
                const t_commrec      *cr,
                Lincs                *li)
 {
-    int          natoms, nflexcon;
+    int          natoms;
     t_blocka     at2con;
     t_iatom     *iatom;
     int          i, ncc_alloc_old, ncon_tot;
@@ -2024,8 +2024,7 @@ void set_lincs(const t_idef         &idef,
     {
         natoms = md.homenr;
     }
-    at2con = make_at2con(0, natoms, idef.il, idef.iparams, bDynamics,
-                         &nflexcon);
+    at2con = make_at2con(0, natoms, idef.il, idef.iparams, bDynamics);
 
     ncon_tot = idef.il[F_CONSTR].nr/3;
 
@@ -2073,7 +2072,7 @@ void set_lincs(const t_idef         &idef,
         /* With energy minimization, flexible constraints are ignored
          * (and thus minimized, as they should be).
          */
-        ncon_assign -= nflexcon;
+        ncon_assign -= countFlexibleConstraints(idef.il, idef.iparams);
     }
 
     /* Set the target constraint count per task to exactly uniform,