From 3a9e67dd55949e7daaa666e5bc4de95bcfece2c4 Mon Sep 17 00:00:00 2001
From: Berk Hess <hess@kth.se>
Date: Mon, 23 May 2016 14:42:49 +0200
Subject: [PATCH] Fix bug in v-rescale thermostat & replica exchange

Commit 2d0247f6 made the random normally distributed number for the
v-rescale thermostat constant for a run only depending on ld-seed,
and moved a reset for replica exchange inside the multi-exchange
branch.

Fixes #1968.

Change-Id: Iadc38ccadb99c6c1232756fd595843a01e5f3ce8
---
 src/gromacs/mdlib/coupling.cpp |  7 +++----
 src/programs/mdrun/repl_ex.cpp | 17 ++++++++++++++---
 2 files changed, 17 insertions(+), 7 deletions(-)

diff --git a/src/gromacs/mdlib/coupling.cpp b/src/gromacs/mdlib/coupling.cpp
index 308a8a2187..abafb2b8b0 100644
--- a/src/gromacs/mdlib/coupling.cpp
+++ b/src/gromacs/mdlib/coupling.cpp
@@ -1386,7 +1386,6 @@ real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
 
 
 static real vrescale_sumnoises(real                            nn,
-                               gmx_int64_t                     step,
                                gmx::ThreeFry2x64<>            *rng,
                                gmx::NormalDistribution<real>  *normalDist)
 {
@@ -1399,8 +1398,6 @@ static real vrescale_sumnoises(real                            nn,
     real                           r;
     gmx::GammaDistribution<real>   gammaDist(0.5*nn, 1.0);
 
-    rng->restart(step, 0);
-
     if (nn < 2 + ndeg_tol)
     {
         int  nn_int, i;
@@ -1453,11 +1450,13 @@ static real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut,
         factor = 0.0;
     }
 
+    rng.restart(step, 0);
+
     rr = normalDist(rng);
 
     ekin_new =
         kk +
-        (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, step, &rng, &normalDist) + rr*rr)/ndeg - kk) +
+        (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, &rng, &normalDist) + rr*rr)/ndeg - kk) +
         2.0*rr*std::sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
 
     return ekin_new;
diff --git a/src/programs/mdrun/repl_ex.cpp b/src/programs/mdrun/repl_ex.cpp
index ab12fc4f17..a3d89573f6 100644
--- a/src/programs/mdrun/repl_ex.cpp
+++ b/src/programs/mdrun/repl_ex.cpp
@@ -961,15 +961,20 @@ test_for_replica_exchange(FILE                 *fplog,
         pind[i] = re->ind[i];
     }
 
+    rng.restart( step, 0 );
+
     if (bMultiEx)
     {
         /* multiple random switch exchange */
         int nself = 0;
 
-        rng.restart( step, 0 );
 
         for (i = 0; i < re->nex + nself; i++)
         {
+            // For now this is superfluous, but just in case we ever add more
+            // calls in different branches it is safer to always reset the distribution.
+            uniformNreplDist.reset();
+
             /* randomly select a pair  */
             /* in theory, could reduce this by identifying only which switches had a nonneglibible
                probability of occurring (log p > -100) and only operate on those switches */
@@ -1014,7 +1019,10 @@ test_for_replica_exchange(FILE                 *fplog,
                 {
                     prob[0] = exp(-delta);
                 }
-                /* roll a number to determine if accepted */
+                // roll a number to determine if accepted. For now it is superfluous to
+                // reset, but just in case we ever add more calls in different branches
+                // it is safer to always reset the distribution.
+                uniformRealDist.reset();
                 bEx[0] = uniformRealDist(rng) < prob[0];
             }
             re->prob_sum[0] += prob[0];
@@ -1060,7 +1068,10 @@ test_for_replica_exchange(FILE                 *fplog,
                     {
                         prob[i] = exp(-delta);
                     }
-                    /* roll a number to determine if accepted */
+                    // roll a number to determine if accepted. For now it is superfluous to
+                    // reset, but just in case we ever add more calls in different branches
+                    // it is safer to always reset the distribution.
+                    uniformRealDist.reset();
                     bEx[i] = uniformRealDist(rng) < prob[i];
                 }
                 re->prob_sum[i] += prob[i];
-- 
2.11.4.GIT