Use more containers and views for constraints code

Lifetime and usage is clearer if we don't use raw pointers.

Refs #2423

Change-Id: I7d76e30acfd209baf2256efbd3c40290ed38caab

diff --git a/src/gromacs/domdec/domdec_constraints.cpp b/src/gromacs/domdec/domdec_constraints.cpp
index 0fe94bfc951f5da24216edb4e871bbca89d8edb3..fffb4028f8e1c1bd8b64d7162483e6505e8db18f 100644 (file)
--- a/src/gromacs/domdec/domdec_constraints.cpp
+++ b/src/gromacs/domdec/domdec_constraints.cpp
@@ -330,7 +330,7 @@ static void atoms_to_settles(gmx_domdec_t *dd,
 static void atoms_to_constraints(gmx_domdec_t *dd,
                                  const gmx_mtop_t *mtop,
                                  const int *cginfo,
-                                 const t_blocka *at2con_mt, int nrec,
+                                 gmx::ArrayRef<const t_blocka> at2con_mt, int nrec,
                                  t_ilist *ilc_local,
                                  ind_req_t *ireq)
 {
@@ -443,14 +443,14 @@ int dd_make_local_constraints(gmx_domdec_t *dd, int at_start,
                               gmx::Constraints *constr, int nrec,
                               t_ilist *il_local)
 {
-    gmx_domdec_constraints_t   *dc;
-    t_ilist                    *ilc_local, *ils_local;
-    ind_req_t                  *ireq;
-    const t_blocka             *at2con_mt;
-    const int                 **at2settle_mt;
-    gmx_hash_t                 *ga2la_specat;
+    gmx_domdec_constraints_t     *dc;
+    t_ilist                      *ilc_local, *ils_local;
+    ind_req_t                    *ireq;
+    gmx::ArrayRef<const t_blocka> at2con_mt;
+    const int                   **at2settle_mt;
+    gmx_hash_t                   *ga2la_specat;
     int at_end, i, j;
-    t_iatom                    *iap;
+    t_iatom                      *iap;
 
     // This code should not be called unless this condition is true,
     // because that's the only time init_domdec_constraints is
@@ -482,7 +482,7 @@ int dd_make_local_constraints(gmx_domdec_t *dd, int at_start,
     else
     {
         // Currently unreachable
-        at2con_mt = nullptr;
+        at2con_mt = gmx::EmptyArrayRef();
         ireq      = nullptr;
     }
 
@@ -512,14 +512,14 @@ int dd_make_local_constraints(gmx_domdec_t *dd, int at_start,
         /* Do the constraints, if present, on the first thread.
          * Do the settles on all other threads.
          */
-        t0_set = ((at2con_mt != nullptr && dc->nthread > 1) ? 1 : 0);
+        t0_set = ((!at2con_mt.empty() && dc->nthread > 1) ? 1 : 0);
 
 #pragma omp parallel for num_threads(dc->nthread) schedule(static)
         for (thread = 0; thread < dc->nthread; thread++)
         {
             try
             {
-                if (at2con_mt && thread == 0)
+                if (!at2con_mt.empty() && thread == 0)
                 {
                     atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec,
                                          ilc_local, ireq);

diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp
index ad2be05bbae50a133e44f5d14398dfdbbfa2c8a2..c8e0df3cb1bbcc11fe136ea826f29f571090bcd8 100644 (file)
--- a/src/gromacs/mdlib/constr.cpp
+++ b/src/gromacs/mdlib/constr.cpp
@@ -77,6 +77,7 @@
 #include "gromacs/topology/ifunc.h"
 #include "gromacs/topology/mtop_lookup.h"
 #include "gromacs/topology/mtop_util.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/exceptions.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/gmxassert.h"
@@ -127,33 +128,31 @@ class Constraints::Impl
                    tensor               *vir,
                    ConstraintVariable    econq);
         //! The total number of constraints.
-        int                ncon_tot = 0;
+        int                   ncon_tot = 0;
         //! The number of flexible constraints.
-        int                nflexcon = 0;
-        //! The size of at2con = number of moltypes.
-        int                n_at2con_mt = 0;
-        //! A list of atoms to constraints.
-        t_blocka          *at2con_mt = nullptr;
+        int                   nflexcon = 0;
+        //! A list of atoms to constraints for each moleculetype.
+        std::vector<t_blocka> at2con_mt;
         //! The size of at2settle = number of moltypes
-        int                n_at2settle_mt = 0;
+        int                   n_at2settle_mt = 0;
         //! A list of atoms to settles.
-        int              **at2settle_mt = nullptr;
+        int                 **at2settle_mt = nullptr;
         //! Whether any SETTLES cross charge-group boundaries.
-        bool               bInterCGsettles = false;
+        bool                  bInterCGsettles = false;
         //! LINCS data.
-        Lincs             *lincsd = nullptr;
+        Lincs                *lincsd = nullptr;
         //! SHAKE data.
-        shakedata         *shaked = nullptr;
+        shakedata            *shaked = nullptr;
         //! SETTLE data.
-        settledata        *settled = nullptr;
+        settledata           *settled = nullptr;
         //! The maximum number of warnings.
-        int                maxwarn = 0;
+        int                   maxwarn = 0;
         //! The number of warnings for LINCS.
-        int                warncount_lincs = 0;
+        int                   warncount_lincs = 0;
         //! The number of warnings for SETTLE.
-        int                warncount_settle = 0;
+        int                   warncount_settle = 0;
         //! The essential dynamics data.
-        gmx_edsam_t        ed = nullptr;
+        gmx_edsam_t           ed = nullptr;
 
         //! Thread-local virial contribution.
         tensor            *vir_r_m_dr_th = {0};
@@ -685,15 +684,15 @@ Constraints::Impl::apply(bool                  bLog,
     return bOK;
 }
 
-real *Constraints::rmsdData() const
+ArrayRef<real> Constraints::rmsdData() const
 {
     if (impl_->lincsd)
     {
-        return lincs_rmsd_data(impl_->lincsd);
+        return lincs_rmsdData(impl_->lincsd);
     }
     else
     {
-        return nullptr;
+        return EmptyArrayRef();
     }
 }
 
@@ -886,6 +885,26 @@ Constraints::setConstraints(const gmx_localtop_t &top,
     impl_->setConstraints(top, md);
 }
 
+/*! \brief Makes a per-moleculetype container of mappings from atom
+ * indices to constraint indices.
+ *
+ * Note that flexible constraints are only enabled with a dynamical integrator. */
+static std::vector<t_blocka>
+makeAtomToConstraintMappings(const gmx_mtop_t            &mtop,
+                             FlexibleConstraintTreatment  flexibleConstraintTreatment)
+{
+    std::vector<t_blocka> mapping;
+    mapping.reserve(mtop.moltype.size());
+    for (const gmx_moltype_t &moltype : mtop.moltype)
+    {
+        mapping.push_back(make_at2con(moltype.atoms.nr,
+                                      moltype.ilist,
+                                      mtop.ffparams.iparams,
+                                      flexibleConstraintTreatment));
+    }
+    return mapping;
+}
+
 Constraints::Constraints(const gmx_mtop_t     &mtop,
                          const t_inputrec     &ir,
                          FILE                 *log,
@@ -941,15 +960,8 @@ Constraints::Impl::Impl(const gmx_mtop_t     &mtop_p,
     nflexcon = 0;
     if (numConstraints > 0)
     {
-        n_at2con_mt = mtop.moltype.size();
-        snew(at2con_mt, n_at2con_mt);
-        for (int mt = 0; mt < static_cast<int>(mtop.moltype.size()); mt++)
-        {
-            at2con_mt[mt] = make_at2con(mtop.moltype[mt].atoms.nr,
-                                        mtop.moltype[mt].ilist,
-                                        mtop.ffparams.iparams,
-                                        flexibleConstraintTreatment(EI_DYNAMICS(ir_p.eI)));
-        }
+        at2con_mt = makeAtomToConstraintMappings(mtop,
+                                                 flexibleConstraintTreatment(EI_DYNAMICS(ir.eI)));
 
         for (const gmx_molblock_t &molblock : mtop.molblock)
         {
@@ -1067,7 +1079,8 @@ void Constraints::saveEdsamPointer(gmx_edsam_t ed)
     impl_->ed = ed;
 }
 
-const t_blocka *Constraints::atom2constraints_moltype() const
+const ArrayRef<const t_blocka>
+Constraints::atom2constraints_moltype() const
 {
     return impl_->at2con_mt;
 }

diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h
index a8c830d3ef9441f6562d37943c42a7c441fea702..9cb5efbe9cd286679406579cb42657bdaca7a03c 100644 (file)
--- a/src/gromacs/mdlib/constr.h
+++ b/src/gromacs/mdlib/constr.h
@@ -50,6 +50,7 @@
 
 #include "gromacs/math/vectypes.h"
 #include "gromacs/topology/idef.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/basedefinitions.h"
 #include "gromacs/utility/classhelpers.h"
 #include "gromacs/utility/gmxassert.h"
@@ -175,14 +176,14 @@ class Constraints
         //! Links the essentialdynamics and constraint code.
         void saveEdsamPointer(gmx_edsam *ed);
         //! Getter for use by domain decomposition.
-        const t_blocka *atom2constraints_moltype() const;
+        const ArrayRef<const t_blocka> atom2constraints_moltype() const;
         //! Getter for use by domain decomposition.
         const int **atom2settle_moltype() const;
 
         /*! \brief Return the data for reduction for determining
-         * constraint RMS relative deviations, or nullptr when not
-         * supported for any active constraints. */
-        real *rmsdData() const;
+         * constraint RMS relative deviations, or an empty ArrayRef
+         * when not supported for any active constraints. */
+        ArrayRef<real> rmsdData() const;
         /*! \brief Return the RMSD of the constraints when available. */
         real rmsd() const;
 

diff --git a/src/gromacs/mdlib/lincs.cpp b/src/gromacs/mdlib/lincs.cpp
index 0023d195162dcd285df60d15039ed703619af7aa..0cc682d139a30f687335761150fc3a89f6f5d128 100644 (file)
--- a/src/gromacs/mdlib/lincs.cpp
+++ b/src/gromacs/mdlib/lincs.cpp
@@ -74,6 +74,7 @@
 #include "gromacs/simd/vector_operations.h"
 #include "gromacs/topology/block.h"
 #include "gromacs/topology/mtop_util.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/basedefinitions.h"
 #include "gromacs/utility/bitmask.h"
 #include "gromacs/utility/cstringutil.h"
@@ -208,9 +209,9 @@ class Lincs
         real           *tmp4;
         /*! @} */
         //! The Lagrange multipliers times -1.
-        real           *mlambda;
+        real               *mlambda;
         //! Storage for the constraint RMS relative deviation output.
-        real            rmsd_data[3];
+        std::array<real, 2> rmsdData;
 };
 
 /*! \brief Define simd_width for memory allocation used for SIMD code */
@@ -224,16 +225,16 @@ static const int simd_width = 1;
    AlignedAllocator, which currently forces 128 byte alignment. */
 static const int align_bytes = 128;
 
-real *lincs_rmsd_data(Lincs *lincsd)
+ArrayRef<real> lincs_rmsdData(Lincs *lincsd)
 {
-    return lincsd->rmsd_data;
+    return lincsd->rmsdData;
 }
 
 real lincs_rmsd(const Lincs *lincsd)
 {
-    if (lincsd->rmsd_data[0] > 0)
+    if (lincsd->rmsdData[0] > 0)
     {
-        return std::sqrt(lincsd->rmsd_data[1]/lincsd->rmsd_data[0]);
+        return std::sqrt(lincsd->rmsdData[1]/lincsd->rmsdData[0]);
     }
     else
     {
@@ -1385,7 +1386,7 @@ static void set_lincs_matrix(Lincs *li, real *invmass, real lambda)
 
 //! Finds all triangles of atoms that share constraints to a central atom.
 static int count_triangle_constraints(const t_ilist  *ilist,
-                                      const t_blocka *at2con)
+                                      const t_blocka &at2con)
 {
     int      ncon1, ncon_tot;
     int      c0, a00, a01, n1, c1, a10, a11, ac1, n2, c2, a20, a21;
@@ -1406,9 +1407,9 @@ static int count_triangle_constraints(const t_ilist  *ilist,
         iap       = constr_iatomptr(ncon1, ia1, ia2, c0);
         a00       = iap[1];
         a01       = iap[2];
-        for (n1 = at2con->index[a01]; n1 < at2con->index[a01+1]; n1++)
+        for (n1 = at2con.index[a01]; n1 < at2con.index[a01+1]; n1++)
         {
-            c1 = at2con->a[n1];
+            c1 = at2con.a[n1];
             if (c1 != c0)
             {
                 iap = constr_iatomptr(ncon1, ia1, ia2, c1);
@@ -1422,9 +1423,9 @@ static int count_triangle_constraints(const t_ilist  *ilist,
                 {
                     ac1 = a10;
                 }
-                for (n2 = at2con->index[ac1]; n2 < at2con->index[ac1+1]; n2++)
+                for (n2 = at2con.index[ac1]; n2 < at2con.index[ac1+1]; n2++)
                 {
-                    c2 = at2con->a[n2];
+                    c2 = at2con.a[n2];
                     if (c2 != c0 && c2 != c1)
                     {
                         iap = constr_iatomptr(ncon1, ia1, ia2, c2);
@@ -1449,7 +1450,7 @@ static int count_triangle_constraints(const t_ilist  *ilist,
 
 //! Finds sequences of sequential constraints.
 static bool more_than_two_sequential_constraints(const t_ilist  *ilist,
-                                                 const t_blocka *at2con)
+                                                 const t_blocka &at2con)
 {
     t_iatom  *ia1, *ia2, *iap;
     int       ncon1, ncon_tot, c;
@@ -1469,8 +1470,8 @@ static bool more_than_two_sequential_constraints(const t_ilist  *ilist,
         a1  = iap[1];
         a2  = iap[2];
         /* Check if this constraint has constraints connected at both atoms */
-        if (at2con->index[a1+1] - at2con->index[a1] > 1 &&
-            at2con->index[a2+1] - at2con->index[a2] > 1)
+        if (at2con.index[a1+1] - at2con.index[a1] > 1 &&
+            at2con.index[a2+1] - at2con.index[a2] > 1)
         {
             bMoreThanTwoSequentialConstraints = TRUE;
         }
@@ -1486,7 +1487,7 @@ static int int_comp(const void *a, const void *b)
 }
 
 Lincs *init_lincs(FILE *fplog, const gmx_mtop_t &mtop,
-                  int nflexcon_global, const t_blocka *at2con,
+                  int nflexcon_global, ArrayRef<const t_blocka> at2con,
                   bool bPLINCS, int nIter, int nProjOrder)
 {
     Lincs                *li;
@@ -1526,10 +1527,10 @@ Lincs *init_lincs(FILE *fplog, const gmx_mtop_t &mtop,
 
         li->ncg_triangle +=
             molb.nmol*
-            count_triangle_constraints(molt.ilist, &at2con[molb.type]);
+            count_triangle_constraints(molt.ilist, at2con[molb.type]);
 
         if (!bMoreThanTwoSeq &&
-            more_than_two_sequential_constraints(molt.ilist, &at2con[molb.type]))
+            more_than_two_sequential_constraints(molt.ilist, at2con[molb.type]))
         {
             bMoreThanTwoSeq = TRUE;
         }
@@ -2434,8 +2435,7 @@ bool constrain_lincs(FILE *fplog, bool bLog, bool bEner,
     {
         if (bLog || bEner)
         {
-            lincsd->rmsd_data[0] = 0;
-            lincsd->rmsd_data[1] = 0;
+            lincsd->rmsdData = {{0}};
         }
 
         return bOK;
@@ -2530,14 +2530,12 @@ bool constrain_lincs(FILE *fplog, bool bLog, bool bEner,
         {
             cconerr(lincsd, xprime, pbc,
                     &ncons_loc, &p_ssd, &p_max, &p_imax);
-            lincsd->rmsd_data[0] = ncons_loc;
-            lincsd->rmsd_data[1] = p_ssd;
+            lincsd->rmsdData[0] = ncons_loc;
+            lincsd->rmsdData[1] = p_ssd;
         }
         else
         {
-            lincsd->rmsd_data[0] = 0;
-            lincsd->rmsd_data[1] = 0;
-            lincsd->rmsd_data[2] = 0;
+            lincsd->rmsdData = {{0}};
         }
         if (bLog && fplog && lincsd->nc > 0)
         {

diff --git a/src/gromacs/mdlib/lincs.h b/src/gromacs/mdlib/lincs.h
index 173c5587dcbea39f516c12c88a1fbddf778b43d7..da813df8482977fd793f2b35c6cc0071a7b189e9 100644 (file)
--- a/src/gromacs/mdlib/lincs.h
+++ b/src/gromacs/mdlib/lincs.h
@@ -47,6 +47,7 @@
 #include <cstdio>
 
 #include "gromacs/math/vectypes.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/basedefinitions.h"
 #include "gromacs/utility/real.h"
 
@@ -69,14 +70,14 @@ enum class ConstraintVariable : int;
 class Lincs;
 
 /*! \brief Return the data for determining constraint RMS relative deviations. */
-real *lincs_rmsd_data(Lincs *lincsd);
+ArrayRef<real> lincs_rmsdData(Lincs *lincsd);
 
 /*! \brief Return the RMSD of the constraint. */
 real lincs_rmsd(const Lincs *lincsd);
 
 /*! \brief Initializes and returns the lincs data struct. */
 Lincs *init_lincs(FILE *fplog, const gmx_mtop_t &mtop,
-                  int nflexcon_global, const t_blocka *at2con,
+                  int nflexcon_global, ArrayRef<const t_blocka> at2con,
                   bool bPLINCS, int nIter, int nProjOrder);
 
 /*! \brief Initialize lincs stuff */

diff --git a/src/gromacs/mdlib/rbin.cpp b/src/gromacs/mdlib/rbin.cpp
index c6354f22f1a837adc9f9971b680be62ce97aa9f3..1e0a209e61112b85d3b5f192560168bdebeaeef4 100644 (file)
--- a/src/gromacs/mdlib/rbin.cpp
+++ b/src/gromacs/mdlib/rbin.cpp
@@ -40,6 +40,7 @@
 #include "rbin.h"
 
 #include "gromacs/gmxlib/network.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/smalloc.h"
 
 t_bin *mk_bin(void)
@@ -66,7 +67,7 @@ void reset_bin(t_bin *b)
     b->nreal = 0;
 }
 
-int add_binr(t_bin *b, int nr, real r[])
+int add_binr(t_bin *b, int nr, const real r[])
 {
 #define MULT 4
     int     i, rest, index;
@@ -99,6 +100,11 @@ int add_binr(t_bin *b, int nr, real r[])
     return index;
 }
 
+int add_binr(t_bin *b, gmx::ArrayRef<const real> r)
+{
+    return add_binr(b, r.size(), r.data());
+}
+
 int add_bind(t_bin *b, int nr, double r[])
 {
 #define MULT 4
@@ -151,6 +157,11 @@ void extract_binr(t_bin *b, int index, int nr, real r[])
     }
 }
 
+void extract_binr(t_bin *b, int index, gmx::ArrayRef<real> r)
+{
+    extract_binr(b, index, r.size(), r.data());
+}
+
 void extract_bind(t_bin *b, int index, int nr, double r[])
 {
     int     i;

diff --git a/src/gromacs/mdlib/rbin.h b/src/gromacs/mdlib/rbin.h
index 86b0bbd5b686ce1ea8b3e7633713f75cd64ab51c..7e3a42f66c61183dd5a4d81d249e34f42ea02b3c 100644 (file)
--- a/src/gromacs/mdlib/rbin.h
+++ b/src/gromacs/mdlib/rbin.h
@@ -37,6 +37,7 @@
 #ifndef GMX_MDLIB_RBIN_H
 #define GMX_MDLIB_RBIN_H
 
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/real.h"
 
 struct t_commrec;
@@ -56,7 +57,8 @@ void destroy_bin(t_bin *b);
 void reset_bin(t_bin *b);
 /* Reset number of entries to zero */
 
-int add_binr(t_bin *b, int nr, real r[]);
+int add_binr(t_bin *b, int nr, const real r[]);
+int add_binr(t_bin *b, gmx::ArrayRef<const real> r);
 int add_bind(t_bin *b, int nr, double r[]);
 /* Add reals to the bin. Returns index */
 
@@ -64,6 +66,7 @@ void sum_bin(t_bin *b, const t_commrec *cr);
 /* Globally sum the reals in the bin */
 
 void extract_binr(t_bin *b, int index, int nr, real r[]);
+void extract_binr(t_bin *b, int index, gmx::ArrayRef<real> r);
 void extract_bind(t_bin *b, int index, int nr, double r[]);
 /* Extract values from the bin, starting from index (see add_bin) */
 

diff --git a/src/gromacs/mdlib/stat.cpp b/src/gromacs/mdlib/stat.cpp
index 478a013831f481c97fad796690fb6a468db25bf0..a84a151015f7e19d86c8bbe7a3b2edd36d2f1ff8 100644 (file)
--- a/src/gromacs/mdlib/stat.cpp
+++ b/src/gromacs/mdlib/stat.cpp
@@ -155,7 +155,6 @@ void global_stat(gmx_global_stat_t gs,
     int        inn[egNR];
     real       copyenerd[F_NRE];
     int        nener, j;
-    real      *rmsd_data = nullptr;
     double     nb;
     gmx_bool   bVV, bTemp, bEner, bPres, bConstrVir, bEkinAveVel, bReadEkin;
     bool       checkNumberOfBondedInteractions = flags & CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS;
@@ -227,16 +226,16 @@ void global_stat(gmx_global_stat_t gs,
         ifv = add_binr(rb, DIM*DIM, fvir[0]);
     }
 
-
+    gmx::ArrayRef<real> rmsdData;
     if (bEner)
     {
         ie  = add_binr(rb, nener, copyenerd);
         if (constr)
         {
-            rmsd_data = constr->rmsdData();
-            if (rmsd_data)
+            rmsdData = constr->rmsdData();
+            if (!rmsdData.empty())
             {
-                irmsd = add_binr(rb, 2, rmsd_data);
+                irmsd = add_binr(rb, 2, rmsdData.data());
             }
         }
         if (!inputrecNeedMutot(inputrec))
@@ -331,9 +330,9 @@ void global_stat(gmx_global_stat_t gs,
     if (bEner)
     {
         extract_binr(rb, ie, nener, copyenerd);
-        if (rmsd_data)
+        if (!rmsdData.empty())
         {
-            extract_binr(rb, irmsd, 2, rmsd_data);
+            extract_binr(rb, irmsd, rmsdData);
         }
         if (!inputrecNeedMutot(inputrec))
         {