Use ArrayRef in signatures of constraints and some pulling
authorJoe Jordan <ejjordan12@gmail.com>
Thu, 25 Mar 2021 06:11:17 +0000 (06:11 +0000)
committerMark Abraham <mark.j.abraham@gmail.com>
Thu, 25 Mar 2021 06:11:17 +0000 (06:11 +0000)
In order to use ArrayRef in constr, settle, shake, and lincs, some
pull code also needed to be refactored. Part of preparation for
refactoring of t_mdatoms.

22 files changed:
src/gromacs/applied_forces/awh/awh.cpp
src/gromacs/applied_forces/awh/awh.h
src/gromacs/domdec/mdsetup.cpp
src/gromacs/gmxpreprocess/grompp.cpp
src/gromacs/gmxpreprocess/readir.h
src/gromacs/gmxpreprocess/readpull.cpp
src/gromacs/mdlib/constr.cpp
src/gromacs/mdlib/constr.h
src/gromacs/mdlib/lincs.cpp
src/gromacs/mdlib/lincs.h
src/gromacs/mdlib/settle.cpp
src/gromacs/mdlib/settle.h
src/gromacs/mdlib/shake.cpp
src/gromacs/mdlib/shake.h
src/gromacs/mdlib/sim_util.cpp
src/gromacs/mdlib/tests/constrtestrunners.cpp
src/gromacs/mdlib/tests/settletestrunners.cpp
src/gromacs/mdlib/tests/shake.cpp
src/gromacs/mdrun/md.cpp
src/gromacs/pulling/pull.cpp
src/gromacs/pulling/pull.h
src/gromacs/pulling/pullutil.cpp

index d987b09f2ceabe2131174fdbfce88eb8c8907ff8..4a5e5aeeb4c335c34f1ec6923513257670ce74ad 100644 (file)
@@ -289,7 +289,7 @@ bool Awh::isOutputStep(int64_t step) const
 }
 
 real Awh::applyBiasForcesAndUpdateBias(PbcType                pbcType,
-                                       const real*            masses,
+                                       ArrayRef<const real>   masses,
                                        ArrayRef<const double> neighborLambdaEnergies,
                                        ArrayRef<const double> neighborLambdaDhdl,
                                        const matrix           box,
index c1d3632b0214078b71f7bfda1dd0889be37ffa4a..752953739495318be731fb2529a1f94f6313a4f3 100644 (file)
@@ -186,7 +186,7 @@ public:
      * \returns the potential energy for the bias.
      */
     real applyBiasForcesAndUpdateBias(PbcType                pbcType,
-                                      const real*            masses,
+                                      ArrayRef<const real>   masses,
                                       ArrayRef<const double> neighborLambdaEnergies,
                                       ArrayRef<const double> neighborLambdaDhdl,
                                       const matrix           box,
index ccf793b0bd7088f25841b0422862ae151c61c746..0baae937e79f05204e5e1410c47b64e49f3b78a1 100644 (file)
@@ -52,6 +52,7 @@
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/topology/mtop_util.h"
 #include "gromacs/topology/topology.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
 
@@ -105,7 +106,7 @@ void mdAlgorithmsSetupAtomData(const t_commrec*     cr,
              numHomeAtoms,
              mdAtoms);
 
-    auto mdatoms = mdAtoms->mdatoms();
+    t_mdatoms* mdatoms = mdAtoms->mdatoms();
     if (usingDomDec)
     {
         dd_sort_local_top(*cr->dd, mdatoms, top);
@@ -150,11 +151,12 @@ void mdAlgorithmsSetupAtomData(const t_commrec*     cr,
         constr->setConstraints(top,
                                mdatoms->nr,
                                mdatoms->homenr,
-                               mdatoms->massT,
-                               mdatoms->invmass,
+                               gmx::arrayRefFromArray(mdatoms->massT, mdatoms->nr),
+                               gmx::arrayRefFromArray(mdatoms->invmass, mdatoms->nr),
                                mdatoms->nMassPerturbed != 0,
                                mdatoms->lambda,
-                               mdatoms->cFREEZE);
+                               mdatoms->cFREEZE ? gmx::arrayRefFromArray(mdatoms->cFREEZE, mdatoms->nr)
+                                                : gmx::ArrayRef<const unsigned short>());
     }
 }
 
index a6d17f12b00899c4ffa78d123697e3149696ae13..85f6e376cc52a6cf14ee207ec476816ada8116a1 100644 (file)
@@ -2391,12 +2391,8 @@ int gmx_grompp(int argc, char* argv[])
 
     if (ir->bPull)
     {
-        pull = set_pull_init(ir,
-                             sys,
-                             state.x.rvec_array(),
-                             state.box,
-                             state.lambda[FreeEnergyPerturbationCouplingType::Mass],
-                             wi);
+        pull = set_pull_init(
+                ir, sys, state.x, state.box, state.lambda[FreeEnergyPerturbationCouplingType::Mass], wi);
     }
 
     /* Modules that supply external potential for pull coordinates
index dfec8859a6c1a117c9ff012113afce84b286627b..da326072f561a153e057750f08ba40cad8a29059 100644 (file)
@@ -165,7 +165,12 @@ void checkPullCoords(gmx::ArrayRef<const t_pull_group> pullGroups,
                      gmx::ArrayRef<const t_pull_coord> pullCoords);
 /* Process the pull coordinates after reading the pull groups */
 
-pull_t* set_pull_init(t_inputrec* ir, const gmx_mtop_t& mtop, rvec* x, matrix box, real lambda, warninp_t wi);
+pull_t* set_pull_init(t_inputrec*                    ir,
+                      const gmx_mtop_t&              mtop,
+                      gmx::ArrayRef<const gmx::RVec> x,
+                      matrix                         box,
+                      real                           lambda,
+                      warninp_t                      wi);
 /* Prints the initial pull group distances in x.
  * If requested, adds the current distance to the initial reference location.
  * Returns the pull_t pull work struct. This should be passed to finish_pull()
index 5af3b4b9e6c93cd5a3eb36f20ce13b6d246b3b48..2f8e417c8ab121995c78c25c120a406215a416e2 100644 (file)
@@ -549,7 +549,12 @@ void checkPullCoords(gmx::ArrayRef<const t_pull_group> pullGroups, gmx::ArrayRef
     }
 }
 
-pull_t* set_pull_init(t_inputrec* ir, const gmx_mtop_t& mtop, rvec* x, matrix box, real lambda, warninp_t wi)
+pull_t* set_pull_init(t_inputrec*                    ir,
+                      const gmx_mtop_t&              mtop,
+                      gmx::ArrayRef<const gmx::RVec> x,
+                      matrix                         box,
+                      real                           lambda,
+                      warninp_t                      wi)
 {
     pull_t* pull_work;
     t_pbc   pbc;
@@ -573,18 +578,14 @@ pull_t* set_pull_init(t_inputrec* ir, const gmx_mtop_t& mtop, rvec* x, matrix bo
 
     if (pull->bSetPbcRefToPrevStepCOM)
     {
-        initPullComFromPrevStep(nullptr, pull_work, md->massT, &pbc, x);
+        initPullComFromPrevStep(nullptr, pull_work, gmx::arrayRefFromArray(md->massT, md->nr), &pbc, x);
     }
-    pull_calc_coms(nullptr, pull_work, md->massT, &pbc, t_start, x, nullptr);
+    pull_calc_coms(nullptr, pull_work, gmx::arrayRefFromArray(md->massT, md->nr), &pbc, t_start, x, {});
 
     for (int g = 0; g < pull->ngroup; g++)
     {
-        bool groupObeysPbc = pullCheckPbcWithinGroup(
-                *pull_work,
-                gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec*>(x), mtop.natoms),
-                pbc,
-                g,
-                c_pullGroupSmallGroupThreshold);
+        bool groupObeysPbc =
+                pullCheckPbcWithinGroup(*pull_work, x, pbc, g, c_pullGroupSmallGroupThreshold);
         if (!groupObeysPbc)
         {
             char buf[STRLEN];
@@ -616,12 +617,7 @@ pull_t* set_pull_init(t_inputrec* ir, const gmx_mtop_t& mtop, rvec* x, matrix bo
         }
         if (groupObeysPbc)
         {
-            groupObeysPbc = pullCheckPbcWithinGroup(
-                    *pull_work,
-                    gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec*>(x), mtop.natoms),
-                    pbc,
-                    g,
-                    c_pullGroupPbcMargin);
+            groupObeysPbc = pullCheckPbcWithinGroup(*pull_work, x, pbc, g, c_pullGroupPbcMargin);
             if (!groupObeysPbc)
             {
                 char buf[STRLEN];
index 1cb575ad6a3bf809f88b59aee73248da5aabd047..d8cddba26ed7dede05b2f0fe9f7a0136b3e364e7 100644 (file)
@@ -59,6 +59,7 @@
 #include "gromacs/fileio/gmxfio.h"
 #include "gromacs/fileio/pdbio.h"
 #include "gromacs/gmxlib/nrnb.h"
+#include "gromacs/math/arrayrefwithpadding.h"
 #include "gromacs/math/utilities.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/gmx_omp_nthreads.h"
@@ -75,6 +76,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"
@@ -109,14 +111,14 @@ public:
          int                   numConstraints,
          int                   numSettles);
     ~Impl();
-    void setConstraints(gmx_localtop_t* top,
-                        int             numAtoms,
-                        int             numHomeAtoms,
-                        real*           masses,
-                        real*           inverseMasses,
-                        bool            hasMassPerturbedAtoms,
-                        real            lambda,
-                        unsigned short* cFREEZE);
+    void setConstraints(gmx_localtop_t*                     top,
+                        int                                 numAtoms,
+                        int                                 numHomeAtoms,
+                        gmx::ArrayRef<const real>           masses,
+                        gmx::ArrayRef<const real>           inverseMasses,
+                        bool                                hasMassPerturbedAtoms,
+                        real                                lambda,
+                        gmx::ArrayRef<const unsigned short> cFREEZE);
     bool apply(bool                      bLog,
                bool                      bEner,
                int64_t                   step,
@@ -169,15 +171,15 @@ public:
     //! Number of local atoms.
     int numHomeAtoms_ = 0;
     //! Atoms masses.
-    real* masses_;
+    gmx::ArrayRef<const real> masses_;
     //! Inverse masses.
-    real* inverseMasses_;
+    gmx::ArrayRef<const real> inverseMasses_;
     //! If there are atoms with perturbed mass (for FEP).
     bool hasMassPerturbedAtoms_;
     //! FEP lambda value.
     real lambda_;
     //! Freeze groups data
-    unsigned short* cFREEZE_;
+    gmx::ArrayRef<const unsigned short> cFREEZE_;
     //! Whether we need to do pbc for handling bonds.
     bool pbcHandlingRequired_ = false;
 
@@ -768,9 +770,9 @@ bool Constraints::Impl::apply(bool                      bLog,
                             cr,
                             ir.delta_t,
                             t,
-                            as_rvec_array(x.unpaddedArrayRef().data()),
-                            as_rvec_array(xprime.unpaddedArrayRef().data()),
-                            as_rvec_array(v.unpaddedArrayRef().data()),
+                            x.unpaddedArrayRef(),
+                            xprime.unpaddedArrayRef(),
+                            v.unpaddedArrayRef(),
                             constraintsVirial);
         }
         if (ed && delta_step > 0)
@@ -788,7 +790,7 @@ bool Constraints::Impl::apply(bool                      bLog,
     wallcycle_stop(wcycle, ewcCONSTR);
 
     const bool haveVelocities = (!v.empty() || econq == ConstraintVariable::Velocities);
-    if (haveVelocities && cFREEZE_)
+    if (haveVelocities && !cFREEZE_.empty())
     {
         /* Set the velocities of frozen dimensions to zero */
         ArrayRef<RVec> vRef;
@@ -993,14 +995,14 @@ static std::vector<int> make_at2settle(int natoms, const InteractionList& ilist)
     return at2s;
 }
 
