Use ArrayRef(WithPadding) in constraint code
authorBerk Hess <hess@kth.se>
Fri, 17 Jan 2020 09:59:18 +0000 (10:59 +0100)
committerPaul Bauer <paul.bauer.q@gmail.com>
Mon, 20 Jan 2020 15:33:53 +0000 (16:33 +0100)
This change is only refactoring, except for adding padding to two
buffers used for flexible constraints. This avoid potential illegal
memory access one element beyond the buffer with certain SIMD
implementations. Note the this memory was not used nor modified.

Change-Id: I385b3007a27888d15741e737ccfcf3e3a4369d1e

14 files changed:
src/gromacs/domdec/domdec.h
src/gromacs/domdec/domdec_constraints.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/tests/constrtestrunners.cpp
src/gromacs/mdlib/tests/settletestrunners.cpp
src/gromacs/mdlib/update.cpp
src/gromacs/mdrun/minimize.cpp
src/gromacs/mdrun/shellfc.cpp
src/gromacs/modularsimulator/constraintelement.cpp

index b7d4fb56ddd45aa099c7e7f8370ed71f13f90a1b..439e86b5e7fb0fb608c2061894eb1866dcf50033 100644 (file)
@@ -237,7 +237,11 @@ void dd_move_f_vsites(struct gmx_domdec_t* dd, rvec* f, rvec* fshift);
 void dd_clear_f_vsites(struct gmx_domdec_t* dd, rvec* f);
 
 /*! \brief Move x0 and also x1 if x1!=NULL. bX1IsCoord tells if to do PBC on x1 */
-void dd_move_x_constraints(struct gmx_domdec_t* dd, const matrix box, rvec* x0, rvec* x1, gmx_bool bX1IsCoord);
+void dd_move_x_constraints(struct gmx_domdec_t*     dd,
+                           const matrix             box,
+                           gmx::ArrayRef<gmx::RVec> x0,
+                           gmx::ArrayRef<gmx::RVec> x1,
+                           gmx_bool                 bX1IsCoord);
 
 /*! \brief Communicates the coordinates involved in virtual sites */
 void dd_move_x_vsites(struct gmx_domdec_t* dd, const matrix box, rvec* x);
index 5bbe432cc23055f943ac2ed4d973cbd7ee91498b..ec39ab14941270368e8a22974d17f7811ec6bdf3 100644 (file)
@@ -103,11 +103,16 @@ struct gmx_domdec_constraints_t
     //! @endcond
 };
 
-void dd_move_x_constraints(gmx_domdec_t* dd, const matrix box, rvec* x0, rvec* x1, gmx_bool bX1IsCoord)
+void dd_move_x_constraints(gmx_domdec_t*            dd,
+                           const matrix             box,
+                           gmx::ArrayRef<gmx::RVec> x0,
+                           gmx::ArrayRef<gmx::RVec> x1,
+                           gmx_bool                 bX1IsCoord)
 {
     if (dd->constraint_comm)
     {
-        dd_move_x_specat(dd, dd->constraint_comm, box, x0, x1, bX1IsCoord);
+        dd_move_x_specat(dd, dd->constraint_comm, box, as_rvec_array(x0.data()),
+                         as_rvec_array(x1.data()), bX1IsCoord);
 
         ddReopenBalanceRegionCpu(dd);
     }
index e17685d74c2ed56f3f1fa101d5a8909baf61aaa8..14defa9ab07ab10cec0c7da249e770dc58518b77 100644 (file)
@@ -82,7 +82,6 @@
 #include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/listoflists.h"
 #include "gromacs/utility/pleasecite.h"
-#include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/txtdump.h"
 
 namespace gmx
@@ -114,20 +113,20 @@ public:
          int                   numSettles);
     ~Impl();
     void setConstraints(const gmx_localtop_t& top, const t_mdatoms& md);
-    bool apply(bool               bLog,
-               bool               bEner,
-               int64_t            step,
-               int                delta_step,
-               real               step_scaling,
-               rvec*              x,
-               rvec*              xprime,
-               rvec*              min_proj,
-               const matrix       box,
-               real               lambda,
-               real*              dvdlambda,
-               rvec*              v,
-               tensor*            vir,
-               ConstraintVariable econq);
+    bool apply(bool                      bLog,
+               bool                      bEner,
+               int64_t                   step,
+               int                       delta_step,
+               real                      step_scaling,
+               ArrayRefWithPadding<RVec> x,
+               ArrayRefWithPadding<RVec> xprime,
+               ArrayRef<RVec>            min_proj,
+               const matrix              box,
+               real                      lambda,
+               real*                     dvdlambda,
+               ArrayRefWithPadding<RVec> v,
+               tensor*                   vir,
+               ConstraintVariable        econq);
     //! The total number of constraints.
     int ncon_tot = 0;
     //! The number of flexible constraints.
@@ -207,7 +206,7 @@ bool Constraints::havePerturbedConstraints() const
 }
 
 //! Clears constraint quantities for atoms in nonlocal region.
-static void clear_constraint_quantity_nonlocal(gmx_domdec_t* dd, rvec* q)
+static void clear_constraint_quantity_nonlocal(gmx_domdec_t* dd, ArrayRef<RVec> q)
 {
     int nonlocal_at_start, nonlocal_at_end, at;
 
@@ -232,14 +231,14 @@ void too_many_constraint_warnings(int eConstrAlg, int warncount)
 }
 
 //! Writes out coordinates.
