Make do_constrain_first() independent of t_state

This changes the do_constrain_first() function to directly take the
required state data (natoms, x, v, box, lambda[efptBONDED]) as input
instead of a pointer to the full t_state object. This makes subsequent
changes to the t_state object easier.

This also adds some trivial const qualifier for the box matrix in the
constraining functions.

Change-Id: I28b58df45481549cd93076334b0778e23e228154

diff --git a/src/gromacs/domdec/domdec.h b/src/gromacs/domdec/domdec.h
index 948d6ec55660560d801570207d59ec7c99542851..8c74c4204dcbaa6bc7e9550c70136d88b756e748 100644 (file)
--- a/src/gromacs/domdec/domdec.h
+++ b/src/gromacs/domdec/domdec.h
@@ -244,7 +244,7 @@ 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, matrix box,
+void dd_move_x_constraints(struct gmx_domdec_t *dd, const matrix box,
                            rvec *x0, rvec *x1, gmx_bool bX1IsCoord);
 
 /*! \brief Communicates the coordinates involved in virtual sites */

diff --git a/src/gromacs/domdec/domdec_constraints.cpp b/src/gromacs/domdec/domdec_constraints.cpp
index 5e3d382ca113b5218fc7253aa13c7a6ce4ae212c..8005f7a8fd6bbad25acaf155956a6e14786fff9a 100644 (file)
--- a/src/gromacs/domdec/domdec_constraints.cpp
+++ b/src/gromacs/domdec/domdec_constraints.cpp
@@ -99,7 +99,7 @@ struct gmx_domdec_constraints_t
     //! @endcond
 };
 
