Implement velocity terms for virtual sites

[alexxy/gromacs.git] / src / gromacs / mdrun / md.cpp

diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp
index ab2d4b77b17c44d251803baea510de6f03b8d529..cf64bf19c6771a63dd4a72434827e3c203b02f97 100644 (file)
--- a/src/gromacs/mdrun/md.cpp
+++ b/src/gromacs/mdrun/md.cpp
@@ -551,7 +551,7 @@ void gmx::LegacySimulator::do_md()
             /* Set the velocities of vsites, shells and frozen atoms to zero */
             for (i = 0; i < mdatoms->homenr; i++)
             {
-                if (mdatoms->ptype[i] == eptVSite || mdatoms->ptype[i] == eptShell)
+                if (mdatoms->ptype[i] == eptShell)
                 {
                     clear_rvec(v[i]);
                 }
@@ -949,11 +949,22 @@ void gmx::LegacySimulator::do_md()
             stateGpu->waitCoordinatesReadyOnHost(AtomLocality::Local);
         }
 
+        // We only need to calculate virtual velocities if we are writing them in the current step
+        const bool needVirtualVelocitiesThisStep =
+                (vsite != nullptr)
+                && (do_per_step(step, ir->nstvout) || checkpointHandler->isCheckpointingStep());
+
         if (vsite != nullptr)
         {
             // Virtual sites need to be updated before domain decomposition and forces are calculated
             wallcycle_start(wcycle, ewcVSITECONSTR);
-            vsite->construct(state->x, ir->delta_t, state->v, state->box);
+            // md-vv calculates virtual velocities once it has full-step real velocities
+            vsite->construct(state->x,
+                             state->v,
+                             state->box,
+                             (!EI_VV(inputrec->eI) && needVirtualVelocitiesThisStep)
+                                     ? VSiteOperation::PositionsAndVelocities
+                                     : VSiteOperation::Positions);
             wallcycle_stop(wcycle, ewcVSITECONSTR);
         }
 
@@ -1242,6 +1253,13 @@ void gmx::LegacySimulator::do_md()
                                  mdlog,
                                  fplog,
                                  wcycle);
+            if (vsite != nullptr && needVirtualVelocitiesThisStep)
+            {
+                // Positions were calculated earlier
+                wallcycle_start(wcycle, ewcVSITECONSTR);
+                vsite->construct(state->x, state->v, state->box, VSiteOperation::Velocities);
+                wallcycle_stop(wcycle, ewcVSITECONSTR);
+            }
         }
 
         /* ########  END FIRST UPDATE STEP  ############## */