-void Constraints::Impl::setConstraints(gmx_localtop_t* top,
-                                       int             numAtoms,
-                                       int             numHomeAtoms,
-                                       real*           masses,
-                                       real*           inverseMasses,
-                                       const bool      hasMassPerturbedAtoms,
-                                       const real      lambda,
-                                       unsigned short* cFREEZE)
+void Constraints::Impl::setConstraints(gmx_localtop_t*                     top,
+                                       int                                 numAtoms,
+                                       int                                 numHomeAtoms,
+                                       gmx::ArrayRef<const real>           masses,
+                                       gmx::ArrayRef<const real>           inverseMasses,
+                                       const bool                          hasMassPerturbedAtoms,
+                                       const real                          lambda,
+                                       gmx::ArrayRef<const unsigned short> cFREEZE)
 {
     numAtoms_              = numAtoms;
     numHomeAtoms_          = numHomeAtoms;
@@ -1050,14 +1052,14 @@ void Constraints::Impl::setConstraints(gmx_localtop_t* top,
     }
 }
 
-void Constraints::setConstraints(gmx_localtop_t* top,
-                                 const int       numAtoms,
-                                 const int       numHomeAtoms,
-                                 real*           masses,
-                                 real*           inverseMasses,
-                                 const bool      hasMassPerturbedAtoms,
-                                 const real      lambda,
-                                 unsigned short* cFREEZE)
+void Constraints::setConstraints(gmx_localtop_t*                     top,
+                                 const int                           numAtoms,
+                                 const int                           numHomeAtoms,
+                                 gmx::ArrayRef<const real>           masses,
+                                 gmx::ArrayRef<const real>           inverseMasses,
+                                 const bool                          hasMassPerturbedAtoms,
+                                 const real                          lambda,
+                                 gmx::ArrayRef<const unsigned short> cFREEZE)
 {
     impl_->setConstraints(
             top, numAtoms, numHomeAtoms, masses, inverseMasses, hasMassPerturbedAtoms, lambda, cFREEZE);
index eb51ffb602e7a74743099f5137664064c16809dc..03eb4772ac0b2f581117f0a01b3a0d9d0e24999e 100644 (file)
@@ -132,14 +132,14 @@ public:
      *
      * \todo Make this a callback that is called automatically
      * once a new domain has been made. */
-    void setConstraints(gmx_localtop_t* top,
-                        int             numAtoms,
-                        int             numHomeAtoms,
-                        real*           masses,
-                        real*           inverseMasses,
-                        bool            hasMassPerturbedAtoms,
-                        real            lambda,
-                        unsigned short* cFREEZE);
+    void setConstraints(gmx_localtop_t*                     top,
+                        int                                 numAtoms,
+                        int                                 numHomeAtoms,
+                        gmx::ArrayRef<const real>           masses,
+                        gmx::ArrayRef<const real>           inverseMasses,
+                        bool                                hasMassPerturbedAtoms,
+                        real                                lambda,
+                        gmx::ArrayRef<const unsigned short> cFREEZE);
 
     /*! \brief Applies constraints to coordinates.
      *
index b5fe54be92ee94ebe18c38b6b53cc707de496279..1bbfd8330e955c636c3c5dceb2e4b535f954f83d 100644 (file)
@@ -349,10 +349,10 @@ static void lincs_update_atoms_noind(int                            ncons,
                                      real                           preFactor,
                                      gmx::ArrayRef<const real>      fac,
                                      gmx::ArrayRef<const gmx::RVec> r,
-                                     const real*                    invmass,
+                                     gmx::ArrayRef<const real>      invmass,
                                      rvec*                          x)
 {
-    if (invmass != nullptr)
+    if (!invmass.empty())
     {
         for (int b = 0; b < ncons; b++)
         {
@@ -398,10 +398,10 @@ static void lincs_update_atoms_ind(gmx::ArrayRef<const int>       ind,
                                    real                           preFactor,
                                    gmx::ArrayRef<const real>      fac,
                                    gmx::ArrayRef<const gmx::RVec> r,
-                                   const real*                    invmass,
+                                   gmx::ArrayRef<const real>      invmass,
                                    rvec*                          x)
 {
-    if (invmass != nullptr)
+    if (!invmass.empty())
     {
         for (int b : ind)
         {
@@ -447,7 +447,7 @@ static void lincs_update_atoms(Lincs*                         li,
                                real                           preFactor,
                                gmx::ArrayRef<const real>      fac,
                                gmx::ArrayRef<const gmx::RVec> r,
-                               const real*                    invmass,
+                               gmx::ArrayRef<const real>      invmass,
                                rvec*                          x)
 {
     if (li->ntask == 1)
@@ -574,7 +574,7 @@ static void do_lincsp(ArrayRefWithPadding<const RVec> xPadded,
                       t_pbc*                          pbc,
                       Lincs*                          lincsd,
                       int                             th,
-                      const real*                     invmass,
+                      ArrayRef<const real>            invmass,
                       ConstraintVariable              econq,
                       bool                            bCalcDHDL,
                       bool                            bCalcVir,
@@ -713,7 +713,7 @@ static void do_lincsp(ArrayRefWithPadding<const RVec> xPadded,
                        1.0,
                        sol,
                        r,
-                       (econq != ConstraintVariable::Force) ? invmass : nullptr,
+                       (econq != ConstraintVariable::Force) ? invmass : gmx::ArrayRef<real>(),
                        as_rvec_array(fp.data()));
 
     if (bCalcDHDL)
@@ -958,7 +958,7 @@ static void do_lincs(ArrayRefWithPadding<const RVec> xPadded,
                      t_pbc*                          pbc,
                      Lincs*                          lincsd,
                      int                             th,
-                     const real*                     invmass,
+                     ArrayRef<const real>            invmass,
                      const t_commrec*                cr,
                      bool                            bCalcDHDL,
                      real                            wangle,
@@ -1217,7 +1217,11 @@ static void do_lincs(ArrayRefWithPadding<const RVec> xPadded,
 }
 
 /*! \brief Sets the elements in the LINCS matrix for task task. */
-static void set_lincs_matrix_task(Lincs* li, Task* li_task, const real* invmass, int* ncc_triangle, int* nCrossTaskTriangles)
+static void set_lincs_matrix_task(Lincs*               li,
+                                  Task*                li_task,
+                                  ArrayRef<const real> invmass,
+                                  int*                 ncc_triangle,
+                                  int*                 nCrossTaskTriangles)
 {
     /* Construct the coupling coefficient matrix blmf */
     li_task->ntriangle   = 0;
@@ -1305,7 +1309,7 @@ static void set_lincs_matrix_task(Lincs* li, Task* li_task, const real* invmass,
 }
 
 /*! \brief Sets the elements in the LINCS matrix. */
-static void set_lincs_matrix(Lincs* li, const real* invmass, real lambda)
+static void set_lincs_matrix(Lincs* li, ArrayRef<const real> invmass, real lambda)
 {
     const real invsqrt2 = 0.7071067811865475244;
 
@@ -1854,7 +1858,7 @@ static void set_matrix_indices(Lincs* li, const Task& li_task, const ListOfLists
 
 void set_lincs(const InteractionDefinitions& idef,
                const int                     numAtoms,
-               const real*                   invmass,
+               ArrayRef<const real>          invmass,
                const real                    lambda,
                bool                          bDynamics,
                const t_commrec*              cr,
@@ -2268,7 +2272,7 @@ bool constrain_lincs(bool                            computeRmsd,
                      const t_inputrec&               ir,
                      int64_t                         step,
                      Lincs*                          lincsd,
-                     const real*                     invmass,
+                     ArrayRef<const real>            invmass,
                      const t_commrec*                cr,
                      const gmx_multisim_t*           ms,
                      ArrayRefWithPadding<const RVec> xPadded,
index 0b00370c018719e946f92c5f214b8d946579c0c2..fcaf91a9854cf9355dcda3e9d57ceff3f68552ce 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -61,9 +61,10 @@ struct t_pbc;
 
 namespace gmx
 {
-
 template<typename>
 class ArrayRefWithPadding;
+template<typename>
+class ArrayRef;
 enum class ConstraintVariable : int;
 class Lincs;
 template<typename>
@@ -91,7 +92,7 @@ void done_lincs(Lincs* li);
 /*! \brief Initialize lincs stuff */
 void set_lincs(const InteractionDefinitions& idef,
                int                           numAtoms,
-               const real*                   invmass,
+               ArrayRef<const real>          invmass,
                real                          lambda,
                bool                          bDynamics,
                const t_commrec*              cr,
@@ -104,7 +105,7 @@ bool constrain_lincs(bool                            computeRmsd,
                      const t_inputrec&               ir,
                      int64_t                         step,
                      Lincs*                          lincsd,
-                     const real*                     invmass,
+                     ArrayRef<const real>            invmass,
                      const t_commrec*                cr,
                      const gmx_multisim_t*           ms,
                      ArrayRefWithPadding<const RVec> x,
index 43e21eb486b58675d8765545f22295ff11e4474b..0ccadf81e95e1161a55798c611bc28a085655d03 100644 (file)
@@ -4,7 +4,7 @@
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -57,7 +57,6 @@
 #include "gromacs/math/invertmatrix.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/constr.h"
-#include "gromacs/mdtypes/mdatom.h"
 #include "gromacs/pbcutil/ishift.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/pbcutil/pbc_simd.h"
@@ -66,6 +65,7 @@
 #include "gromacs/topology/idef.h"
 #include "gromacs/topology/ifunc.h"
 #include "gromacs/topology/mtop_util.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/gmxassert.h"
 
@@ -196,10 +196,10 @@ SettleData::SettleData(const gmx_mtop_t& mtop) :
     parametersAllMasses1_ = settleParameters(1.0, 1.0, 1.0, 1.0, dOH, dHH);
 }
 
-void SettleData::setConstraints(const InteractionList& il_settle,
-                                const int              numHomeAtoms,
-                                const real*            masses,
-                                const real*            inverseMasses)
+void SettleData::setConstraints(const InteractionList&    il_settle,
+                                const int                 numHomeAtoms,
+                                gmx::ArrayRef<const real> masses,
+                                gmx::ArrayRef<const real> inverseMasses)
 {
 #if GMX_SIMD_HAVE_REAL
     const int pack_size = GMX_SIMD_REAL_WIDTH;
index 4714bd669fcd4b9a2d78d254ec6dc54fb1dd87fb..86bf81fbf12478561e3563ca6b515a0354c67ec1 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -121,10 +121,10 @@ public:
     SettleData(const gmx_mtop_t& mtop);
 
     //! Sets the constraints from the interaction list and the masses
-    void setConstraints(const InteractionList& il_settle,
-                        int                    numHomeAtoms,
-                        const real*            masses,
-                        const real*            inverseMasses);
+    void setConstraints(const InteractionList&    il_settle,
+                        int                       numHomeAtoms,
+                        gmx::ArrayRef<const real> masses,
+                        gmx::ArrayRef<const real> inverseMasses);
 
     //! Returns settle parameters for constraining coordinates and forces
     const SettleParameters& parametersMassWeighted() const { return parametersMassWeighted_; }
index 80e2deee7eb95e48c6c968e901efa39dc29367ac..47ef4167f233eb43f289ade1391e75879411d6d2 100644 (file)
@@ -61,6 +61,7 @@
 #include "gromacs/mdtypes/md_enums.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/topology/invblock.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/smalloc.h"
 
@@ -282,7 +283,7 @@ void cshake(const int            iatom[],
             ArrayRef<const RVec> initial_displacements,
             ArrayRef<const real> half_of_reduced_mass,
             real                 omega,
-            const real           invmass[],
+            ArrayRef<const real> invmass,
             ArrayRef<const real> distance_squared_tolerance,
             ArrayRef<real>       scaled_lagrange_multiplier,
             int*                 nerror)
@@ -371,7 +372,7 @@ static void crattle(const int            iatom[],
                     ArrayRef<const RVec> rij,
                     ArrayRef<const real> m2,
                     real                 omega,
-                    const real           invmass[],
+                    ArrayRef<const real> invmass,
                     ArrayRef<const real> distance_squared_tolerance,
                     ArrayRef<real>       scaled_lagrange_multiplier,
                     int*                 nerror,
@@ -440,7 +441,7 @@ static void crattle(const int            iatom[],
 //! Applies SHAKE
 static int vec_shakef(FILE*                     fplog,
                       shakedata*                shaked,
-                      const real                invmass[],
+                      ArrayRef<const real>      invmass,
                       int                       ncon,
                       ArrayRef<const t_iparams> ip,
                       const int*                iatom,
@@ -636,7 +637,7 @@ static void check_cons(FILE*                     log,
                        const t_pbc*              pbc,
                        ArrayRef<const t_iparams> ip,
                        const int*                iatom,
-                       const real                invmass[],
+                       ArrayRef<const real>      invmass,
                        ConstraintVariable        econq)
 {
     int  ai, aj;
@@ -699,7 +700,7 @@ static void check_cons(FILE*                     log,
 //! Applies SHAKE.
 static bool bshakef(FILE*                         log,
                     shakedata*                    shaked,
-                    const real                    invmass[],
+                    ArrayRef<const real>          invmass,
                     const InteractionDefinitions& idef,
                     const t_inputrec&             ir,
                     ArrayRef<const RVec>          x_s,
@@ -821,7 +822,7 @@ static bool bshakef(FILE*                         log,
 
 bool constrain_shake(FILE*                         log,
                      shakedata*                    shaked,
-                     const real                    invmass[],
+                     ArrayRef<const real>          invmass,
                      const InteractionDefinitions& idef,
                      const t_inputrec&             ir,
                      ArrayRef<const RVec>          x_s,
index 8559c57553855291095dedf5eadb2247512dc50c..f0723818bc7c899236e15558d90cd1c0ebdd8fa4 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -108,17 +108,17 @@ void make_shake_sblock_dd(shakedata* shaked, const InteractionList& ilcon);
  * sblock[n] to sblock[n+1]. Array sblock should be large enough.
  * Return TRUE when OK, FALSE when shake-error
  */
-bool constrain_shake(FILE*                         log,       /* Log file                      */
-                     shakedata*                    shaked,    /* Total number of atoms */
-                     const real                    invmass[], /* Atomic masses         */
-                     const InteractionDefinitions& idef,      /* The interaction def           */
-                     const t_inputrec&             ir,        /* Input record                  */
-                     ArrayRef<const RVec>          x_s,       /* Coords before update          */
-                     ArrayRef<RVec>                xprime, /* Output coords when constraining x */
-                     ArrayRef<RVec>                vprime, /* Output coords when constraining v */
-                     const t_pbc*                  pbc,    /* PBC information              */
-                     t_nrnb*                       nrnb,   /* Performance measure          */
-                     real                          lambda, /* FEP lambda                   */
+bool constrain_shake(FILE*                         log,     /* Log file                        */
+                     shakedata*                    shaked,  /* Total number of atoms   */
+                     gmx::ArrayRef<const real>     invmass, /* Atomic masses           */
+                     const InteractionDefinitions& idef,    /* The interaction def             */
+                     const t_inputrec&             ir,      /* Input record                    */
+                     ArrayRef<const RVec>          x_s,     /* Coords before update            */
+                     ArrayRef<RVec>                xprime,  /* Output coords when constraining x */
+                     ArrayRef<RVec>                vprime,  /* Output coords when constraining v */
+                     const t_pbc*                  pbc,     /* PBC information              */
+                     t_nrnb*                       nrnb,    /* Performance measure          */
+                     real                          lambda,  /* FEP lambda                   */
                      real*                         dvdlambda,  /* FEP force                    */
                      real                          invdt,      /* 1/delta_t                    */
                      ArrayRef<RVec>                v,          /* Also constrain v if not empty  */
@@ -138,7 +138,7 @@ void cshake(const int            iatom[],
             ArrayRef<const RVec> rij,
             ArrayRef<const real> half_of_reduced_mass,
             real                 omega,
-            const real           invmass[],
+            ArrayRef<const real> invmass,
             ArrayRef<const real> distance_squared_tolerance,
             ArrayRef<real>       scaled_lagrange_multiplier,
             int*                 nerror);
index 6f42594ed010fd8ffb858cc3e04fb4efdaa1310f..424e9844fc4e0e8a138b64e07ac1d978aa7b6699 100644 (file)
@@ -205,12 +205,12 @@ static void pull_potential_wrapper(const t_commrec*               cr,
     dvdl = 0;
     enerd->term[F_COM_PULL] +=
             pull_potential(pull_work,
-                           mdatoms->massT,
+                           gmx::arrayRefFromArray(mdatoms->massT, mdatoms->nr),
                            &pbc,
                            cr,
                            t,
                            lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Restraint)],
-                           as_rvec_array(x.data()),
+                           x,
                            force,
                            &dvdl);
     enerd->dvdl_lin[FreeEnergyPerturbationCouplingType::Restraint] += dvdl;
@@ -680,16 +680,17 @@ static void computeSpecialForces(FILE*                          fplog,
             }
 
             auto& forceWithVirial = (mtsLevel == 0) ? forceWithVirialMtsLevel0 : forceWithVirialMtsLevel1;
-            enerd->term[F_COM_PULL] += awh->applyBiasForcesAndUpdateBias(inputrec.pbcType,
-                                                                         mdatoms->massT,
-                                                                         foreignLambdaDeltaH,
-                                                                         foreignLambdaDhDl,
-                                                                         box,
-                                                                         forceWithVirial,
-                                                                         t,
-                                                                         step,
-                                                                         wcycle,
-                                                                         fplog);
+            enerd->term[F_COM_PULL] += awh->applyBiasForcesAndUpdateBias(
+                    inputrec.pbcType,
+                    gmx::arrayRefFromArray(mdatoms->massT, mdatoms->nr),
+                    foreignLambdaDeltaH,
+                    foreignLambdaDhDl,
+                    box,
+                    forceWithVirial,
+                    t,
+                    step,
+                    wcycle,
+                    fplog);
         }
     }
 
index d8b0bfc53c80ef3c0236671ea3d3666c2e62ce05..32ebfbeb5541b175845c4013938f785ca877c53f 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -86,7 +86,7 @@ void ShakeConstraintsRunner::applyConstraints(ConstraintsTestData* testData, t_p
     make_shake_sblock_serial(&shaked, testData->idef_.get(), testData->numAtoms_);
     bool success = constrain_shake(nullptr,
                                    &shaked,
-                                   testData->invmass_.data(),
+                                   testData->invmass_,
                                    *testData->idef_,
                                    testData->ir_,
                                    testData->x_,
@@ -141,7 +141,7 @@ void LincsConstraintsRunner::applyConstraints(ConstraintsTestData* testData, t_p
                         testData->ir_.nProjOrder);
     set_lincs(*testData->idef_,
               testData->numAtoms_,
-              testData->invmass_.data(),
+              testData->invmass_,
               testData->lambda_,
               EI_DYNAMICS(testData->ir_.eI),
               &cr,
@@ -152,7 +152,7 @@ void LincsConstraintsRunner::applyConstraints(ConstraintsTestData* testData, t_p
                                    testData->ir_,
                                    0,
                                    lincsd,
-                                   testData->invmass_.data(),
+                                   testData->invmass_,
                                    &cr,
                                    &ms,
                                    testData->x_.arrayRefWithPadding(),
index 6dc80bf66e4f6716315f7c41da73b846219d48a6..c494341cc2b836967ee24222176d6fc2d1e4a910 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -65,10 +65,8 @@ void SettleHostTestRunner::applySettle(SettleTestData*    testData,
 {
     SettleData settled(testData->mtop_);
 
-    settled.setConstraints(testData->idef_->il[F_SETTLE],
-                           testData->numAtoms_,
-                           testData->masses_.data(),
-                           testData->inverseMasses_.data());
+    settled.setConstraints(
+            testData->idef_->il[F_SETTLE], testData->numAtoms_, testData->masses_, testData->inverseMasses_);
 
     bool errorOccured;
     int  numThreads  = 1;
index 50ef3e0fb2c70960603eb9dcb431c5eb491d6d4e..7f58294382fc033278734f35bfb22a0d68a67778 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2014,2015,2017,2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2017,2018,2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -205,7 +205,7 @@ public:
                initialDisplacements,
                halfOfReducedMasses,
                omega_,
-               inverseMasses.data(),
+               inverseMasses,
                distanceSquaredTolerances,
                lagrangianValues,
                &numErrors);
index 77d3641a6db22983f2f0021d84de954a329f1b15..5584b60232d7e89090ea58cacc193c2144edd0e5 100644 (file)
@@ -515,8 +515,13 @@ void gmx::LegacySimulator::do_md()
         EnergyData::initializeEnergyHistory(startingBehavior, observablesHistory, &energyOutput);
     }
 
-    preparePrevStepPullCom(
-            ir, pull_work, mdatoms->massT, state, state_global, cr, startingBehavior != StartingBehavior::NewSimulation);
+    preparePrevStepPullCom(ir,
+                           pull_work,
+                           gmx::arrayRefFromArray(mdatoms->massT, mdatoms->nr),
+                           state,
+                           state_global,
+                           cr,
+                           startingBehavior != StartingBehavior::NewSimulation);
 
     // TODO: Remove this by converting AWH into a ForceProvider
     auto awh = prepareAwhModule(fplog,
index ae29d66801edcef995d39581182a6325b24fc219..5a5538ccb9d81d4f8387804aed9e63894b6cc11f 100644 (file)
@@ -37,7 +37,6 @@
  */
 #include "gmxpre.h"
 
-#include "gromacs/utility/stringutil.h"
 #include "pull.h"
 
 #include "config.h"
@@ -71,6 +70,7 @@
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/topology/mtop_lookup.h"
 #include "gromacs/topology/topology.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/exceptions.h"
 #include "gromacs/utility/fatalerror.h"
@@ -80,6 +80,7 @@
 #include "gromacs/utility/real.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/strconvert.h"
+#include "gromacs/utility/stringutil.h"
 
 #include "pull_internal.h"
 
@@ -175,13 +176,13 @@ double pull_conversion_factor_internal2userinput(const t_pull_coord* pcrd)
 }
 
 /* Apply forces in a mass weighted fashion for part of the pull group */
-static void apply_forces_grp_part(const pull_group_work_t* pgrp,
-                                  int                      ind_start,
-                                  int                      ind_end,
-                                  const real*              masses,
-                                  const dvec               f_pull,
-                                  int                      sign,
-                                  rvec*                    f)
+static void apply_forces_grp_part(const pull_group_work_t*  pgrp,
+                                  int                       ind_start,
+                                  int                       ind_end,
+                                  gmx::ArrayRef<const real> masses,
+                                  const dvec                f_pull,
+                                  int                       sign,
+                                  rvec*                     f)
 {
     double inv_wm = pgrp->mwscale;
 
@@ -203,12 +204,12 @@ static void apply_forces_grp_part(const pull_group_work_t* pgrp,
 }
 
 /* Apply forces in a mass weighted fashion */
-static void apply_forces_grp(const pull_group_work_t* pgrp,
-                             const real*              masses,
-                             const dvec               f_pull,
-                             int                      sign,
-                             rvec*                    f,
-                             int                      nthreads)
+static void apply_forces_grp(const pull_group_work_t*  pgrp,
+                             gmx::ArrayRef<const real> masses,
+                             const dvec                f_pull,
+                             int                       sign,
+                             rvec*                     f,
+                             int                       nthreads)
 {
     auto localAtomIndices = pgrp->atomSet.localIndex();
 
@@ -243,13 +244,13 @@ static void apply_forces_grp(const pull_group_work_t* pgrp,
 }
 
 /* Apply forces in a mass weighted fashion to a cylinder group */
-static void apply_forces_cyl_grp(const pull_group_work_t* pgrp,
-                                 const double             dv_corr,
-                                 const real*              masses,
-                                 const dvec               f_pull,
-                                 double                   f_scal,
-                                 int                      sign,
-                                 rvec*                    f,
+static void apply_forces_cyl_grp(const pull_group_work_t*  pgrp,
+                                 const double              dv_corr,
+                                 gmx::ArrayRef<const real> masses,
+                                 const dvec                f_pull,
+                                 double                    f_scal,
+                                 int                       sign,
+                                 rvec*                     f,
                                  int gmx_unused nthreads)
 {
     double inv_wm = pgrp->mwscale;
@@ -290,10 +291,10 @@ static void apply_forces_cyl_grp(const pull_group_work_t* pgrp,
 /* Apply torque forces in a mass weighted fashion to the groups that define
  * the pull vector direction for pull coordinate pcrd.
  */
-static void apply_forces_vec_torque(const struct pull_t*     pull,
-                                    const pull_coord_work_t* pcrd,
-                                    const real*              masses,
-                                    rvec*                    f)
+static void apply_forces_vec_torque(const struct pull_t*      pull,
+                                    const pull_coord_work_t*  pcrd,
+                                    gmx::ArrayRef<const real> masses,
+                                    rvec*                     f)
 {
     const PullCoordSpatialData& spatialData = pcrd->spatialData;
 
@@ -326,7 +327,7 @@ static void apply_forces_vec_torque(const struct pull_t*     pull,
 static void apply_forces_coord(struct pull_t*               pull,
                                int                          coord,
                                const PullCoordVectorForces& forces,
-                               const real*                  masses,
+                               gmx::ArrayRef<const real>    masses,
                                rvec*                        f)
 {
     /* Here it would be more efficient to use one large thread-parallel
@@ -817,8 +818,14 @@ void clear_pull_forces(pull_t* pull)
 }
 
 /* Apply constraint using SHAKE */
-static void
-do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaster, tensor vir, double dt, double t)
+static void do_constraint(struct pull_t*           pull,
+                          t_pbc*                   pbc,
+                          gmx::ArrayRef<gmx::RVec> x,
+                          gmx::ArrayRef<gmx::RVec> v,
+                          gmx_bool                 bMaster,
+                          tensor                   vir,
+                          double                   dt,
+                          double                   t)
 {
 
     dvec*    r_ij;   /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
@@ -1111,7 +1118,7 @@ do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaste
 
     /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
 
-    if (v)
+    if (!v.empty())
     {
         invdt = 1 / dt;
     }
@@ -1147,7 +1154,7 @@ do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaste
             {
                 x[ii][m] += tmp[m];
             }
-            if (v)
+            if (!v.empty())
             {
                 for (m = 0; m < DIM; m++)
                 {
@@ -1532,11 +1539,11 @@ static void check_external_potential_registration(const struct pull_t* pull)
  * potential energy is added either to the pull term or to a term
  * specific to the external module.
  */
-void apply_external_pull_coord_force(struct pull_t*        pull,
-                                     int                   coord_index,
-                                     double                coord_force,
-                                     const real*           masses,
-                                     gmx::ForceWithVirial* forceWithVirial)
+void apply_external_pull_coord_force(struct pull_t*            pull,
+                                     int                       coord_index,
+                                     double                    coord_force,
+                                     gmx::ArrayRef<const real> masses,
+                                     gmx::ForceWithVirial*     forceWithVirial)
 {
     pull_coord_work_t* pcrd;
 
@@ -1602,15 +1609,15 @@ static PullCoordVectorForces do_pull_pot_coord(struct pull_t* pull,
     return pullCoordForces;
 }
 
-real pull_potential(struct pull_t*        pull,
-                    const real*           masses,
-                    t_pbc*                pbc,
-                    const t_commrec*      cr,
-                    double                t,
-                    real                  lambda,
-                    const rvec*           x,
-                    gmx::ForceWithVirial* force,
-                    real*                 dvdlambda)
+real pull_potential(struct pull_t*                 pull,
+                    gmx::ArrayRef<const real>      masses,
+                    t_pbc*                         pbc,
+                    const t_commrec*               cr,
+                    double                         t,
+                    real                           lambda,
+                    gmx::ArrayRef<const gmx::RVec> x,
+                    gmx::ForceWithVirial*          force,
+                    real*                          dvdlambda)
 {
     real V = 0;
 
@@ -1626,7 +1633,7 @@ real pull_potential(struct pull_t*        pull,
     {
         real dVdl = 0;
 
-        pull_calc_coms(cr, pull, masses, pbc, t, x, nullptr);
+        pull_calc_coms(cr, pull, masses, pbc, t, x, {});
 
         rvec*      f             = as_rvec_array(force->force_.data());
         matrix     virial        = { { 0 } };
@@ -1666,16 +1673,16 @@ real pull_potential(struct pull_t*        pull,
     return (MASTER(cr) ? V : 0.0);
 }
 
-void pull_constraint(struct pull_t*   pull,
-                     const real*      masses,
-                     t_pbc*           pbc,
-                     const t_commrec* cr,
-                     double           dt,
-                     double           t,
-                     rvec*            x,
-                     rvec*            xp,
-                     rvec*            v,
-                     tensor           vir)
+void pull_constraint(struct pull_t*            pull,
+                     gmx::ArrayRef<const real> masses,
+                     t_pbc*                    pbc,
+                     const t_commrec*          cr,
+                     double                    dt,
+                     double                    t,
+                     gmx::ArrayRef<gmx::RVec>  x,
+                     gmx::ArrayRef<gmx::RVec>  xp,
+                     gmx::ArrayRef<gmx::RVec>  v,
+                     tensor                    vir)
 {
     assert(pull != nullptr);
 
@@ -2400,13 +2407,13 @@ static void destroy_pull(struct pull_t* pull)
     delete pull;
 }
 
-void preparePrevStepPullCom(const t_inputrec* ir,
-                            pull_t*           pull_work,
-                            const real*       masses,
-                            t_state*          state,
-                            const t_state*    state_global,
-                            const t_commrec*  cr,
-                            bool              startingFromCheckpoint)
+void preparePrevStepPullCom(const t_inputrec*         ir,
+                            pull_t*                   pull_work,
+                            gmx::ArrayRef<const real> masses,
+                            t_state*                  state,
+                            const t_state*            state_global,
+                            const t_commrec*          cr,
+                            bool                      startingFromCheckpoint)
 {
     if (!ir->pull || !ir->pull->bSetPbcRefToPrevStepCOM)
     {
@@ -2432,7 +2439,8 @@ void preparePrevStepPullCom(const t_inputrec* ir,
     {
         t_pbc pbc;
         set_pbc(&pbc, ir->pbcType, state->box);
-        initPullComFromPrevStep(cr, pull_work, masses, &pbc, state->x.rvec_array());
+        initPullComFromPrevStep(
+                cr, pull_work, masses, &pbc, state->x.arrayRefWithPadding().unpaddedArrayRef());
         updatePrevStepPullCom(pull_work, state);
     }
 }
index a91784cf835d8e96915e2d18b1460afc7b8b57fa..4016bed9faeaa1874b2b063132f4f9a619f5ee86 100644 (file)
@@ -72,6 +72,8 @@ class t_state;
 
 namespace gmx
 {
+template<typename>
+class ArrayRef;
 class ForceWithVirial;
 class LocalAtomSetManager;
 } // namespace gmx
@@ -155,12 +157,11 @@ void register_external_pull_potential(struct pull_t* pull, int coord_index, cons
  * \param[in]     masses           Atoms masses.
  * \param[in,out] forceWithVirial  Force and virial buffers.
  */
-void apply_external_pull_coord_force(struct pull_t*        pull,
-                                     int                   coord_index,
-                                     double                coord_force,
-                                     const real*           masses,
-                                     gmx::ForceWithVirial* forceWithVirial);
-
+void apply_external_pull_coord_force(struct pull_t*            pull,
+                                     int                       coord_index,
+                                     double                    coord_force,
+                                     gmx::ArrayRef<const real> masses,
+                                     gmx::ForceWithVirial*     forceWithVirial);
 
 /*! \brief Set the all the pull forces to zero.
  *
@@ -183,15 +184,15 @@ void clear_pull_forces(struct pull_t* pull);
  *
  * \returns The pull potential energy.
  */
-real pull_potential(struct pull_t*        pull,
-                    const real*           masses,
-                    struct t_pbc*         pbc,
-                    const t_commrec*      cr,
-                    double                t,
-                    real                  lambda,
-                    const rvec*           x,
-                    gmx::ForceWithVirial* force,
-                    real*                 dvdlambda);
+real pull_potential(struct pull_t*                 pull,
+                    gmx::ArrayRef<const real>      masses,
+                    struct t_pbc*                  pbc,
+                    const t_commrec*               cr,
+                    double                         t,
+                    real                           lambda,
+                    gmx::ArrayRef<const gmx::RVec> x,
+                    gmx::ForceWithVirial*          force,
+                    real*                          dvdlambda);
 
 
 /*! \brief Constrain the coordinates xp in the directions in x
@@ -208,16 +209,16 @@ real pull_potential(struct pull_t*        pull,
  * \param[in,out] v      Velocities, which may get a pull correction.
  * \param[in,out] vir    The virial, which, if != NULL, gets a pull correction.
  */
-void pull_constraint(struct pull_t*   pull,
-                     const real*      masses,
-                     struct t_pbc*    pbc,
-                     const t_commrec* cr,
-                     double           dt,
-                     double           t,
-                     rvec*            x,
-                     rvec*            xp,
-                     rvec*            v,
-                     tensor           vir);
+void pull_constraint(struct pull_t*            pull,
+                     gmx::ArrayRef<const real> masses,
+                     struct t_pbc*             pbc,
+                     const t_commrec*          cr,
+                     double                    dt,
+                     double                    t,
+                     gmx::ArrayRef<gmx::RVec>  x,
+                     gmx::ArrayRef<gmx::RVec>  xp,
+                     gmx::ArrayRef<gmx::RVec>  v,
+                     tensor                    vir);
 
 
 /*! \brief Make a selection of the home atoms for all pull groups.
@@ -266,13 +267,13 @@ void finish_pull(struct pull_t* pull);
  * \param[in,out] xp   Updated x, can be NULL.
  *
  */
-void pull_calc_coms(const t_commrec* cr,
-                    pull_t*          pull,
-                    const real*      masses,
-                    t_pbc*           pbc,
-                    double           t,
-                    const rvec       x[],
-                    rvec*            xp);
+void pull_calc_coms(const t_commrec*               cr,
+                    pull_t*                        pull,
+                    gmx::ArrayRef<const real>      masses,
+                    t_pbc*                         pbc,
+                    double                         t,
+                    gmx::ArrayRef<const gmx::RVec> x,
+                    gmx::ArrayRef<gmx::RVec>       xp);
 
 /*! \brief Margin for checking pull group PBC distances compared to half the box size */
 static constexpr real c_pullGroupPbcMargin = 0.9;
@@ -297,7 +298,7 @@ static constexpr real c_pullGroupSmallGroupThreshold = 0.5;
  * \param[in] pbcMargin  The minimum margin (as a fraction) to half the box size
  * \returns -1 when all groups obey PBC or the first group index that fails PBC
  */
-int pullCheckPbcWithinGroups(const pull_t& pull, const rvec* x, const t_pbc& pbc, real pbcMargin);
+int pullCheckPbcWithinGroups(const pull_t& pull, gmx::ArrayRef<const gmx::RVec> x, const t_pbc& pbc, real pbcMargin);
 
 /*! \brief Checks whether a specific group that uses a reference atom is within PBC restrictions
  *
@@ -376,13 +377,13 @@ void updatePrevStepPullCom(struct pull_t* pull, t_state* state);
  * \param[in] cr                     Struct for communication info.
  * \param[in] startingFromCheckpoint Is the simulation starting from a checkpoint?
  */
-void preparePrevStepPullCom(const t_inputrec* ir,
-                            pull_t*           pull_work,
-                            const real*       masses,
-                            t_state*          state,
-                            const t_state*    state_global,
-                            const t_commrec*  cr,
-                            bool              startingFromCheckpoint);
+void preparePrevStepPullCom(const t_inputrec*         ir,
+                            pull_t*                   pull_work,
+                            gmx::ArrayRef<const real> masses,
+                            t_state*                  state,
+                            const t_state*            state_global,
+                            const t_commrec*          cr,
+                            bool                      startingFromCheckpoint);
 
 /*! \brief Initializes the COM of the previous step (set to initial COM)
  *
@@ -392,6 +393,10 @@ void preparePrevStepPullCom(const t_inputrec* ir,
  * \param[in] pbc      Information struct about periodicity.
  * \param[in] x        The local positions.
  */
-void initPullComFromPrevStep(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc* pbc, const rvec x[]);
+void initPullComFromPrevStep(const t_commrec*               cr,
+                             pull_t*                        pull,
+                             gmx::ArrayRef<const real>      masses,
+                             t_pbc*                         pbc,
+                             gmx::ArrayRef<const gmx::RVec> x);
 
 #endif
index c4985a8f80ef618dd92193f5a8ec83e7edbfc7e7..b9dbb2d68ea509637955801aa9147d6036774254 100644 (file)
 #include "gromacs/math/units.h"
 #include "gromacs/math/utilities.h"
 #include "gromacs/math/vec.h"
+#include "gromacs/math/vectypes.h"
 #include "gromacs/mdtypes/commrec.h"
 #include "gromacs/mdtypes/inputrec.h"
 #include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/mdtypes/mdatom.h"
 #include "gromacs/mdtypes/state.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/pulling/pull.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/futil.h"
 #include "gromacs/utility/gmxassert.h"
@@ -119,7 +120,7 @@ static void pullAllReduce(const t_commrec* cr, pull_comm_t* comm, int n, T* data
 /* Copies the coordinates of the PBC atom of pgrp to x_pbc.
  * When those coordinates are not available on this rank, clears x_pbc.
  */
-static void setPbcAtomCoords(const pull_group_work_t& pgrp, const rvec* x, rvec x_pbc)
+static void setPbcAtomCoords(const pull_group_work_t& pgrp, gmx::ArrayRef<const gmx::RVec> x, rvec x_pbc)
 {
     if (pgrp.pbcAtomSet != nullptr)
     {
@@ -140,7 +141,10 @@ static void setPbcAtomCoords(const pull_group_work_t& pgrp, const rvec* x, rvec
     }
 }
 
-static void pull_set_pbcatoms(const t_commrec* cr, struct pull_t* pull, const rvec* x, gmx::ArrayRef<gmx::RVec> x_pbc)
+static void pull_set_pbcatoms(const t_commrec*               cr,
+                              struct pull_t*                 pull,
+                              gmx::ArrayRef<const gmx::RVec> x,
+                              gmx::ArrayRef<gmx::RVec>       x_pbc)
 {
     int numPbcAtoms = 0;
     for (size_t g = 0; g < pull->group.size(); g++)
@@ -167,8 +171,12 @@ static void pull_set_pbcatoms(const t_commrec* cr, struct pull_t* pull, const rv
     }
 }
 
-static void
-make_cyl_refgrps(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc* pbc, double t, const rvec* x)
+static void make_cyl_refgrps(const t_commrec*               cr,
+                             pull_t*                        pull,
+                             gmx::ArrayRef<const real>      masses,
+                             t_pbc*                         pbc,
+                             double                         t,
+                             gmx::ArrayRef<const gmx::RVec> x)
 {
     pull_comm_t* comm = &pull->comm;
 
@@ -380,15 +388,15 @@ static double atan2_0_2pi(double y, double x)
     return a;
 }
 
-static void sum_com_part(const pull_group_work_t* pgrp,
-                         int                      ind_start,
-                         int                      ind_end,
-                         const rvec*              x,
-                         const rvec*              xp,
-                         const real*              mass,
-                         const t_pbc*             pbc,
-                         const rvec               x_pbc,
-                         ComSums*                 sum_com)
+static void sum_com_part(const pull_group_work_t*       pgrp,
+                         int                            ind_start,
+                         int                            ind_end,
+                         gmx::ArrayRef<const gmx::RVec> x,
+                         gmx::ArrayRef<const gmx::RVec> xp,
+                         gmx::ArrayRef<const real>      mass,
+                         const t_pbc*                   pbc,
+                         const rvec                     x_pbc,
+                         ComSums*                       sum_com)
 {
     double sum_wm   = 0;
     double sum_wwm  = 0;
@@ -421,7 +429,7 @@ static void sum_com_part(const pull_group_work_t* pgrp,
             {
                 sum_wmx[d] += wm * x[ii][d];
             }
-            if (xp)
+            if (!xp.empty())
             {
                 for (int d = 0; d < DIM; d++)
                 {
@@ -439,7 +447,7 @@ static void sum_com_part(const pull_group_work_t* pgrp,
             {
                 sum_wmx[d] += wm * dx[d];
             }
-            if (xp)
+            if (!xp.empty())
             {
                 /* For xp add the difference between xp and x to dx,
                  * such that we use the same periodic image,
@@ -456,21 +464,21 @@ static void sum_com_part(const pull_group_work_t* pgrp,
     sum_com->sum_wm  = sum_wm;
     sum_com->sum_wwm = sum_wwm;
     copy_dvec(sum_wmx, sum_com->sum_wmx);
-    if (xp)
+    if (!xp.empty())
     {
         copy_dvec(sum_wmxp, sum_com->sum_wmxp);
     }
 }
 
-static void sum_com_part_cosweight(const pull_group_work_t* pgrp,
-                                   int                      ind_start,
-                                   int                      ind_end,
-                                   int                      cosdim,
-                                   real                     twopi_box,
-                                   const rvec*              x,
-                                   const rvec*              xp,
-                                   const real*              mass,
-                                   ComSums*                 sum_com)
+static void sum_com_part_cosweight(const pull_group_work_t*       pgrp,
+                                   int                            ind_start,
+                                   int                            ind_end,
+                                   int                            cosdim,
+                                   real                           twopi_box,
+                                   gmx::ArrayRef<const gmx::RVec> x,
+                                   gmx::ArrayRef<const gmx::RVec> xp,
+                                   gmx::ArrayRef<const real>      mass,
+                                   ComSums*                       sum_com)
 {
     /* Cosine weighting geometry */
     double sum_cm  = 0;
@@ -496,7 +504,7 @@ static void sum_com_part_cosweight(const pull_group_work_t* pgrp,
         sum_csm += static_cast<double>(cw * sw * m);
         sum_ssm += static_cast<double>(sw * sw * m);
 
-        if (xp != nullptr)
+        if (!xp.empty())
         {
             real cw = std::cos(xp[ii][cosdim] * twopi_box);
             real sw = std::sin(xp[ii][cosdim] * twopi_box);
@@ -515,7 +523,13 @@ static void sum_com_part_cosweight(const pull_group_work_t* pgrp,
 }
 
 /* calculates center of mass of selection index from all coordinates x */
-void pull_calc_coms(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc* pbc, double t, const rvec x[], rvec* xp)
+void pull_calc_coms(const t_commrec*               cr,
+                    pull_t*                        pull,
+                    gmx::ArrayRef<const real>      masses,
+                    t_pbc*                         pbc,
+                    double                         t,
+                    gmx::ArrayRef<const gmx::RVec> x,
+                    gmx::ArrayRef<gmx::RVec>       xp)
 {
     real         twopi_box = 0;
     pull_comm_t* comm;
@@ -605,9 +619,9 @@ void pull_calc_coms(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc
                 if (pgrp->params.ind.size() == 1 && pgrp->atomSet.numAtomsLocal() == 1
                     && masses[pgrp->atomSet.localIndex()[0]] == 0)
                 {
-                    GMX_ASSERT(xp == nullptr,
+                    GMX_ASSERT(xp.empty(),
                                "We should not have groups with zero mass with constraints, i.e. "
-                               "xp!=NULL");
+                               "xp not empty");
 
                     /* Copy the single atom coordinate */
                     for (int d = 0; d < DIM; d++)
@@ -743,14 +757,14 @@ void pull_calc_coms(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc
                 for (m = 0; m < DIM; m++)
                 {
                     pgrp->x[m] = comBuffer[0][m] * pgrp->mwscale;
-                    if (xp)
+                    if (!xp.empty())
                     {
                         pgrp->xp[m] = comBuffer[1][m] * pgrp->mwscale;
                     }
                     if (pgrp->epgrppbc == epgrppbcREFAT || pgrp->epgrppbc == epgrppbcPREVSTEPCOM)
                     {
                         pgrp->x[m] += comm->pbcAtomBuffer[g][m];
-                        if (xp)
+                        if (!xp.empty())
                         {
                             pgrp->xp[m] += comm->pbcAtomBuffer[g][m];
                         }
@@ -784,7 +798,7 @@ void pull_calc_coms(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc
                     pgrp->localWeights[i] = csw * std::cos(twopi_box * x[ii][pull->cosdim])
                                             + snw * std::sin(twopi_box * x[ii][pull->cosdim]);
                 }
-                if (xp)
+                if (!xp.empty())
                 {
                     csw                    = comBuffer[2][0];
                     snw                    = comBuffer[2][1];
@@ -808,12 +822,12 @@ void pull_calc_coms(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc
 using BoolVec = gmx::BasicVector<bool>;
 
 /* Returns whether the pull group obeys the PBC restrictions */
-static bool pullGroupObeysPbcRestrictions(const pull_group_work_t& group,
-                                          const BoolVec&           dimUsed,
-                                          const rvec*              x,
-                                          const t_pbc&             pbc,
-                                          const gmx::RVec&         x_pbc,
-                                          const real               pbcMargin)
+static bool pullGroupObeysPbcRestrictions(const pull_group_work_t&       group,
+                                          const BoolVec&                 dimUsed,
+                                          gmx::ArrayRef<const gmx::RVec> x,
+                                          const t_pbc&                   pbc,
+                                          const gmx::RVec&               x_pbc,
+                                          const real                     pbcMargin)
 {
     /* Determine which dimensions are relevant for PBC */
     BoolVec dimUsesPbc       = { false, false, false };
@@ -896,7 +910,7 @@ static bool pullGroupObeysPbcRestrictions(const pull_group_work_t& group,
     return true;
 }
 
-int pullCheckPbcWithinGroups(const pull_t& pull, const rvec* x, const t_pbc& pbc, real pbcMargin)
+int pullCheckPbcWithinGroups(const pull_t& pull, gmx::ArrayRef<const gmx::RVec> x, const t_pbc& pbc, real pbcMargin)
 {
     if (pbc.pbcType == PbcType::No)
     {
@@ -976,7 +990,7 @@ bool pullCheckPbcWithinGroup(const pull_t&                  pull,
     }
 
     return (pullGroupObeysPbcRestrictions(
-            group, dimUsed, as_rvec_array(x.data()), pbc, pull.comm.pbcAtomBuffer[groupNr], pbcMargin));
+            group, dimUsed, x, pbc, pull.comm.pbcAtomBuffer[groupNr], pbcMargin));
 }
 
 void setPrevStepPullComFromState(struct pull_t* pull, const t_state* state)
@@ -1019,7 +1033,11 @@ void allocStatePrevStepPullCom(t_state* state, const pull_t* pull)
     }
 }
 
-void initPullComFromPrevStep(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc* pbc, const rvec x[])
+void initPullComFromPrevStep(const t_commrec*               cr,
+                             pull_t*                        pull,
+                             gmx::ArrayRef<const real>      masses,
+                             t_pbc*                         pbc,
+                             gmx::ArrayRef<const gmx::RVec> x)
 {
     pull_comm_t* comm   = &pull->comm;
     size_t       ngroup = pull->group.size();
@@ -1048,7 +1066,7 @@ void initPullComFromPrevStep(const t_commrec* cr, pull_t* pull, const real* mass
                        "Groups with no atoms, or only one atom, should not "
                        "use the COM from the previous step as reference.");
 
-            rvec x_pbc = { 0, 0, 0 };
+            gmx::RVec x_pbc = { 0, 0, 0 };
             copy_rvec(comm->pbcAtomBuffer[g], x_pbc);
 
             if (debug)
@@ -1068,7 +1086,7 @@ void initPullComFromPrevStep(const t_commrec* cr, pull_t* pull, const real* mass
 
             if (pgrp->atomSet.numAtomsLocal() <= c_pullMaxNumLocalAtomsSingleThreaded)
             {
-                sum_com_part(pgrp, 0, pgrp->atomSet.numAtomsLocal(), x, nullptr, masses, pbc, x_pbc, &comSumsTotal);
+                sum_com_part(pgrp, 0, pgrp->atomSet.numAtomsLocal(), x, {}, masses, pbc, x_pbc, &comSumsTotal);
             }
             else
             {
@@ -1077,8 +1095,7 @@ void initPullComFromPrevStep(const t_commrec* cr, pull_t* pull, const real* mass
                 {
                     int ind_start = (pgrp->atomSet.numAtomsLocal() * (t + 0)) / pull->nthreads;
                     int ind_end   = (pgrp->atomSet.numAtomsLocal() * (t + 1)) / pull->nthreads;
-                    sum_com_part(
-                            pgrp, ind_start, ind_end, x, nullptr, masses, pbc, x_pbc, &pull->comSums[t]);
+                    sum_com_part(pgrp, ind_start, ind_end, x, {}, masses, pbc, x_pbc, &pull->comSums[t]);
                 }
 
                 /* Reduce the thread contributions to sum_com[0] */