Implementation of a stochastic velocity rescaling thermostat

Bussi, et al. proposed a Gibbs sampling type thermostat which
couples the system to a heat bath by solving an auxiliary
Fokker-Planck equation for the kinetic energy in order to obtain
a (globaly) velocity scaling factor. The solution is to a stochastic
differential equation so the scaling factor depends on the current
kinetic energy as well as a stochastic component.

This algorithm has a two distinct advantages over Langevin dynamics:

1) In general, only 2 random variables are required at each
  rescaling, down from 3*Natoms for Langevin. Even for a box of
  216 TIP3P waters on my laptop, this leads to a ~5% speed-up!

2) Rescaling velocities, as opposed to adding noise to the vectors,
  is much less disruptive of the underlying dynamics and leaves
  holonomic constraints intact. One can thus avoid over-application
  of expensive RATTLE calls during the Ornstein-Uhlenbeck steps.

Possible issues with the present implementation:

1) The rescaling routines currently assume that the number of
  degrees of freedom is very large such that the underlying Gamma
  distribution can be treated as a Gaussian distribution. I haven't
  investigated this in-depth, but clearly it works for O(10^3)
  degrees of freedom.

2) The equations in the paper utilize a time parameter, which we
  call the stochVelRescalingPeriod, however, this could easily
  be converted to a frequency and save a division operation every
  few steps. Furthermore, using the standard NAMD idiom, the number
  of steps between rescalings is the stochVelRescalingFreq, even
  though that is also technically a period (in timesteps). It's
  debatable if this is a nomenclature we would like to encourage.

3) We don't currently permit different time factors per atom, as
  tCouple and Langevin do. It's not clear how recommendable this
  might be and it will change the implementation and performance,
  but other codes support it.

4) There are currently no defaults, although it would be quite
  reasonable to suggest

  stochVelRescalingPeriod 1.0 ;# in ps, suitable for water
  stochVelRescalingFreq   10  ;# the default pairlist frequency

5) It is generally straightforward to compute heat transfer during
  the rescaling, which can be (nearly trivially) used to compute
  the "shadow work" performed by the integrator (a quantity related
  to integrator error). This may be a very easy and general way to
  use this integrator for nonequilibrium MD/MC, but it is not done
  at present.

Change-Id: I6fbea09971a41f46f6f670977f1deb30c06a427d

diff --git a/src/Broadcasts.h b/src/Broadcasts.h
index 4740495..cc644d1 100644 (file)
--- a/src/Broadcasts.h
+++ b/src/Broadcasts.h
@@ -45,6 +45,7 @@ enum {
   positionRescaleFactor2Tag,
   // End multigrator
   tcoupleCoefficientTag,
+  stochRescaleCoefficientTag,
   minimizeCoefficientTag,
   momentumCorrectionTag,
 #if USE_BARRIER
@@ -72,6 +73,7 @@ struct ControllerBroadcasts
   SimpleBroadcastObject<Tensor> positionRescaleFactor2;
   // End multigrator
   SimpleBroadcastObject<BigReal> tcoupleCoefficient;
+  SimpleBroadcastObject<BigReal> stochRescaleCoefficient;
   SimpleBroadcastObject<BigReal> minimizeCoefficient;
   SimpleBroadcastObject<Vector> momentumCorrection;
 #if USE_BARRIER
@@ -80,12 +82,12 @@ struct ControllerBroadcasts
   SimpleBroadcastObject<int> scriptBarrier;
   SimpleBroadcastObject<int> traceBarrier;
   SimpleBroadcastObject<Vector> accelMDRescaleFactor;
-  SimpleBroadcastObject<BigReal> adaptTemperature; 
+  SimpleBroadcastObject<BigReal> adaptTemperature;
 #ifdef MEASURE_NAMD_WITH_PAPI
   SimpleBroadcastObject<int> papiMeasureBarrier;
 #endif
 
-  ControllerBroadcasts(const LDObjHandle *ldObjPtr = 0) : 
+  ControllerBroadcasts(const LDObjHandle *ldObjPtr = 0) :
     velocityRescaleFactor(velocityRescaleFactorTag, ldObjPtr),
     positionRescaleFactor(positionRescaleFactorTag, ldObjPtr),
     // For multigrator
@@ -95,13 +97,14 @@ struct ControllerBroadcasts
     positionRescaleFactor2(positionRescaleFactor2Tag, ldObjPtr),
     // End multigrator
     tcoupleCoefficient(tcoupleCoefficientTag, ldObjPtr),
+    stochRescaleCoefficient(stochRescaleCoefficientTag, ldObjPtr),
     minimizeCoefficient(minimizeCoefficientTag, ldObjPtr),
     momentumCorrection(momentumCorrectionTag, ldObjPtr),
 #if USE_BARRIER
     cycleBarrier(cycleBarrierTag, ldObjPtr),
 #endif
     accelMDRescaleFactor(accelMDRescaleFactorTag, ldObjPtr),
-    adaptTemperature(adaptTemperatureTag, ldObjPtr), 
+    adaptTemperature(adaptTemperatureTag, ldObjPtr),
     scriptBarrier(scriptBarrierTag, ldObjPtr),
 #ifdef MEASURE_NAMD_WITH_PAPI
        papiMeasureBarrier(papiMeasureTag, ldObjPtr),

diff --git a/src/Controller.C b/src/Controller.C
index 1cea978..be93f0e 100644 (file)
--- a/src/Controller.C
+++ b/src/Controller.C
@@ -183,7 +183,13 @@ Controller::Controller(NamdState *s) :
     random = new Random(simParams->randomSeed);
     random->split(0,PatchMap::Object()->numPatches()+1);
 
-    rescaleVelocities_sumTemps = 0;  rescaleVelocities_numTemps = 0;
+    rescaleVelocities_sumTemps = 0;  
+    rescaleVelocities_numTemps = 0;
+    if (simParams->stochRescaleOn) {
+      stochRescale_count = 0;
+      stochRescaleTimefactor =\
+        exp(-simParams->stochRescaleFreq*simParams->dt*0.001/simParams->stochRescalePeriod);
+    }
     berendsenPressure_avg = 0; berendsenPressure_count = 0;
     // strainRate tensor is symmetric to avoid rotation
     langevinPiston_strainRate =
@@ -468,6 +474,7 @@ void Controller::integrate(int scriptTask) {
         adaptTempUpdate(step);
         rescaleVelocities(step);
        tcoupleVelocities(step);
+        stochRescaleVelocities(step);
        berendsenPressure(step);
        langevinPiston1(step);
         rescaleaccelMD(step);
@@ -1318,6 +1325,39 @@ void Controller::tcoupleVelocities(int step)
   }
 }
 
