Clean up make_at2con

Added an enum to control the flexible constraint treatment.
Removed start (was always zero), simplified indexing and added
documentation.

Change-Id: I08494dd553e85f8ff4aa39303899639583fdc7ac

diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp
index cab849e7c094b55615e4dfb0fc7d3b539e57aff9..abebccec3c011a96d0621ffd69df144e6f60e244 100644 (file)
--- a/src/gromacs/gmxpreprocess/grompp.cpp
+++ b/src/gromacs/gmxpreprocess/grompp.cpp
@@ -1385,9 +1385,10 @@ static bool haveDecoupledModeInMol(const gmx_moltype_t *molt,
 
     const t_atom * atom = molt->atoms.atom;
 
-    t_blocka       atomToConstraints = gmx::make_at2con(0, molt->atoms.nr,
-                                                        molt->ilist, iparams,
-                                                        FALSE);
+    t_blocka       atomToConstraints =
+        gmx::make_at2con(molt->atoms.nr,
+                         molt->ilist, iparams,
+                         gmx::FlexibleConstraintTreatment::Exclude);
 
     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 a17ed8ac31c5aae0487a614e0690bc05feb7b656..9448b953b541c65158ab7dfd715dd9309626ed7c 100644 (file)
--- a/src/gromacs/mdlib/constr.cpp
+++ b/src/gromacs/mdlib/constr.cpp
@@ -79,6 +79,7 @@
 #include "gromacs/topology/mtop_util.h"
 #include "gromacs/utility/exceptions.h"
 #include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/pleasecite.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/txtdump.h"
@@ -708,76 +709,82 @@ real Constraints::rmsd() const
     }
 }
 
-t_blocka make_at2con(int start, int natoms,
-                     const t_ilist *ilist, const t_iparams *iparams,
-                     bool bDynamics)
+FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegrator)
 {
-    int      *count;
-    t_blocka  at2con;
+    if (haveDynamicsIntegrator)
+    {
+        return FlexibleConstraintTreatment::Include;
+    }
+    else
+    {
+        return FlexibleConstraintTreatment::Exclude;
+    }
+}
+
+t_blocka make_at2con(int                          numAtoms,
+                     const t_ilist               *ilists,
+                     const t_iparams             *iparams,
+                     FlexibleConstraintTreatment  flexibleConstraintTreatment)
+{
+    GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || iparams != nullptr, "With flexible constraint detection we need valid iparams");
+
+    std::vector<int> count(numAtoms);
 
-    snew(count, natoms);
     for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
     {
-        int      ncon    = ilist[ftype].nr/3;
-        t_iatom *ia      = ilist[ftype].iatoms;
-        for (int con = 0; con < ncon; con++)
+        const t_ilist &ilist  = ilists[ftype];
+        const int      stride = 1 + NRAL(ftype);
+        for (int i = 0; i < ilist.nr; i += stride)
         {
-            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)
+            if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include ||
+                !isConstraintFlexible(iparams, ilist.iatoms[i]))
             {
-                for (int i = 1; i < 3; i++)
+                for (int j = 1; j < 3; j++)
                 {
-                    int a = ia[i] - start;
+                    int a = ilist.iatoms[i + j];
                     count[a]++;
                 }
             }
-            ia += 3;
         }
     }
 
-    at2con.nr           = natoms;
-    at2con.nalloc_index = at2con.nr+1;
+    t_blocka at2con;
+    at2con.nr           = numAtoms;
+    at2con.nalloc_index = at2con.nr + 1;
     snew(at2con.index, at2con.nalloc_index);
     at2con.index[0] = 0;
-    for (int a = 0; a < natoms; a++)
+    for (int a = 0; a < numAtoms; a++)
     {
-        at2con.index[a+1] = at2con.index[a] + count[a];
-        count[a]          = 0;
+        at2con.index[a + 1] = at2con.index[a] + count[a];
+        count[a]            = 0;
     }
-    at2con.nra      = at2con.index[natoms];
+    at2con.nra      = at2con.index[at2con.nr];
     at2con.nalloc_a = at2con.nra;
     snew(at2con.a, at2con.nalloc_a);
 
     /* The F_CONSTRNC constraints have constraint numbers
      * that continue after the last F_CONSTR constraint.
      */
-    int con_tot = 0;
+    int numConstraints = 0;
     for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
     {
-        int      ncon = ilist[ftype].nr/3;
-        t_iatom *ia   = ilist[ftype].iatoms;
-        for (int con = 0; con < ncon; con++)
+        const t_ilist &ilist  = ilists[ftype];
+        const int      stride = 1 + NRAL(ftype);
+        for (int i = 0; i < ilist.nr; i += stride)
         {
-            bool bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
-                             iparams[ia[0]].constr.dB == 0);
-            if (bDynamics || !bFlexCon)
+            if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include ||
+                !isConstraintFlexible(iparams, ilist.iatoms[i]))
             {
-                for (int i = 1; i < 3; i++)
+                for (int j = 1; j < 3; j++)
                 {
-                    int a = ia[i] - start;
-                    at2con.a[at2con.index[a]+count[a]++] = con_tot;
+                    int a = ilist.iatoms[i + j];
+                    at2con.a[at2con.index[a] + count[a]++] = numConstraints;
                 }
             }