-static void write_constr_pdb(const char*       fn,
-                             const char*       title,
-                             const gmx_mtop_t& mtop,
-                             int               start,
-                             int               homenr,
-                             const t_commrec*  cr,
-                             const rvec        x[],
-                             const matrix      box)
+static void write_constr_pdb(const char*          fn,
+                             const char*          title,
+                             const gmx_mtop_t&    mtop,
+                             int                  start,
+                             int                  homenr,
+                             const t_commrec*     cr,
+                             ArrayRef<const RVec> x,
+                             const matrix         box)
 {
     char          fname[STRLEN];
     FILE*         out;
@@ -294,15 +293,15 @@ static void write_constr_pdb(const char*       fn,
 }
 
 //! Writes out domain contents to help diagnose crashes.
-static void dump_confs(FILE*             log,
-                       int64_t           step,
-                       const gmx_mtop_t& mtop,
-                       int               start,
-                       int               homenr,
-                       const t_commrec*  cr,
-                       const rvec        x[],
-                       rvec              xprime[],
-                       const matrix      box)
+static void dump_confs(FILE*                log,
+                       int64_t              step,
+                       const gmx_mtop_t&    mtop,
+                       int                  start,
+                       int                  homenr,
+                       const t_commrec*     cr,
+                       ArrayRef<const RVec> x,
+                       ArrayRef<const RVec> xprime,
+                       const matrix         box)
 {
     char buf[STRLEN], buf2[22];
 
@@ -323,39 +322,39 @@ static void dump_confs(FILE*             log,
     fprintf(stderr, "Wrote pdb files with previous and current coordinates\n");
 }
 
-bool Constraints::apply(bool               bLog,
-                        bool               bEner,
-                        int64_t            step,
-                        int                delta_step,
-                        real               step_scaling,
-                        rvec*              x,
-                        rvec*              xprime,
-                        rvec*              min_proj,
-                        const matrix       box,
-                        real               lambda,
-                        real*              dvdlambda,
-                        rvec*              v,
-                        tensor*            vir,
-                        ConstraintVariable econq)
+bool Constraints::apply(bool                      bLog,
+                        bool                      bEner,
+                        int64_t                   step,
+                        int                       delta_step,
+                        real                      step_scaling,
+                        ArrayRefWithPadding<RVec> x,
+                        ArrayRefWithPadding<RVec> xprime,
+                        ArrayRef<RVec>            min_proj,
+                        const matrix              box,
+                        real                      lambda,
+                        real*                     dvdlambda,
+                        ArrayRefWithPadding<RVec> v,
+                        tensor*                   vir,
+                        ConstraintVariable        econq)
 {
-    return impl_->apply(bLog, bEner, step, delta_step, step_scaling, x, xprime, min_proj, box,
-                        lambda, dvdlambda, v, vir, econq);
+    return impl_->apply(bLog, bEner, step, delta_step, step_scaling, std::move(x), std::move(xprime),
+                        min_proj, box, lambda, dvdlambda, std::move(v), vir, econq);
 }
 
-bool Constraints::Impl::apply(bool               bLog,
-                              bool               bEner,
-                              int64_t            step,
-                              int                delta_step,
-                              real               step_scaling,
-                              rvec*              x,
-                              rvec*              xprime,
-                              rvec*              min_proj,
-                              const matrix       box,
-                              real               lambda,
-                              real*              dvdlambda,
-                              rvec*              v,
-                              tensor*            vir,
-                              ConstraintVariable econq)
+bool Constraints::Impl::apply(bool                      bLog,
+                              bool                      bEner,
+                              int64_t                   step,
+                              int                       delta_step,
+                              real                      step_scaling,
+                              ArrayRefWithPadding<RVec> x,
+                              ArrayRefWithPadding<RVec> xprime,
+                              ArrayRef<RVec>            min_proj,
+                              const matrix              box,
+                              real                      lambda,
+                              real*                     dvdlambda,
+                              ArrayRefWithPadding<RVec> v,
+                              tensor*                   vir,
+                              ConstraintVariable        econq)
 {
     bool   bOK, bDump;
     int    start, homenr;
@@ -444,23 +443,24 @@ bool Constraints::Impl::apply(bool               bLog,
      */
     if (cr->dd)
     {
-        dd_move_x_constraints(cr->dd, box, x, xprime, econq == ConstraintVariable::Positions);
+        dd_move_x_constraints(cr->dd, box, x.unpaddedArrayRef(), xprime.unpaddedArrayRef(),
+                              econq == ConstraintVariable::Positions);
 
-        if (v != nullptr)
+        if (!v.empty())
         {
             /* We need to initialize the non-local components of v.
              * We never actually use these values, but we do increment them,
              * so we should avoid uninitialized variables and overflows.
              */
-            clear_constraint_quantity_nonlocal(cr->dd, v);
+            clear_constraint_quantity_nonlocal(cr->dd, v.unpaddedArrayRef());
         }
     }
 
     if (lincsd != nullptr)
     {
         bOK = constrain_lincs(bLog || bEner, ir, step, lincsd, md, cr, ms, x, xprime, min_proj, box,
-                              pbc_null, lambda, dvdlambda, invdt, v, vir != nullptr, vir_r_m_dr,
-                              econq, nrnb, maxwarn, &warncount_lincs);
+                              pbc_null, lambda, dvdlambda, invdt, v.unpaddedArrayRef(),
+                              vir != nullptr, vir_r_m_dr, econq, nrnb, maxwarn, &warncount_lincs);
         if (!bOK && maxwarn < INT_MAX)
         {
             if (log != nullptr)
@@ -474,8 +474,11 @@ bool Constraints::Impl::apply(bool               bLog,
 
     if (shaked != nullptr)
     {
-        bOK = constrain_shake(log, shaked, md.invmass, *idef, ir, x, xprime, min_proj, nrnb, lambda,
-                              dvdlambda, invdt, v, vir != nullptr, vir_r_m_dr, maxwarn < INT_MAX, econq);
+        bOK = constrain_shake(
+                log, shaked, md.invmass, *idef, ir, as_rvec_array(x.unpaddedArrayRef().data()),
+                as_rvec_array(xprime.unpaddedArrayRef().data()), as_rvec_array(min_proj.data()),
+                nrnb, lambda, dvdlambda, invdt, as_rvec_array(v.unpaddedArrayRef().data()),
+                vir != nullptr, vir_r_m_dr, maxwarn < INT_MAX, econq);
 
         if (!bOK && maxwarn < INT_MAX)
         {
@@ -505,14 +508,14 @@ bool Constraints::Impl::apply(bool               bLog,
                             clear_mat(vir_r_m_dr_th[th]);
                         }
 
-                        csettle(settled, nth, th, pbc_null, x[0], xprime[0], invdt, v ? v[0] : nullptr,
-                                vir != nullptr, th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th],
+                        csettle(settled, nth, th, pbc_null, x, xprime, invdt, v, vir != nullptr,
+                                th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th],
                                 th == 0 ? &bSettleErrorHasOccurred0 : &bSettleErrorHasOccurred[th]);
                     }
                     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
                 }
                 inc_nrnb(nrnb, eNR_SETTLE, nsettle);
-                if (v != nullptr)
+                if (!v.empty())
                 {
                     inc_nrnb(nrnb, eNR_CONSTR_V, nsettle * 3);
                 }
@@ -553,8 +556,8 @@ bool Constraints::Impl::apply(bool               bLog,
                         {
                             settle_proj(settled, econq, end_th - start_th,
                                         settle->iatoms + start_th * (1 + NRAL(F_SETTLE)), pbc_null,
-                                        x, xprime, min_proj, calcvir_atom_end,
-                                        th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th]);
+                                        x.unpaddedArrayRef(), xprime.unpaddedArrayRef(), min_proj,
+                                        calcvir_atom_end, th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th]);
                         }
                     }
                     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
@@ -644,7 +647,8 @@ bool Constraints::Impl::apply(bool               bLog,
 
     if (bDump)
     {
-        dump_confs(log, step, mtop, start, homenr, cr, x, xprime, box);
+        dump_confs(log, step, mtop, start, homenr, cr, x.unpaddedArrayRef(),
+                   xprime.unpaddedArrayRef(), box);
     }
 
     if (econq == ConstraintVariable::Positions)
@@ -660,19 +664,24 @@ bool Constraints::Impl::apply(bool               bLog,
                 t = ir.init_t;
             }
             set_pbc(&pbc, ir.pbcType, box);
-            pull_constraint(pull_work, &md, &pbc, cr, ir.delta_t, t, x, xprime, v, *vir);
+            pull_constraint(pull_work, &md, &pbc, cr, ir.delta_t, t,
+                            as_rvec_array(x.unpaddedArrayRef().data()),
+                            as_rvec_array(xprime.unpaddedArrayRef().data()),
+                            as_rvec_array(v.unpaddedArrayRef().data()), *vir);
         }
         if (ed && delta_step > 0)
         {
             /* apply the essential dynamics constraints here */
-            do_edsam(&ir, step, cr, xprime, v, box, ed);
+            do_edsam(&ir, step, cr, as_rvec_array(xprime.unpaddedArrayRef().data()),
+                     as_rvec_array(v.unpaddedArrayRef().data()), box, ed);
         }
     }
     wallcycle_stop(wcycle, ewcCONSTR);
 
-    if (v != nullptr && md.cFREEZE)
+    if (!v.empty() && md.cFREEZE)
     {
         /* Set the velocities of frozen dimensions to zero */
+        ArrayRef<RVec> vRef = v.unpaddedArrayRef();
 
         int gmx_unused numThreads = gmx_omp_nthreads_get(emntUpdate);
 
@@ -685,7 +694,7 @@ bool Constraints::Impl::apply(bool               bLog,
             {
                 if (ir.opts.nFreeze[freezeGroup][d])
                 {
-                    v[i][d] = 0;
+                    vRef[i][d] = 0;
                 }
             }
         }
@@ -1136,15 +1145,8 @@ void do_constrain_first(FILE*                     fplog,
     int64_t step;
     real    dt = ir->delta_t;
     real    dvdl_dum;
-    rvec*   savex;
 
-    auto xRvec = as_rvec_array(x.paddedArrayRef().data());
-    auto vRvec = as_rvec_array(v.paddedArrayRef().data());
-
-    /* We need to allocate one element extra, since we might use
-     * (unaligned) 4-wide SIMD loads to access rvec entries.
-     */
-    snew(savex, natoms + 1);
+    PaddedVector<RVec> savex(natoms);
 
     start = 0;
     end   = md->homenr;
@@ -1163,14 +1165,14 @@ void do_constrain_first(FILE*                     fplog,
     dvdl_dum = 0;
 
     /* constrain the current position */
-    constr->apply(TRUE, FALSE, step, 0, 1.0, xRvec, xRvec, nullptr, box, lambda, &dvdl_dum, nullptr,
-                  nullptr, gmx::ConstraintVariable::Positions);
+    constr->apply(TRUE, FALSE, step, 0, 1.0, x, x, ArrayRef<RVec>(), box, lambda, &dvdl_dum,
+                  ArrayRefWithPadding<RVec>(), nullptr, gmx::ConstraintVariable::Positions);
     if (EI_VV(ir->eI))
     {
         /* constrain the inital velocity, and save it */
         /* also may be useful if we need the ekin from the halfstep for velocity verlet */
-        constr->apply(TRUE, FALSE, step, 0, 1.0, xRvec, vRvec, vRvec, box, lambda, &dvdl_dum,
-                      nullptr, nullptr, gmx::ConstraintVariable::Velocities);
+        constr->apply(TRUE, FALSE, step, 0, 1.0, x, v, v.unpaddedArrayRef(), box, lambda, &dvdl_dum,
+                      ArrayRefWithPadding<RVec>(), nullptr, gmx::ConstraintVariable::Velocities);
     }
     /* constrain the inital velocities at t-dt/2 */
     if (EI_STATE_VELOCITY(ir->eI) && ir->eI != eiVV)
@@ -1196,8 +1198,8 @@ void do_constrain_first(FILE*                     fplog,
             fprintf(fplog, "\nConstraining the coordinates at t0-dt (step %s)\n", gmx_step_str(step, buf));
         }
         dvdl_dum = 0;
-        constr->apply(TRUE, FALSE, step, -1, 1.0, xRvec, savex, nullptr, box, lambda, &dvdl_dum,
-                      vRvec, nullptr, gmx::ConstraintVariable::Positions);
+        constr->apply(TRUE, FALSE, step, -1, 1.0, x, savex.arrayRefWithPadding(), ArrayRef<RVec>(),
+                      box, lambda, &dvdl_dum, v, nullptr, gmx::ConstraintVariable::Positions);
 
         for (i = start; i < end; i++)
         {
@@ -1208,7 +1210,6 @@ void do_constrain_first(FILE*                     fplog,
             }
         }
     }
-    sfree(savex);
 }
 
 } // namespace gmx
index 9824e43a977d6c78f1a60295b0010e3fc552d4f1..c06d7d3576f3cf6dde8cf3130766aa9936ae6e2f 100644 (file)
@@ -49,6 +49,7 @@
 
 #include <cstdio>
 
+#include "gromacs/math/arrayrefwithpadding.h"
 #include "gromacs/math/vectypes.h"
 #include "gromacs/topology/idef.h"
 #include "gromacs/utility/arrayref.h"
@@ -170,20 +171,20 @@ public:
      *
      * /note x is non-const, because non-local atoms need to be communicated.
      */
-    bool apply(bool               bLog,
-               bool               bEner,
-               int64_t            step,
-               int                delta_step,
-               real               step_scaling,
-               rvec*              x,
-               rvec*              xprime,
-               rvec*              min_proj,
-               const matrix       box,
-               real               lambda,
-               real*              dvdlambda,
-               rvec*              v,
-               tensor*            vir,
-               ConstraintVariable econq);
+    bool apply(bool                      bLog,
+               bool                      bEner,
+               int64_t                   step,
+               int                       delta_step,
+               real                      step_scaling,
+               ArrayRefWithPadding<RVec> x,
+               ArrayRefWithPadding<RVec> xprime,
+               ArrayRef<RVec>            min_proj,
+               const matrix              box,
+               real                      lambda,
+               real*                     dvdlambda,
+               ArrayRefWithPadding<RVec> v,
+               tensor*                   vir,
+               ConstraintVariable        econq);
     //! Links the essentialdynamics and constraint code.
     void saveEdsamPointer(gmx_edsam* ed);
     //! Getter for use by domain decomposition.
index 98fdc2e03c04bea63cfbd928990de56d222c079f..b425d12eafdd4429b045e3144ffcc2bc1b963038 100644 (file)
@@ -569,17 +569,17 @@ static void gmx_simdcall calc_dr_x_f_simd(int                           b0,
 #endif // GMX_SIMD_HAVE_REAL
 
 /*! \brief LINCS projection, works on derivatives of the coordinates. */
-static void do_lincsp(const rvec*        x,
-                      rvec*              f,
-                      rvec*              fp,
-                      t_pbc*             pbc,
-                      Lincs*             lincsd,
-                      int                th,
-                      real*              invmass,
-                      ConstraintVariable econq,
-                      bool               bCalcDHDL,
-                      bool               bCalcVir,
-                      tensor             rmdf)
+static void do_lincsp(ArrayRefWithPadding<const RVec> xPadded,
+                      ArrayRefWithPadding<RVec>       fPadded,
+                      ArrayRef<RVec>                  fp,
+                      t_pbc*                          pbc,
+                      Lincs*                          lincsd,
+                      int                             th,
+                      real*                           invmass,
+                      ConstraintVariable              econq,
+                      bool                            bCalcDHDL,
+                      bool                            bCalcVir,
+                      tensor                          rmdf)
 {
     const int b0 = lincsd->task[th].b0;
     const int b1 = lincsd->task[th].b1;
@@ -608,6 +608,9 @@ static void do_lincsp(const rvec*        x,
     gmx::ArrayRef<real> rhs2 = lincsd->tmp2;
     gmx::ArrayRef<real> sol  = lincsd->tmp3;
 
+    const rvec* x = as_rvec_array(xPadded.paddedArrayRef().data());
+    rvec*       f = as_rvec_array(fPadded.paddedArrayRef().data());
+
 #if GMX_SIMD_HAVE_REAL
     /* This SIMD code does the same as the plain-C code after the #else.
      * The only difference is that we always call pbc code, as with SIMD
@@ -706,8 +709,8 @@ static void do_lincsp(const rvec*        x,
     /* When constraining forces, we should not use mass weighting,
      * so we pass invmass=NULL, which results in the use of 1 for all atoms.
      */
-    lincs_update_atoms(lincsd, th, 1.0, sol, r,
-                       (econq != ConstraintVariable::Force) ? invmass : nullptr, fp);
+    lincs_update_atoms(lincsd, th, 1.0, sol, r, (econq != ConstraintVariable::Force) ? invmass : nullptr,
+                       as_rvec_array(fp.data()));
 
     if (bCalcDHDL)
     {
@@ -945,22 +948,26 @@ static void gmx_simdcall calc_dist_iter_simd(int                           b0,
 #endif // GMX_SIMD_HAVE_REAL
 
 //! Implements LINCS constraining.
-static void do_lincs(const rvec*      x,
-                     rvec*            xp,
-                     const matrix     box,
-                     t_pbc*           pbc,
-                     Lincs*           lincsd,
-                     int              th,
-                     const real*      invmass,
-                     const t_commrec* cr,
-                     bool             bCalcDHDL,
-                     real             wangle,
-                     bool*            bWarn,
-                     real             invdt,
-                     rvec* gmx_restrict v,
-                     bool               bCalcVir,
-                     tensor             vir_r_m_dr)
+static void do_lincs(ArrayRefWithPadding<const RVec> xPadded,
+                     ArrayRefWithPadding<RVec>       xpPadded,
+                     const matrix                    box,
+                     t_pbc*                          pbc,
+                     Lincs*                          lincsd,
+                     int                             th,
+                     const real*                     invmass,
+                     const t_commrec*                cr,
+                     bool                            bCalcDHDL,
+                     real                            wangle,
+                     bool*                           bWarn,
+                     real                            invdt,
+                     ArrayRef<RVec>                  vRef,
+                     bool                            bCalcVir,
+                     tensor                          vir_r_m_dr)
 {
+    const rvec* x        = as_rvec_array(xPadded.paddedArrayRef().data());
+    rvec*       xp       = as_rvec_array(xpPadded.paddedArrayRef().data());
+    rvec* gmx_restrict v = as_rvec_array(vRef.data());
+
     const int b0 = lincsd->task[th].b0;
     const int b1 = lincsd->task[th].b1;
 
@@ -1098,7 +1105,7 @@ static void do_lincs(const rvec*      x,
                 /* Communicate the corrected non-local coordinates */
                 if (DOMAINDECOMP(cr))
                 {
-                    dd_move_x_constraints(cr->dd, box, xp, nullptr, FALSE);
+                    dd_move_x_constraints(cr->dd, box, xpPadded.unpaddedArrayRef(), ArrayRef<RVec>(), FALSE);
                 }
             }
 #pragma omp barrier
@@ -2124,8 +2131,8 @@ void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_
 
 //! Issues a warning when LINCS constraints cannot be satisfied.
 static void lincs_warning(gmx_domdec_t*                 dd,
-                          const rvec*                   x,
-                          rvec*                         xprime,
+                          ArrayRef<const RVec>          x,
+                          ArrayRef<const RVec>          xprime,
                           t_pbc*                        pbc,
                           int                           ncons,
                           gmx::ArrayRef<const AtomPair> atoms,
@@ -2192,7 +2199,7 @@ struct LincsDeviations
 };
 
 //! Determine how well the constraints have been satisfied.
-static LincsDeviations makeLincsDeviations(const Lincs& lincsd, const rvec* x, const t_pbc* pbc)
+static LincsDeviations makeLincsDeviations(const Lincs& lincsd, ArrayRef<const RVec> x, const t_pbc* pbc)
 {
     LincsDeviations                result;
     const ArrayRef<const AtomPair> atoms  = lincsd.atoms;
@@ -2241,28 +2248,28 @@ static LincsDeviations makeLincsDeviations(const Lincs& lincsd, const rvec* x, c
     return result;
 }
 
-bool constrain_lincs(bool                  computeRmsd,
-                     const t_inputrec&     ir,
-                     int64_t               step,
-                     Lincs*                lincsd,
-                     const t_mdatoms&      md,
-                     const t_commrec*      cr,
-                     const gmx_multisim_t* ms,
-                     const rvec*           x,
-                     rvec*                 xprime,
-                     rvec*                 min_proj,
-                     const matrix          box,
-                     t_pbc*                pbc,
-                     real                  lambda,
-                     real*                 dvdlambda,
-                     real                  invdt,
-                     rvec*                 v,
-                     bool                  bCalcVir,
-                     tensor                vir_r_m_dr,
-                     ConstraintVariable    econq,
-                     t_nrnb*               nrnb,
-                     int                   maxwarn,
-                     int*                  warncount)
+bool constrain_lincs(bool                            computeRmsd,
+                     const t_inputrec&               ir,
+                     int64_t                         step,
+                     Lincs*                          lincsd,
+                     const t_mdatoms&                md,
+                     const t_commrec*                cr,
+                     const gmx_multisim_t*           ms,
+                     ArrayRefWithPadding<const RVec> xPadded,
+                     ArrayRefWithPadding<RVec>       xprimePadded,
+                     ArrayRef<RVec>                  min_proj,
+                     const matrix                    box,
+                     t_pbc*                          pbc,
+                     real                            lambda,
+                     real*                           dvdlambda,
+                     real                            invdt,
+                     ArrayRef<RVec>                  v,
+                     bool                            bCalcVir,
+                     tensor                          vir_r_m_dr,
+                     ConstraintVariable              econq,
+                     t_nrnb*                         nrnb,
+                     int                             maxwarn,
+                     int*                            warncount)
 {
     bool bOK = TRUE;
 
@@ -2282,6 +2289,9 @@ bool constrain_lincs(bool                  computeRmsd,
         return bOK;
     }
 
+    ArrayRef<const RVec> x      = xPadded.unpaddedArrayRef();
+    ArrayRef<RVec>       xprime = xprimePadded.unpaddedArrayRef();
+
     if (econq == ConstraintVariable::Positions)
     {
         /* We can't use bCalcDHDL here, since NULL can be passed for dvdlambda
@@ -2355,8 +2365,9 @@ bool constrain_lincs(bool                  computeRmsd,
 
                 clear_mat(lincsd->task[th].vir_r_m_dr);
 
-                do_lincs(x, xprime, box, pbc, lincsd, th, md.invmass, cr, bCalcDHDL, ir.LincsWarnAngle,
-                         &bWarn, invdt, v, bCalcVir, th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr);
+                do_lincs(xPadded, xprimePadded, box, pbc, lincsd, th, md.invmass, cr, bCalcDHDL,
+                         ir.LincsWarnAngle, &bWarn, invdt, v, bCalcVir,
+                         th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr);
             }
             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
         }
@@ -2433,8 +2444,8 @@ bool constrain_lincs(bool                  computeRmsd,
             {
                 int th = gmx_omp_get_thread_num();
 
-                do_lincsp(x, xprime, min_proj, pbc, lincsd, th, md.invmass, econq, bCalcDHDL,
-                          bCalcVir, th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr);
+                do_lincsp(xPadded, xprimePadded, min_proj, pbc, lincsd, th, md.invmass, econq,
+                          bCalcDHDL, bCalcVir, th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr);
             }
             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
         }
@@ -2475,7 +2486,7 @@ bool constrain_lincs(bool                  computeRmsd,
     {
         inc_nrnb(nrnb, eNR_LINCSMAT, lincsd->nOrder * lincsd->ncc_triangle);
     }
-    if (v)
+    if (!v.empty())
     {
         inc_nrnb(nrnb, eNR_CONSTR_V, lincsd->nc_real * 2);
     }
index ef43c3a5865bad17c99e917c760f80a1468c065e..649c58d661bd7b2af9214d9f2de157fe52a27c70 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020, 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.
@@ -46,6 +46,7 @@
 
 #include <cstdio>
 
+#include "gromacs/math/arrayrefwithpadding.h"
 #include "gromacs/math/vectypes.h"
 #include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/basedefinitions.h"
@@ -93,28 +94,28 @@ void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_
 /*! \brief Applies LINCS constraints.
  *
  * \returns true if the constraining succeeded. */
-bool constrain_lincs(bool                  computeRmsd,
-                     const t_inputrec&     ir,
-                     int64_t               step,
-                     Lincs*                lincsd,
-                     const t_mdatoms&      md,
-                     const t_commrec*      cr,
-                     const gmx_multisim_t* ms,
-                     const rvec*           x,
-                     rvec*                 xprime,
-                     rvec*                 min_proj,
-                     const matrix          box,
-                     t_pbc*                pbc,
-                     real                  lambda,
-                     real*                 dvdlambda,
-                     real                  invdt,
-                     rvec*                 v,
-                     bool                  bCalcVir,
-                     tensor                vir_r_m_dr,
-                     ConstraintVariable    econq,
-                     t_nrnb*               nrnb,
-                     int                   maxwarn,
-                     int*                  warncount);
+bool constrain_lincs(bool                            computeRmsd,
+                     const t_inputrec&               ir,
+                     int64_t                         step,
+                     Lincs*                          lincsd,
+                     const t_mdatoms&                md,
+                     const t_commrec*                cr,
+                     const gmx_multisim_t*           ms,
+                     ArrayRefWithPadding<const RVec> x,
+                     ArrayRefWithPadding<RVec>       xprime,
+                     ArrayRef<RVec>                  min_proj,
+                     const matrix                    box,
+                     t_pbc*                          pbc,
+                     real                            lambda,
+                     real*                           dvdlambda,
+                     real                            invdt,
+                     ArrayRef<RVec>                  v,
+                     bool                            bCalcVir,
+                     tensor                          vir_r_m_dr,
+                     ConstraintVariable              econq,
+                     t_nrnb*                         nrnb,
+                     int                             maxwarn,
+                     int*                            warncount);
 
 } // namespace gmx
 
index f8f183bd46e8502c4400fa186e64e4824459e7bd..2fa116a051d38530f129c9f17a0d8e5103b2120e 100644 (file)
@@ -306,16 +306,16 @@ void settle_set_constraints(settledata* settled, const t_ilist* il_settle, const
     }
 }
 
-void settle_proj(settledata*        settled,
-                 ConstraintVariable econq,
-                 int                nsettle,
-                 const t_iatom      iatoms[],
-                 const t_pbc*       pbc,
-                 const rvec         x[],
-                 rvec*              der,
-                 rvec*              derp,
-                 int                calcvir_atom_end,
-                 tensor             vir_r_m_dder)
+void settle_proj(settledata*          settled,
+                 ConstraintVariable   econq,
+                 int                  nsettle,
+                 const t_iatom        iatoms[],
+                 const t_pbc*         pbc,
+                 ArrayRef<const RVec> x,
+                 ArrayRef<RVec>       der,
+                 ArrayRef<RVec>       derp,
+                 int                  calcvir_atom_end,
+                 tensor               vir_r_m_dder)
 {
     /* Settle for projection out constraint components
      * of derivatives of the coordinates.
@@ -812,18 +812,22 @@ static void settleTemplateWrapper(settledata* settled,
     }
 }
 
-void csettle(settledata*  settled,
-             int          nthread,
-             int          thread,
-             const t_pbc* pbc,
-             const real   x[],
-             real         xprime[],
-             real         invdt,
-             real*        v,
-             bool         bCalcVirial,
-             tensor       vir_r_m_dr,
-             bool*        bErrorHasOccurred)
+void csettle(settledata*                     settled,
+             int                             nthread,
+             int                             thread,
+             const t_pbc*                    pbc,
+             ArrayRefWithPadding<const RVec> x,
+             ArrayRefWithPadding<RVec>       xprime,
+             real                            invdt,
+             ArrayRefWithPadding<RVec>       v,
+             bool                            bCalcVirial,
+             tensor                          vir_r_m_dr,
+             bool*                           bErrorHasOccurred)
 {
+    const real* xPtr      = as_rvec_array(x.paddedArrayRef().data())[0];
+    real*       xprimePtr = as_rvec_array(xprime.paddedArrayRef().data())[0];
+    real*       vPtr      = as_rvec_array(v.paddedArrayRef().data())[0];
+
 #if GMX_SIMD_HAVE_REAL
     if (settled->bUseSimd)
     {
@@ -832,8 +836,8 @@ void csettle(settledata*  settled,
         set_pbc_simd(pbc, pbcSimd);
 
         settleTemplateWrapper<SimdReal, SimdBool, GMX_SIMD_REAL_WIDTH, const real*>(
-                settled, nthread, thread, pbcSimd, x, xprime, invdt, v, bCalcVirial, vir_r_m_dr,
-                bErrorHasOccurred);
+                settled, nthread, thread, pbcSimd, xPtr, xprimePtr, invdt, vPtr, bCalcVirial,
+                vir_r_m_dr, bErrorHasOccurred);
     }
     else
 #endif
@@ -852,9 +856,9 @@ void csettle(settledata*  settled,
             pbcNonNull = &pbcNo;
         }
 
-        settleTemplateWrapper<real, bool, 1, const t_pbc*>(settled, nthread, thread, pbcNonNull, x,
-                                                           xprime, invdt, v, bCalcVirial,
-                                                           vir_r_m_dr, bErrorHasOccurred);
+        settleTemplateWrapper<real, bool, 1, const t_pbc*>(settled, nthread, thread, pbcNonNull,
+                                                           &xPtr[0], &xprimePtr[0], invdt, &vPtr[0],
+                                                           bCalcVirial, vir_r_m_dr, bErrorHasOccurred);
     }
 }
 
index f503ac8481edce84c34316dddd775e18b25afd9f..7cc0a89beab97ef9b0ae971375f4b68f21f021d8 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020, 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.
@@ -44,6 +44,7 @@
 #ifndef GMX_MDLIB_SETTLE_H
 #define GMX_MDLIB_SETTLE_H
 
+#include "gromacs/math/arrayrefwithpadding.h"
 #include "gromacs/topology/idef.h"
 
 struct gmx_cmap_t;
@@ -72,32 +73,32 @@ void settle_set_constraints(settledata* settled, const t_ilist* il_settle, const
 /*! \brief Constrain coordinates using SETTLE.
  * Can be called on any number of threads.
  */
-void csettle(settledata*  settled,          /* The SETTLE structure */
-             int          nthread,          /* The number of threads used */
-             int          thread,           /* Our thread index */
-             const t_pbc* pbc,              /* PBC data pointer, can be NULL */
-             const real   x[],              /* Reference coordinates */
-             real         xprime[],         /* New coords, to be settled */
-             real         invdt,            /* 1/delta_t */
-             real*        v,                /* Also constrain v if v!=NULL */
-             bool         bCalcVirial,      /* Calculate the virial contribution */
-             tensor       vir_r_m_dr,       /* sum r x m delta_r */
-             bool*        bErrorHasOccurred /* True if a settle error occurred */
+void csettle(settledata*                     settled,     /* The SETTLE structure */
+             int                             nthread,     /* The number of threads used */
+             int                             thread,      /* Our thread index */
+             const t_pbc*                    pbc,         /* PBC data pointer, can be NULL */
+             ArrayRefWithPadding<const RVec> x,           /* Reference coordinates */
+             ArrayRefWithPadding<RVec>       xprime,      /* New coords, to be settled */
+             real                            invdt,       /* 1/delta_t */
+             ArrayRefWithPadding<RVec>       v,           /* Also constrain v if v!=NULL */
+             bool                            bCalcVirial, /* Calculate the virial contribution */
+             tensor                          vir_r_m_dr,  /* sum r x m delta_r */
+             bool*                           bErrorHasOccurred /* True if a settle error occurred */
 );
 
 /*! \brief Analytical algorithm to subtract the components of derivatives
  * of coordinates working on settle type constraint.
  */
-void settle_proj(settledata*        settled,
-                 ConstraintVariable econq,
-                 int                nsettle,
-                 const t_iatom      iatoms[],
-                 const t_pbc*       pbc, /* PBC data pointer, can be NULL  */
-                 const rvec         x[],
-                 rvec*              der,
-                 rvec*              derp,
-                 int                CalcVirAtomEnd,
-                 tensor             vir_r_m_dder);
+void settle_proj(settledata*          settled,
+                 ConstraintVariable   econq,
+                 int                  nsettle,
+                 const t_iatom        iatoms[],
+                 const t_pbc*         pbc, /* PBC data pointer, can be NULL  */
+                 ArrayRef<const RVec> x,
+                 ArrayRef<RVec>       der,
+                 ArrayRef<RVec>       derp,
+                 int                  CalcVirAtomEnd,
+                 tensor               vir_r_m_dder);
 
 } // namespace gmx
 
index b1f33560927122ab76ff6aa26b1b749850dd05a6..5f974aaf8bd3ae77eafcdfc3ac1c803c31c68888 100644 (file)
@@ -131,11 +131,12 @@ void applyLincs(ConstraintsTestData* testData, t_pbc pbc)
     // Evaluate constraints
     bool success = constrain_lincs(
             false, testData->ir_, 0, lincsd, testData->md_, &testData->cr_, &testData->ms_,
-            as_rvec_array(testData->x_.data()), as_rvec_array(testData->xPrime_.data()),
-            as_rvec_array(testData->xPrime2_.data()), pbc.box, &pbc, testData->md_.lambda,
-            &testData->dHdLambda_, testData->invdt_, as_rvec_array(testData->v_.data()),
-            testData->computeVirial_, testData->virialScaled_, gmx::ConstraintVariable::Positions,
-            &testData->nrnb_, maxwarn, &warncount_lincs);
+            testData->x_.arrayRefWithPadding(), testData->xPrime_.arrayRefWithPadding(),
+            testData->xPrime2_.arrayRefWithPadding().unpaddedArrayRef(), pbc.box, &pbc,
+            testData->md_.lambda, &testData->dHdLambda_, testData->invdt_,
+            testData->v_.arrayRefWithPadding().unpaddedArrayRef(), testData->computeVirial_,
+            testData->virialScaled_, gmx::ConstraintVariable::Positions, &testData->nrnb_, maxwarn,
+            &warncount_lincs);
     EXPECT_TRUE(success) << "Test failed with a false return value in LINCS.";
     EXPECT_EQ(warncount_lincs, 0) << "There were warnings in LINCS.";
     done_lincs(lincsd);
index 88a4c1774e74d111114b5885c6976a8a643c518b..59cc352b87fac927f373f353ad7e8f2fc69654ba 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020, 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.
@@ -70,10 +70,9 @@ void applySettle(SettleTestData*    testData,
     bool errorOccured;
     int  numThreads  = 1;
     int  threadIndex = 0;
-    csettle(settled, numThreads, threadIndex, &pbc,
-            reinterpret_cast<real*>(as_rvec_array(testData->x_.data())),
-            reinterpret_cast<real*>(as_rvec_array(testData->xPrime_.data())), testData->reciprocalTimeStep_,
-            updateVelocities ? reinterpret_cast<real*>(as_rvec_array(testData->v_.data())) : nullptr,
+    csettle(settled, numThreads, threadIndex, &pbc, testData->x_.arrayRefWithPadding(),
+            testData->xPrime_.arrayRefWithPadding(), testData->reciprocalTimeStep_,
+            updateVelocities ? testData->v_.arrayRefWithPadding() : ArrayRefWithPadding<RVec>(),
             calcVirial, testData->virial_, &errorOccured);
     settle_free(settled);
     EXPECT_FALSE(errorOccured) << testDescription;
index 1afa13fecfd06e9afcc901541d46565cd90923ab..275bdf68c159f07567648631fefbd3a214c50ffa 100644 (file)
@@ -1465,9 +1465,10 @@ void constrain_velocities(int64_t step,
         clear_mat(vir_part);
 
         /* Constrain the coordinates upd->xp */
-        constr->apply(do_log, do_ene, step, 1, 1.0, state->x.rvec_array(), state->v.rvec_array(),
-                      state->v.rvec_array(), state->box, state->lambda[efptBONDED], dvdlambda,
-                      nullptr, bCalcVir ? &vir_con : nullptr, ConstraintVariable::Velocities);
+        constr->apply(do_log, do_ene, step, 1, 1.0, state->x.arrayRefWithPadding(),
+                      state->v.arrayRefWithPadding(), state->v.arrayRefWithPadding().unpaddedArrayRef(),
+                      state->box, state->lambda[efptBONDED], dvdlambda, ArrayRefWithPadding<RVec>(),
+                      bCalcVir ? &vir_con : nullptr, ConstraintVariable::Velocities);
 
         if (bCalcVir)
         {
@@ -1498,10 +1499,10 @@ void constrain_coordinates(int64_t step,
         clear_mat(vir_part);
 
         /* Constrain the coordinates upd->xp */
-        constr->apply(do_log, do_ene, step, 1, 1.0, state->x.rvec_array(), upd->xp()->rvec_array(),
-                      nullptr, state->box, state->lambda[efptBONDED], dvdlambda,
-                      as_rvec_array(state->v.data()), bCalcVir ? &vir_con : nullptr,
-                      ConstraintVariable::Positions);
+        constr->apply(do_log, do_ene, step, 1, 1.0, state->x.arrayRefWithPadding(),
+                      upd->xp()->arrayRefWithPadding(), ArrayRef<RVec>(), state->box,
+                      state->lambda[efptBONDED], dvdlambda, state->v.arrayRefWithPadding(),
+                      bCalcVir ? &vir_con : nullptr, ConstraintVariable::Positions);
 
         if (bCalcVir)
         {
@@ -1564,9 +1565,10 @@ void update_sd_second_half(int64_t step,
         wallcycle_stop(wcycle, ewcUPDATE);
 
         /* Constrain the coordinates upd->xp for half a time step */
-        constr->apply(do_log, do_ene, step, 1, 0.5, state->x.rvec_array(), upd->xp()->rvec_array(),
-                      nullptr, state->box, state->lambda[efptBONDED], dvdlambda,
-                      as_rvec_array(state->v.data()), nullptr, ConstraintVariable::Positions);
+        constr->apply(do_log, do_ene, step, 1, 0.5, state->x.arrayRefWithPadding(),
+                      upd->xp()->arrayRefWithPadding(), ArrayRef<RVec>(), state->box,
+                      state->lambda[efptBONDED], dvdlambda, state->v.arrayRefWithPadding(), nullptr,
+                      ConstraintVariable::Positions);
     }
 }
 
index 593d432b37000b8bb3830c2f9f879a904e764ced..64da82baacff911034f4cecd9cf3da12b1987ee5 100644 (file)
 #include "legacysimulator.h"
 #include "shellfc.h"
 
+using gmx::ArrayRef;
 using gmx::MdrunScheduleWorkload;
+using gmx::RVec;
 
 //! Utility structure for manipulating states during EM
 typedef struct
@@ -455,8 +457,9 @@ static void init_em(FILE*                fplog,
         {
             /* Constrain the starting coordinates */
             dvdl_constr = 0;
-            constr->apply(TRUE, TRUE, -1, 0, 1.0, ems->s.x.rvec_array(), ems->s.x.rvec_array(),
-                          nullptr, ems->s.box, ems->s.lambda[efptFEP], &dvdl_constr, nullptr,
+            constr->apply(TRUE, TRUE, -1, 0, 1.0, ems->s.x.arrayRefWithPadding(),
+                          ems->s.x.arrayRefWithPadding(), ArrayRef<RVec>(), ems->s.box,
+                          ems->s.lambda[efptFEP], &dvdl_constr, gmx::ArrayRefWithPadding<RVec>(),
                           nullptr, gmx::ConstraintVariable::Positions);
         }
     }
@@ -680,9 +683,10 @@ static bool do_em_step(const t_commrec*                   cr,
     if (constr)
     {
         dvdl_constr = 0;
-        validStep = constr->apply(TRUE, TRUE, count, 0, 1.0, s1->x.rvec_array(), s2->x.rvec_array(),
-                                  nullptr, s2->box, s2->lambda[efptBONDED], &dvdl_constr, nullptr,
-                                  nullptr, gmx::ConstraintVariable::Positions);
+        validStep   = constr->apply(
+                TRUE, TRUE, count, 0, 1.0, s1->x.arrayRefWithPadding(), s2->x.arrayRefWithPadding(),
+                ArrayRef<RVec>(), s2->box, s2->lambda[efptBONDED], &dvdl_constr,
+                gmx::ArrayRefWithPadding<RVec>(), nullptr, gmx::ConstraintVariable::Positions);
 
         if (cr->nnodes > 1)
         {
@@ -887,11 +891,11 @@ void EnergyEvaluator::run(em_state_t* ems, rvec mu_tot, tensor vir, tensor pres,
     if (constr)
     {
         /* Project out the constraint components of the force */
-        dvdl_constr  = 0;
-        rvec* f_rvec = ems->f.rvec_array();
-        constr->apply(FALSE, FALSE, count, 0, 1.0, ems->s.x.rvec_array(), f_rvec, f_rvec,
-                      ems->s.box, ems->s.lambda[efptBONDED], &dvdl_constr, nullptr, &shake_vir,
-                      gmx::ConstraintVariable::ForceDispl);
+        dvdl_constr = 0;
+        auto f      = ems->f.arrayRefWithPadding();
+        constr->apply(FALSE, FALSE, count, 0, 1.0, ems->s.x.arrayRefWithPadding(), f,
+                      f.unpaddedArrayRef(), ems->s.box, ems->s.lambda[efptBONDED], &dvdl_constr,
+                      gmx::ArrayRefWithPadding<RVec>(), &shake_vir, gmx::ConstraintVariable::ForceDispl);
         enerd->term[F_DVDL_CONSTR] += dvdl_constr;
         m_add(force_vir, shake_vir, vir);
     }
index a864ef0d481378c4eb0f27d3a65f4a3873042f8a..094674d983927f18c4acf7ab6b117a28d8498a77 100644 (file)
@@ -115,8 +115,8 @@ struct gmx_shellfc_t
     /* Flexible constraint working data */
     std::vector<RVec>       acc_dir;                /* Acceleration direction for flexcon        */
     gmx::PaddedVector<RVec> x_old;                  /* Old coordinates for flexcon               */
-    std::vector<RVec>       adir_xnold;             /* Work space for init_adir                  */
-    std::vector<RVec>       adir_xnew;              /* Work space for init_adir                  */
+    gmx::PaddedVector<RVec> adir_xnold;             /* Work space for init_adir                  */
+    gmx::PaddedVector<RVec> adir_xnew;              /* Work space for init_adir                  */
     std::int64_t            numForceEvaluations;    /* Total number of force evaluations         */
     int                     numConvergedIterations; /* Total number of iterations that converged */
 };
@@ -872,8 +872,8 @@ static void init_adir(gmx_shellfc_t*            shfc,
     {
         n = end;
     }
-    shfc->adir_xnold.resize(n);
-    shfc->adir_xnew.resize(n);
+    shfc->adir_xnold.resizeWithPadding(n);
+    shfc->adir_xnew.resizeWithPadding(n);
     rvec* xnold = as_rvec_array(shfc->adir_xnold.data());
     rvec* xnew  = as_rvec_array(shfc->adir_xnew.data());
     rvec* x_old = as_rvec_array(xOld.paddedArrayRef().data());
@@ -902,10 +902,12 @@ static void init_adir(gmx_shellfc_t*            shfc,
             }
         }
     }
-    constr->apply(FALSE, FALSE, step, 0, 1.0, x, xnold, nullptr, box, lambda[efptBONDED],
-                  &(dvdlambda[efptBONDED]), nullptr, nullptr, gmx::ConstraintVariable::Positions);
-    constr->apply(FALSE, FALSE, step, 0, 1.0, x, xnew, nullptr, box, lambda[efptBONDED],
-                  &(dvdlambda[efptBONDED]), nullptr, nullptr, gmx::ConstraintVariable::Positions);
+    constr->apply(FALSE, FALSE, step, 0, 1.0, xCurrent, shfc->adir_xnold.arrayRefWithPadding(),
+                  ArrayRef<RVec>(), box, lambda[efptBONDED], &(dvdlambda[efptBONDED]),
+                  ArrayRefWithPadding<RVec>(), nullptr, gmx::ConstraintVariable::Positions);
+    constr->apply(FALSE, FALSE, step, 0, 1.0, xCurrent, shfc->adir_xnew.arrayRefWithPadding(),
+                  ArrayRef<RVec>(), box, lambda[efptBONDED], &(dvdlambda[efptBONDED]),
+                  ArrayRefWithPadding<RVec>(), nullptr, gmx::ConstraintVariable::Positions);
 
     for (n = 0; n < end; n++)
     {
@@ -918,9 +920,9 @@ static void init_adir(gmx_shellfc_t*            shfc,
     }
 
     /* Project the acceleration on the old bond directions */
-    constr->apply(FALSE, FALSE, step, 0, 1.0, x_old, xnew, as_rvec_array(acc_dir.data()), box,
-                  lambda[efptBONDED], &(dvdlambda[efptBONDED]), nullptr, nullptr,
-                  gmx::ConstraintVariable::Deriv_FlexCon);
+    constr->apply(FALSE, FALSE, step, 0, 1.0, xOld, shfc->adir_xnew.arrayRefWithPadding(), acc_dir,
+                  box, lambda[efptBONDED], &(dvdlambda[efptBONDED]), ArrayRefWithPadding<RVec>(),
+                  nullptr, gmx::ConstraintVariable::Deriv_FlexCon);
 }
 
 void relax_shell_flexcon(FILE*                               fplog,
index 90f993f00523a22204dfe60e6bae0d9e9f2acfe2..4a0b7fe545b6c456b0cfa0a1881919bb7d039a7f 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2019, by the GROMACS development team, led by
+ * Copyright (c) 2019,2020, 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.
@@ -125,7 +125,10 @@ void ConstraintsElement<variable>::apply(Step step, bool calculateVirial, bool w
 {
     tensor vir_con;
 
-    rvec *x, *xprime, *min_proj, *v;
+    ArrayRefWithPadding<RVec> x;
+    ArrayRefWithPadding<RVec> xprime;
+    ArrayRef<RVec>            min_proj;
+    ArrayRefWithPadding<RVec> v;
 
     // disabled functionality
     real  lambda    = 0;
@@ -134,16 +137,14 @@ void ConstraintsElement<variable>::apply(Step step, bool calculateVirial, bool w
     switch (variable)
     {
         case ConstraintVariable::Positions:
-            x = as_rvec_array(statePropagatorData_->previousPositionsView().paddedArrayRef().data());
-            xprime   = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
-            min_proj = nullptr;
-            v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
+            x      = statePropagatorData_->previousPositionsView();
+            xprime = statePropagatorData_->positionsView();
+            v      = statePropagatorData_->velocitiesView();
             break;
         case ConstraintVariable::Velocities:
-            x      = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
-            xprime = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
-            min_proj = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
-            v        = nullptr;
+            x        = statePropagatorData_->positionsView();
+            xprime   = statePropagatorData_->velocitiesView();
+            min_proj = statePropagatorData_->velocitiesView().unpaddedArrayRef();
             break;
         default: gmx_fatal(FARGS, "Constraint algorithm not implemented for modular simulator.");
     }