+/*
+* Generate random numbers and broadcast the velocity rescaling factor.
+*/
+void Controller::stochRescaleVelocities(int step)
+{
+  if ( simParams->stochRescaleOn )
+  {
+    ++stochRescale_count;
+    if ( stochRescale_count == simParams->stochRescaleFreq )
+    { 
+      const BigReal stochRescaleTemp = simParams->stochRescaleTemp;
+
+      BigReal coefficient = 1.;
+      if ( temperature > 0. ) 
+      {
+        BigReal R1 = random->gaussian();
+        BigReal gammaShape = 0.5*(numDegFreedom - 1);
+        // R2sum is the sum of (numDegFreedom - 1) squared normal variables, which is
+        // chi-squared distributed. This is in turn a special case of the Gamma
+        // distribution, which converges to a normal distribution in the limit of a
+        // large shape parameter.
+        BigReal R2sum = 2*(gammaShape + sqrt(gammaShape)*random->gaussian()) + R1*R1;
+        BigReal tempfactor = stochRescaleTemp/(temperature*numDegFreedom);
+
+        coefficient = sqrt(stochRescaleTimefactor + (1 - stochRescaleTimefactor)*tempfactor*R2sum
+                  + 2*R1*sqrt(tempfactor*(1 - stochRescaleTimefactor)*stochRescaleTimefactor)); 
+      }
+      broadcast->stochRescaleCoefficient.publish(step,coefficient);
+      stochRescale_count = 0;
+    }
+  }
+}
+
 static char *FORMAT(BigReal X)
 {
   static char tmp_string[25];

diff --git a/src/Controller.h b/src/Controller.h
index 57bfe07..49d1f45 100644 (file)
--- a/src/Controller.h
+++ b/src/Controller.h
@@ -171,6 +171,9 @@ protected:
       int rescaleVelocities_numTemps;
     void reassignVelocities(int);
     void tcoupleVelocities(int);
+    void stochRescaleVelocities(int);
+    int stochRescale_count; 
+    BigReal stochRescaleTimefactor;
     void berendsenPressure(int);
       // Tensor berendsenPressure_avg;
       // int berendsenPressure_count;

diff --git a/src/Sequencer.C b/src/Sequencer.C
index d159ff6..b0264f5 100644 (file)
--- a/src/Sequencer.C
+++ b/src/Sequencer.C
@@ -72,6 +72,7 @@ Sequencer::Sequencer(HomePatch *p) :
     random->split(patch->getPatchID()+1,PatchMap::Object()->numPatches()+1);
 
     rescaleVelocities_numTemps = 0;
+    stochRescale_count = 0;
     berendsenPressure_count = 0;
 //    patch->write_tip4_props();
 }