-void dd_move_x_constraints(gmx_domdec_t *dd, matrix box,
+void dd_move_x_constraints(gmx_domdec_t *dd, const matrix box,
                            rvec *x0, rvec *x1, gmx_bool bX1IsCoord)
 {
     if (dd->constraint_comm)

diff --git a/src/gromacs/essentialdynamics/edsam.cpp b/src/gromacs/essentialdynamics/edsam.cpp
index f4b2b2d0790e451cf661f213521e0b5017912a2b..2b9fd374a66c90f8f15d095520b8d4a91bee1775 100644 (file)
--- a/src/gromacs/essentialdynamics/edsam.cpp
+++ b/src/gromacs/essentialdynamics/edsam.cpp
@@ -1946,7 +1946,7 @@ void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
 }
 
 
-static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
+static inline void ed_unshift_single_coord(const matrix box, const rvec x, const ivec is, rvec xu)
 {
     int tx, ty, tz;
 
@@ -3000,7 +3000,7 @@ void do_edsam(const t_inputrec *ir,
               const t_commrec  *cr,
               rvec              xs[],
               rvec              v[],
-              matrix            box,
+              const matrix      box,
               gmx_edsam        *ed)
 {
     int                i, edinr, iupdate = 500;

diff --git a/src/gromacs/essentialdynamics/edsam.h b/src/gromacs/essentialdynamics/edsam.h
index 04cf6cdcd7867edbd22df0e1b9f83f033209223f..7fd2d28dc9bac790b270438f633b334fa72c2ab0 100644 (file)
--- a/src/gromacs/essentialdynamics/edsam.h
+++ b/src/gromacs/essentialdynamics/edsam.h
@@ -102,7 +102,7 @@ class MDLogger;
  * \param ed                The essential dynamics data.
  */
 void do_edsam(const t_inputrec *ir, int64_t step,
-              const t_commrec *cr, rvec xs[], rvec v[], matrix box, gmx_edsam *ed);
+              const t_commrec *cr, rvec xs[], rvec v[], const matrix box, gmx_edsam *ed);
 
 
 /*! \brief Initializes the essential dynamics and flooding module.

diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp
index 2be5fd3c3a553e9dc360114002286ba2b5660e15..6181835e4ff64c71e5d81b8e13b935a0cf7c6bfe 100644 (file)
--- a/src/gromacs/mdlib/constr.cpp
+++ b/src/gromacs/mdlib/constr.cpp
@@ -122,7 +122,7 @@ class Constraints::Impl
                    rvec                 *x,
                    rvec                 *xprime,
                    rvec                 *min_proj,
-                   matrix                box,
+                   const matrix          box,
                    real                  lambda,
                    real                 *dvdlambda,
                    rvec                 *v,
@@ -238,7 +238,7 @@ void too_many_constraint_warnings(int eConstrAlg, int warncount)
 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[], matrix box)
+                             const rvec x[], const matrix box)
 {
     char          fname[STRLEN];
     FILE         *out;
@@ -295,7 +295,7 @@ static void write_constr_pdb(const char *fn, const char *title,
 //! 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[], matrix box)
+                       const rvec x[], rvec xprime[], const matrix box)
 {
     char  buf[STRLEN], buf2[22];
 
@@ -327,7 +327,7 @@ Constraints::apply(bool                  bLog,
                    rvec                 *x,
                    rvec                 *xprime,
                    rvec                 *min_proj,
-                   matrix                box,
+                   const matrix          box,
                    real                  lambda,
                    real                 *dvdlambda,
                    rvec                 *v,
@@ -359,7 +359,7 @@ Constraints::Impl::apply(bool                  bLog,
                          rvec                 *x,
                          rvec                 *xprime,
                          rvec                 *min_proj,
-                         matrix                box,
+                         const matrix          box,
                          real                  lambda,
                          real                 *dvdlambda,
                          rvec                 *v,
@@ -1273,7 +1273,11 @@ bool inter_charge_group_settles(const gmx_mtop_t &mtop)
 
 void do_constrain_first(FILE *fplog, gmx::Constraints *constr,
                         const t_inputrec *ir, const t_mdatoms *md,
-                        t_state *state)
+                        int natoms,
+                        ArrayRefWithPadding<RVec> x,
+                        ArrayRefWithPadding<RVec> v,
+                        const matrix box,
+                        real lambda)
 {
     int             i, m, start, end;
     int64_t         step;
@@ -1281,10 +1285,13 @@ void do_constrain_first(FILE *fplog, gmx::Constraints *constr,
     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, state->natoms + 1);
+    snew(savex, natoms + 1);
 
     start = 0;
     end   = md->homenr;
@@ -1307,9 +1314,9 @@ void do_constrain_first(FILE *fplog, gmx::Constraints *constr,
     /* constrain the current position */
     constr->apply(TRUE, FALSE,
                   step, 0, 1.0,
-                  state->x.rvec_array(), state->x.rvec_array(), nullptr,
-                  state->box,
-                  state->lambda[efptBONDED], &dvdl_dum,
+                  xRvec, xRvec, nullptr,
+                  box,
+                  lambda, &dvdl_dum,
                   nullptr, nullptr, gmx::ConstraintVariable::Positions);
     if (EI_VV(ir->eI))
     {
@@ -1317,24 +1324,24 @@ void do_constrain_first(FILE *fplog, gmx::Constraints *constr,
         /* also may be useful if we need the ekin from the halfstep for velocity verlet */
         constr->apply(TRUE, FALSE,
                       step, 0, 1.0,
-                      state->x.rvec_array(), state->v.rvec_array(), state->v.rvec_array(),
-                      state->box,
-                      state->lambda[efptBONDED], &dvdl_dum,
+                      xRvec, vRvec, vRvec,
+                      box,
+                      lambda, &dvdl_dum,
                       nullptr, nullptr, gmx::ConstraintVariable::Velocities);
     }
     /* constrain the inital velocities at t-dt/2 */
     if (EI_STATE_VELOCITY(ir->eI) && ir->eI != eiVV)
     {
-        auto x = makeArrayRef(state->x).subArray(start, end);
-        auto v = makeArrayRef(state->v).subArray(start, end);
+        auto subX = x.paddedArrayRef().subArray(start, end);
+        auto subV = v.paddedArrayRef().subArray(start, end);
         for (i = start; (i < end); i++)
         {
             for (m = 0; (m < DIM); m++)
             {
                 /* Reverse the velocity */
-                v[i][m] = -v[i][m];
+                subV[i][m] = -subV[i][m];
                 /* Store the position at t-dt in buf */
-                savex[i][m] = x[i][m] + dt*v[i][m];
+                savex[i][m] = subX[i][m] + dt*subV[i][m];
             }
         }
         /* Shake the positions at t=-dt with the positions at t=0
@@ -1349,17 +1356,17 @@ void do_constrain_first(FILE *fplog, gmx::Constraints *constr,
         dvdl_dum = 0;
         constr->apply(TRUE, FALSE,
                       step, -1, 1.0,
-                      state->x.rvec_array(), savex, nullptr,
-                      state->box,
-                      state->lambda[efptBONDED], &dvdl_dum,
-                      state->v.rvec_array(), nullptr, gmx::ConstraintVariable::Positions);
+                      xRvec, savex, nullptr,
+                      box,
+                      lambda, &dvdl_dum,
+                      vRvec, nullptr, gmx::ConstraintVariable::Positions);
 
         for (i = start; i < end; i++)
         {
             for (m = 0; m < DIM; m++)
             {
                 /* Re-reverse the velocities */
-                v[i][m] = -v[i][m];
+                subV[i][m] = -subV[i][m];
             }
         }
     }

diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h
index 452d2ac36d4d6d9d01259731a99d02c2f8ee1768..8fd7f3003d82395d94f7533dd529107012171d4e 100644 (file)
--- a/src/gromacs/mdlib/constr.h
+++ b/src/gromacs/mdlib/constr.h
@@ -74,6 +74,7 @@ class t_state;
 
 namespace gmx
 {
+template <typename T> class ArrayRefWithPadding;
 
 //! Describes supported flavours of constrained updates.
 enum class ConstraintVariable : int
@@ -172,7 +173,7 @@ class Constraints
                    rvec                 *x,
                    rvec                 *xprime,
                    rvec                 *min_proj,
-                   matrix                box,
+                   const matrix          box,
                    real                  lambda,
                    real                 *dvdlambda,
                    rvec                 *v,
@@ -295,7 +296,11 @@ bool inter_charge_group_settles(const gmx_mtop_t &mtop);
 /*! \brief Constrain the initial coordinates and velocities */
 void do_constrain_first(FILE *log, gmx::Constraints *constr,
                         const t_inputrec *inputrec, const t_mdatoms *md,
-                        t_state *state);
+                        int natoms,
+                        ArrayRefWithPadding<RVec> x,
+                        ArrayRefWithPadding<RVec> v,
+                        const matrix box,
+                        real lambda);
 
 }  // namespace gmx
 