-            con_tot++;
-            ia += 3;
+            numConstraints++;
         }
     }
 
-    sfree(count);
-
     return at2con;
 }
 
@@ -938,10 +945,10 @@ 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++)
         {
-            at2con_mt[mt] = make_at2con(0, mtop.moltype[mt].atoms.nr,
+            at2con_mt[mt] = make_at2con(mtop.moltype[mt].atoms.nr,
                                         mtop.moltype[mt].ilist,
                                         mtop.ffparams.iparams,
-                                        EI_DYNAMICS(ir.eI));
+                                        flexibleConstraintTreatment(EI_DYNAMICS(ir_p.eI)));
         }
 
         for (const gmx_molblock_t &molblock : mtop.molblock)

diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h
index 7b331bad08a69e56aad36c35969bc3d9f39e1cf3..a8c830d3ef9441f6562d37943c42a7c441fea702 100644 (file)
--- a/src/gromacs/mdlib/constr.h
+++ b/src/gromacs/mdlib/constr.h
@@ -49,8 +49,10 @@
 #include <cstdio>
 
 #include "gromacs/math/vectypes.h"
+#include "gromacs/topology/idef.h"
 #include "gromacs/utility/basedefinitions.h"
 #include "gromacs/utility/classhelpers.h"
+#include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/real.h"
 
 struct gmx_edsam;
@@ -65,7 +67,6 @@ struct t_inputrec;
 struct t_mdatoms;
 struct t_nrnb;
 struct t_pbc;
-union t_iparams;
 class t_state;
 
 namespace gmx
@@ -195,16 +196,52 @@ class Constraints
 /*! \brief Generate a fatal error because of too many LINCS/SETTLE warnings. */
 void too_many_constraint_warnings(int eConstrAlg, int warncount);
 
+/*! \brief Returns whether constraint with parameter \p iparamsIndex is a flexible constraint */
+static inline bool isConstraintFlexible(const t_iparams *iparams,
+                                        int              iparamsIndex)
+{
+    GMX_ASSERT(iparams != nullptr, "Need a valid iparams array");
+
+    return (iparams[iparamsIndex].constr.dA == 0 &&
+            iparams[iparamsIndex].constr.dB == 0);
+};
+
 /* The at2con t_blocka struct returned by the routines below
  * contains a list of constraints per atom.
  * The F_CONSTRNC constraints in this structure number consecutively
  * after the F_CONSTR constraints.
  */
 
-/*! \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);
+/*! \brief Tells make_at2con how to treat flexible constraints */
+enum class FlexibleConstraintTreatment
+{
+    Include, //!< Include all flexible constraints
+    Exclude  //!< Exclude all flexible constraints
+};
+
+/*! \brief Returns the flexible constraint treatment depending on whether the integrator is dynamic */
+FlexibleConstraintTreatment
+flexibleConstraintTreatment(bool haveDynamicsIntegrator);
+
+/*! \brief Returns a block struct to go from atoms to constraints
+ *
+ * The block struct will contain constraint indices with lower indices
+ * directly matching the order in F_CONSTR and higher indices matching
+ * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR.
+ *
+ * \param[in]  numAtoms  The number of atoms to construct the list for
+ * \param[in]  ilist     Interaction list, size F_NRE
+ * \param[in]  iparams   Interaction parameters, can be null when flexibleConstraintTreatment=Include
+ * \param[in]  flexibleConstraintTreatment  The flexible constraint treatment, see enum above
+ * \returns a block struct with all constraints for each atom
+ */
+t_blocka make_at2con(int                          numAtoms,
+                     const t_ilist               *ilist,
+                     const t_iparams             *iparams,
+                     FlexibleConstraintTreatment  flexibleConstraintTreatment);
+
+/*! \brief Returns an array of atom to constraints lists for the moltypes */
+const t_blocka *atom2constraints_moltype(const Constraints *constr);
 
 //! Return the number of flexible constraints in the \c ilist and \c iparams.
 int countFlexibleConstraints(const t_ilist   *ilist,

diff --git a/src/gromacs/mdlib/constraintrange.cpp b/src/gromacs/mdlib/constraintrange.cpp
index 381c5a449fb726d8bf4da68a61bd4c482a84db8f..95cf473a8ce2acad822eb309f33424c37ad641dc 100644 (file)
--- a/src/gromacs/mdlib/constraintrange.cpp
+++ b/src/gromacs/mdlib/constraintrange.cpp
@@ -170,8 +170,8 @@ 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));
+    at2con = make_at2con(natoms, molt->ilist, iparams,
+                         flexibleConstraintTreatment(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 ae096002478c97abb131840021ccab10e4c472fe..0023d195162dcd285df60d15039ed703619af7aa 100644 (file)
--- a/src/gromacs/mdlib/lincs.cpp
+++ b/src/gromacs/mdlib/lincs.cpp
@@ -2024,7 +2024,9 @@ void set_lincs(const t_idef         &idef,
     {
         natoms = md.homenr;
     }
-    at2con = make_at2con(0, natoms, idef.il, idef.iparams, bDynamics);
+
+    at2con = make_at2con(natoms, idef.il, idef.iparams,
+                         flexibleConstraintTreatment(bDynamics));
 
     ncon_tot = idef.il[F_CONSTR].nr/3;