Clear vsite velocities for simple integrators

The simple integrator loops (introduced in 69470fc4) do not clear
the velocities of virtual sites. This allows velocities of virtual
sites to slowly increase over time. To prevent this, velocities
of virtual sites are now cleared in a separate loop.

Fixes #2316

Change-Id: I12ff0fae2cd3c45ad4e63bfeccfc8c88505cdb1e

diff --git a/src/gromacs/mdlib/mdatoms.cpp b/src/gromacs/mdlib/mdatoms.cpp
index fe998ff9c53c1811ce9a6ef99fd091495aefefe0..b9633373c6edb248c65c0972bf12254216376d58 100644 (file)
--- a/src/gromacs/mdlib/mdatoms.cpp
+++ b/src/gromacs/mdlib/mdatoms.cpp
@@ -95,6 +95,7 @@ makeMDAtoms(FILE *fp, const gmx_mtop_t &mtop, const t_inputrec &ir,
     double                     totalMassA = 0.0;
     double                     totalMassB = 0.0;
 
+    md->haveVsites = FALSE;
     gmx_mtop_atomloop_block_t  aloop = gmx_mtop_atomloop_block_init(&mtop);
     const t_atom              *atom;
     int                        nmol;
@@ -103,6 +104,11 @@ makeMDAtoms(FILE *fp, const gmx_mtop_t &mtop, const t_inputrec &ir,
         totalMassA += nmol*atom->m;
         totalMassB += nmol*atom->mB;
 
+        if (atom->ptype == eptVSite)
+        {
+            md->haveVsites = TRUE;
+        }
+
         if (ir.efep != efepNO && PERTURBED(*atom))
         {
             md->nPerturbed++;

diff --git a/src/gromacs/mdlib/update.cpp b/src/gromacs/mdlib/update.cpp
index 41a24e62f841969873b3d423b2974de0a5968162..6fc8aa4bbf8e0f181b1f8bf9e563e6b55a721af6 100644 (file)
--- a/src/gromacs/mdlib/update.cpp
+++ b/src/gromacs/mdlib/update.cpp
@@ -154,6 +154,21 @@ static bool isPressureCouplingStep(gmx_int64_t step, const t_inputrec *ir)
             do_per_step(step + ir->nstpcouple - offset, ir->nstpcouple));
 }
 
+/*! \brief Sets the velocities of virtual sites to zero */
+static void clearVsiteVelocities(int                   start,
+                                 int                   nrend,
+                                 const unsigned short *particleType,
+                                 rvec * gmx_restrict   v)
+{
+    for (int a = start; a < nrend; a++)
+    {
+        if (particleType[a] == eptVSite)
+        {
+            clear_rvec(v[a]);
+        }
+    }
+}
+
 /*! \brief Sets the number of different temperature coupling values */
 enum class NumTempScaleValues
 {
@@ -555,6 +570,19 @@ static void do_update_md(int                         start,
     }
     else
     {
+        /* Use a simple and thus more efficient integration loop. */
+        /* The simple loop does not check for particle type (so it can
+         * be vectorized), which means we need to clear the velocities
+         * of virtual sites in advance, when present. Note that vsite
+         * velocities are computed after update and constraints from
+         * their displacement.
+         */
+        if (md->haveVsites)
+        {
+            /* Note: The overhead of this loop is completely neligible */
+            clearVsiteVelocities(start, nrend, md->ptype, v);
+        }
+
         /* We determine if we have a single T-coupling lambda value for all
          * atoms. That allows for better SIMD acceleration in the template.
          * If we do not do temperature coupling (in the run or this step),

diff --git a/src/gromacs/mdtypes/mdatom.h b/src/gromacs/mdtypes/mdatom.h
index f2f28986f1c953061a352d9bc5bfaa56034696ce..41fbe70e4bc91d3ddadac86d8c35385f2e53eb2f 100644 (file)
--- a/src/gromacs/mdtypes/mdatom.h
+++ b/src/gromacs/mdtypes/mdatom.h
@@ -62,6 +62,8 @@ typedef struct t_mdatoms {
     int                    nenergrp;
     //! Do we have multiple center of mass motion removal groups
     gmx_bool               bVCMgrps;
+    //! Do we have any virtual sites?
+    gmx_bool               haveVsites;
     //! Do we have atoms that are frozen along 1 or 2 (not 3) dimensions?
     gmx_bool               havePartiallyFrozenAtoms;
     //! Number of perturbed atoms