diff --git a/src/gromacs/mdlib/lincs.cpp b/src/gromacs/mdlib/lincs.cpp
index 0b6243a79a237711928e165c94fa4980cff0ab36..149c87590cb2200adb37c886e61167c2e4d42fbb 100644 (file)
--- a/src/gromacs/mdlib/lincs.cpp
+++ b/src/gromacs/mdlib/lincs.cpp
@@ -952,7 +952,7 @@ calc_dist_iter_simd(int                           b0,
 #endif // GMX_SIMD_HAVE_REAL
 
 //! Implements LINCS constraining.
-static void do_lincs(const rvec *x, rvec *xp, matrix box, t_pbc *pbc,
+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,
@@ -2296,7 +2296,7 @@ bool constrain_lincs(bool computeRmsd,
                      const t_commrec *cr,
                      const gmx_multisim_t *ms,
                      const rvec *x, rvec *xprime, rvec *min_proj,
-                     matrix box, t_pbc *pbc,
+                     const matrix box, t_pbc *pbc,
                      real lambda, real *dvdlambda,
                      real invdt, rvec *v,
                      bool bCalcVir, tensor vir_r_m_dr,

diff --git a/src/gromacs/mdlib/lincs.h b/src/gromacs/mdlib/lincs.h
index c04741e45fd10c53c205c02f158ecc4ae5e9b103..3405df1648538060678d41edd68fe982c42238a6 100644 (file)
--- a/src/gromacs/mdlib/lincs.h
+++ b/src/gromacs/mdlib/lincs.h
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2018, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019, 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.
@@ -99,7 +99,7 @@ constrain_lincs(bool computeRmsd,
                 const t_commrec *cr,
                 const gmx_multisim_t *ms,
                 const rvec *x, rvec *xprime, rvec *min_proj,
-                matrix box, t_pbc *pbc,
+                const matrix box, t_pbc *pbc,
                 real lambda, real *dvdlambda,
                 real invdt, rvec *v,
                 bool bCalcVir, tensor vir_r_m_dr,

diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp
index 18dc39041dc83743567e2a93e8a0d559e4fe2b82..9e33cc02ddf7b68e8005185f3e94344f863a7738 100644 (file)
--- a/src/gromacs/mdrun/md.cpp
+++ b/src/gromacs/mdrun/md.cpp
@@ -454,7 +454,11 @@ void gmx::LegacySimulator::do_md()
         if (constr)
         {
             /* Constrain the initial coordinates and velocities */
-            do_constrain_first(fplog, constr, ir, mdatoms, state);
+            do_constrain_first(fplog, constr, ir, mdatoms,
+                               state->natoms,
+                               state->x.arrayRefWithPadding(),
+                               state->v.arrayRefWithPadding(),
+                               state->box, state->lambda[efptBONDED]);
         }
         if (vsite)
         {