@@ -346,6 +347,7 @@ void Sequencer::integrate(int scriptTask) {
 
       rescaleVelocities(step);
       tcoupleVelocities(timestep,step);
+      stochRescaleVelocities(timestep,step);
       berendsenPressure(step);
 
       if ( ! commOnly ) {
@@ -1671,6 +1673,30 @@ void Sequencer::tcoupleVelocities(BigReal dt_fs, int step)
   }
 }
 
+/*
+* Rescale velocities with the scale factor sent from the corresponding Controller method.
+*/
+void Sequencer::stochRescaleVelocities(BigReal dt_fs, int step)
+{
+  if ( simParams->stochRescaleOn )
+  {
+    FullAtom *a = patch->atom.begin();
+    int numAtoms = patch->numAtoms;
+    ++stochRescale_count;
+    if ( stochRescale_count == simParams->stochRescaleFreq )
+    {
+      // Blocking receive for the temperature coupling coefficient.
+      BigReal coefficient = broadcast->stochRescaleCoefficient.get(step);
+      DebugM(4, "stochastically rescaling velocities at step " << step << " by " << coefficient << "\n");
+      for ( int i = 0; i < numAtoms; ++i )
+      {
+        a[i].velocity *= coefficient;
+      }
+      stochRescale_count = 0;
+    }
+  }
+}
+
 void Sequencer::saveForce(const int ftag)
 {
   patch->saveForce(ftag);

diff --git a/src/Sequencer.h b/src/Sequencer.h
index f781b13..b02356b 100644 (file)
--- a/src/Sequencer.h
+++ b/src/Sequencer.h
@@ -88,6 +88,8 @@ protected:
     void reinitVelocities(void);
     void rescaleVelocitiesByFactor(BigReal);
     void tcoupleVelocities(BigReal,int);
+    void stochRescaleVelocities(BigReal,int);
+    int stochRescale_count;
     void berendsenPressure(int);
       int berendsenPressure_count;
       int checkpoint_berendsenPressure_count;

diff --git a/src/SimParameters.C b/src/SimParameters.C
index ab76a45..506135d 100644 (file)
--- a/src/SimParameters.C
+++ b/src/SimParameters.C
@@ -1208,6 +1208,23 @@ void SimParameters::config_parser_methods(ParseOptions &opts) {
      "containing the temperature coupling term B(i);\n"
      "default is 'O'", PARSE_STRING);
 
+   opts.optionalB("main", "stochRescale",
+      "Should stochastic velocity rescaling be performed?",
+      &stochRescaleOn, FALSE);
+   opts.require("stochRescale", "stochRescaleTemp",
+      "Temperature for stochastic velocity rescaling",
+       &stochRescaleTemp);
+   opts.range("stochRescaleTemp", NOT_NEGATIVE);
+   opts.units("stochRescaleTemp", N_KELVIN);
+   opts.require("stochRescale", "stochRescalePeriod",
+      "Time scale for stochastic velocity rescaling (ps)",
+       &stochRescalePeriod);
+   opts.range("stochRescalePeriod", POSITIVE);
+   opts.require("stochRescale", "stochRescaleFreq",
+       "Number of steps between stochastic rescalings",
+        &stochRescaleFreq);
+   opts.range("stochRescaleFreq", POSITIVE);
+
    opts.optional("main", "rescaleFreq", "Number of steps between "
     "velocity rescaling", &rescaleFreq);
    opts.range("rescaleFreq", POSITIVE);

diff --git a/src/SimParameters.h b/src/SimParameters.h
index 8839efb..3090d4b 100644 (file)
--- a/src/SimParameters.h
+++ b/src/SimParameters.h
@@ -562,6 +562,11 @@ public:
                                        //  active
        BigReal tCoupleTemp;            //  Temperature for temp coupling
 
+        Bool stochRescaleOn;            //  Flag TRUE-> stochastic velocity rescaling
+        BigReal stochRescaleTemp;       //  Temperature for stochastic velocity rescaling
+        BigReal stochRescalePeriod;     //  Timescale (ps) for stochastic velocity rescaling
+        int stochRescaleFreq;           //  How frequently stochastic velocity rescaling occurs
+
        int rescaleFreq;                //  Velocity rescale frequency
        BigReal rescaleTemp;            //  Temperature to rescale to