From 2e4d9cb3aae2164176291dfb30172a534dde45c6 Mon Sep 17 00:00:00 2001
From: Pascal Merz <pascal.merz@me.com>
Date: Fri, 12 Feb 2021 08:02:58 +0000
Subject: [PATCH] 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
---
 docs/release-notes/2022/major/bugs-fixed.rst | 12 ++++++++++++
 src/gromacs/mdrun/md.cpp                     | 20 ++++++++------------
 src/gromacs/mdrun/mimic.cpp                  | 13 ++++++-------
 src/gromacs/mdrun/rerun.cpp                  |  7 -------
 4 files changed, 26 insertions(+), 26 deletions(-)

diff --git a/docs/release-notes/2022/major/bugs-fixed.rst b/docs/release-notes/2022/major/bugs-fixed.rst
index 32793b8aaf..6a68e9674f 100644
--- 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 2cc3dad638..793dfc37c2 100644
--- 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 970bb40afe..bae23f3d25 100644
--- 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 0857b5b795..d80761bc95 100644
--- 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;
-- 
2.22.0