Move virtual site calculation

The previous implementation would lead to slight inaccuracies when
virtual positions were calculated before shifting positions in
pressure coupling algorithms. This change moves the calculation to the
top of the loop, right before DD and force calculation.

This also removes an unnecessary call to the virtual site construction
in rerun. That call was placed after force calculation and trajectory
writing, so had no influence since coordinates for the next step are
read from file. Since virtual sites are constructed when reading
coordinates from file, the results should have been correct anyway, the
call just had no influence.

This also moves the call to the virtual site construction in MiMiC.
Similar to the rerun case, the current call happened after virtual
sites are used. This call has been moved to right after the coordinates
are received. This is likely to have lead to wrong results when using
MiMiC with virtual sites.

Refs #3866

diff --git a/docs/release-notes/2022/major/bugs-fixed.rst b/docs/release-notes/2022/major/bugs-fixed.rst
index 32793b8aaffed66dfd7129fcb5f05f04709d8062..6a68e9674ff591d47fda2316dd49bf26138eed99 100644 (file)
--- a/docs/release-notes/2022/major/bugs-fixed.rst
+++ b/docs/release-notes/2022/major/bugs-fixed.rst
@@ -1,6 +1,18 @@
 Bugs fixed
 ^^^^^^^^^^
 
+Fixed slight inaccuracies when using virtual sites with pressure coupling
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+Virtual sites were reconstructed after the system was propagated, but before
+scaling due to pressure coupling. For virtual site types which are not a linear
+combination of other atoms, this is not completely correct. Since the scaling
+due to pressure coupling is very small in healthy simulations, the resulting
+inaccuracies are expected to have been extremely minor, and in most cases
+undetectable.
+
+:issue:`3866`
+
 .. Note to developers!
    Please use """"""" to underline the individual entries for fixed issues in the subfolders,
    otherwise the formatting on the webpage is messed up.

diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp
index 2cc3dad638ce910c6cc6eca40df306098e29e4ab..793dfc37c2055045e346f6d5f1b71b80e88afe3e 100644 (file)
--- a/src/gromacs/mdrun/md.cpp
+++ b/src/gromacs/mdrun/md.cpp
@@ -581,11 +581,6 @@ void gmx::LegacySimulator::do_md()
                                state->box,
                                state->lambda[efptBONDED]);
         }
-        if (vsite)
-        {
-            /* Construct the virtual sites for the initial configuration */
-            vsite->construct(state->x, ir->delta_t, {}, state->box);
-        }
     }
 
     if (ir->efep != efepNO)
@@ -954,6 +949,14 @@ void gmx::LegacySimulator::do_md()
             stateGpu->waitCoordinatesReadyOnHost(AtomLocality::Local);
         }
 
+        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);
+            wallcycle_stop(wcycle, ewcVSITECONSTR);
+        }
+
         if (bNS && !(bFirstStep && ir->bContinuation))
         {
             bMasterState = FALSE;
@@ -1564,13 +1567,6 @@ void gmx::LegacySimulator::do_md()
             enerd->term[F_DVDL_CONSTR] += dvdl_constr;
         }
 
-        if (vsite != nullptr)
-        {
-            wallcycle_start(wcycle, ewcVSITECONSTR);
-            vsite->construct(state->x, ir->delta_t, state->v, state->box);
-            wallcycle_stop(wcycle, ewcVSITECONSTR);
-        }
-
         /* ############## IF NOT VV, Calculate globals HERE  ############ */
         /* With Leap-Frog we can skip compute_globals at
          * non-communication steps, but we need to calculate

diff --git a/src/gromacs/mdrun/mimic.cpp b/src/gromacs/mdrun/mimic.cpp
index 970bb40afe3fd60c5ebc37d3ad06fb0559e8fc8c..bae23f3d25fc59dd125229e42f0e67a5d6ac29c8 100644 (file)
--- a/src/gromacs/mdrun/mimic.cpp
+++ b/src/gromacs/mdrun/mimic.cpp
@@ -469,6 +469,12 @@ void gmx::LegacySimulator::do_mimic()
                           "decomposition, "
                           "use a single rank");
             }
+            if (constructVsites)
+            {
+                wallcycle_start(wcycle, ewcVSITECONSTR);
+                vsite->construct(state->x, ir->delta_t, state->v, state->box);
+                wallcycle_stop(wcycle, ewcVSITECONSTR);
+            }
         }
 
         if (DOMAINDECOMP(cr))
@@ -621,13 +627,6 @@ void gmx::LegacySimulator::do_mimic()
 
         stopHandler->setSignal();
 
-        if (vsite != nullptr)
-        {
-            wallcycle_start(wcycle, ewcVSITECONSTR);
-            vsite->construct(state->x, ir->delta_t, state->v, state->box);
-            wallcycle_stop(wcycle, ewcVSITECONSTR);
-        }
-
         {
             const bool          doInterSimSignal = false;
             const bool          doIntraSimSignal = true;

diff --git a/src/gromacs/mdrun/rerun.cpp b/src/gromacs/mdrun/rerun.cpp
index 0857b5b7953dc34e94c83804e7a48c33d0cbbd15..d80761bc950204180bed282f6184b13259316d4e 100644 (file)
--- a/src/gromacs/mdrun/rerun.cpp
+++ b/src/gromacs/mdrun/rerun.cpp
@@ -739,13 +739,6 @@ void gmx::LegacySimulator::do_rerun()
 
         stopHandler->setSignal();
 
-        if (vsite != nullptr)
-        {
-            wallcycle_start(wcycle, ewcVSITECONSTR);
-            vsite->construct(state->x, ir->delta_t, state->v, state->box);
-            wallcycle_stop(wcycle, ewcVSITECONSTR);
-        }
-
         {
             const bool          doInterSimSignal = false;
             const bool          doIntraSimSignal = true;