Moving do_constraints_first(...) to constr.h

This is another patch in the cleaning efforts of sim_util.

The idea is to only have functions relevant to the force
schedules there so that it becomes easy to move it into
its own module with minimal merging pains.

Change-Id: Ic311bd9afe8baddc0589f368ebf464f1699589dc
Related: #2574

diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp
index a888364b74c1cacb82071ca0923b72f65635c278..a56a3ff6de8979fee002799e384ffd9aa9267cca 100644 (file)
--- a/src/gromacs/mdlib/constr.cpp
+++ b/src/gromacs/mdlib/constr.cpp
@@ -1264,4 +1264,99 @@ bool inter_charge_group_settles(const gmx_mtop_t &mtop)
     return bInterCG;
 }
 
+void do_constrain_first(FILE *fplog, gmx::Constraints *constr,
+                        const t_inputrec *ir, const t_mdatoms *md,
+                        t_state *state)
+{
+    int             i, m, start, end;
+    int64_t         step;
+    real            dt = ir->delta_t;
+    real            dvdl_dum;
+    rvec           *savex;
+
+    /* 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);
+
+    start = 0;
+    end   = md->homenr;
+
+    if (debug)
+    {
+        fprintf(debug, "vcm: start=%d, homenr=%d, end=%d\n",
+                start, md->homenr, end);
+    }
+    /* Do a first constrain to reset particles... */
+    step = ir->init_step;
+    if (fplog)
+    {
+        char buf[STEPSTRSIZE];
+        fprintf(fplog, "\nConstraining the starting coordinates (step %s)\n",
+                gmx_step_str(step, buf));
+    }
+    dvdl_dum = 0;
+
+    /* 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,
+                  nullptr, 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,
+                      state->x.rvec_array(), state->v.rvec_array(), state->v.rvec_array(),
+                      state->box,
+                      state->lambda[efptBONDED], &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);
+        for (i = start; (i < end); i++)
+        {
+            for (m = 0; (m < DIM); m++)
+            {
+                /* Reverse the velocity */
+                v[i][m] = -v[i][m];
+                /* Store the position at t-dt in buf */
+                savex[i][m] = x[i][m] + dt*v[i][m];
+            }
+        }
+        /* Shake the positions at t=-dt with the positions at t=0
+         * as reference coordinates.
+         */
+        if (fplog)
+        {
+            char buf[STEPSTRSIZE];
+            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,
+                      state->x.rvec_array(), savex, nullptr,
+                      state->box,
+                      state->lambda[efptBONDED], &dvdl_dum,
+                      state->v.rvec_array(), 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];
+            }
+        }
+    }
+    sfree(savex);
+}
+
 }  // namespace gmx

diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h
index 147ad445615f67e1a461b5adf9e3837c147b818b..e0104e73b56b6ba5adc527d2db3671cb8fbbe9cb 100644 (file)
--- a/src/gromacs/mdlib/constr.h
+++ b/src/gromacs/mdlib/constr.h
@@ -290,6 +290,11 @@ bool inter_charge_group_constraints(const gmx_mtop_t &mtop);
 /*! \brief Returns whether there are inter charge group settles */
 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);
+
 }  // namespace gmx
 
 #endif

diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp
index 5c76eeb8c2de0361cb112b172cbac20b57e73e46..6c5d3995d7770dc787e64788bb70eea2eabd29e1 100644 (file)
--- a/src/gromacs/mdlib/sim_util.cpp
+++ b/src/gromacs/mdlib/sim_util.cpp
@@ -1883,102 +1883,6 @@ void do_force(FILE                                     *fplog,
     ddBalanceRegionHandler.openBeforeForceComputationCpu(DdAllowBalanceRegionReopen::no);
 }
 
-
-void do_constrain_first(FILE *fplog, gmx::Constraints *constr,
-                        const t_inputrec *ir, const t_mdatoms *md,
-                        t_state *state)
-{
-    int             i, m, start, end;
-    int64_t         step;
-    real            dt = ir->delta_t;
-    real            dvdl_dum;
-    rvec           *savex;
-
-    /* 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);
-
-    start = 0;
-    end   = md->homenr;
-
-    if (debug)
-    {
-        fprintf(debug, "vcm: start=%d, homenr=%d, end=%d\n",
-                start, md->homenr, end);
-    }
-    /* Do a first constrain to reset particles... */
-    step = ir->init_step;
-    if (fplog)
-    {
-        char buf[STEPSTRSIZE];
-        fprintf(fplog, "\nConstraining the starting coordinates (step %s)\n",
-                gmx_step_str(step, buf));
-    }
-    dvdl_dum = 0;
-
-    /* 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,
-                  nullptr, 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,
-                      state->x.rvec_array(), state->v.rvec_array(), state->v.rvec_array(),
-                      state->box,
-                      state->lambda[efptBONDED], &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);
-        for (i = start; (i < end); i++)
-        {
-            for (m = 0; (m < DIM); m++)
-            {
-                /* Reverse the velocity */
-                v[i][m] = -v[i][m];
-                /* Store the position at t-dt in buf */
-                savex[i][m] = x[i][m] + dt*v[i][m];
-            }
-        }
-        /* Shake the positions at t=-dt with the positions at t=0
-         * as reference coordinates.
-         */
-        if (fplog)
-        {
-            char buf[STEPSTRSIZE];
-            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,
-                      state->x.rvec_array(), savex, nullptr,
-                      state->box,
-                      state->lambda[efptBONDED], &dvdl_dum,
-                      state->v.rvec_array(), 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];
-            }
-        }
-    }
-    sfree(savex);
-}
-
 void put_atoms_in_box_omp(int ePBC, const matrix box, gmx::ArrayRef<gmx::RVec> x)
 {
     int t, nth;

diff --git a/src/gromacs/mdlib/sim_util.h b/src/gromacs/mdlib/sim_util.h
index c44b967fe1745437edd2ac27d9687377a69c5dba..3782c104e9e37ada2d183f0df80b9b998c2d0de0 100644 (file)
--- a/src/gromacs/mdlib/sim_util.h
+++ b/src/gromacs/mdlib/sim_util.h
@@ -113,10 +113,6 @@ void initialize_lambdas(FILE               *fplog,
                         gmx::ArrayRef<real> lambda,
                         double             *lam0);
 
-void do_constrain_first(FILE *log, gmx::Constraints *constr,
-                        const t_inputrec *inputrec, const t_mdatoms *md,
-                        t_state *state);
-
 /* Routine in sim_util.c */
 
 gmx_bool use_GPU(const nonbonded_verlet_t *nbv);