Finalize stochastic velocity rescaling implementation 49/4049/4
authorradakb <brian.radak@gmail.com>
Tue, 24 Apr 2018 13:09:37 +0000 (09:09 -0400)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Tue, 24 Apr 2018 21:45:51 +0000 (16:45 -0500)
Here we add additional scripting features and increase robustness
of the stochastic velocity rescaling implementation. More doxygen
style comments are also added.

o stochRescaleFreq is no longer a required option and instead
  defaults to the steps per cycle frequency (default 20). This is
  in line with the default behavior in GROMACS.

o stochRescale and stochRescaleTemp can now be set from a Tcl
  script. This permits, for example, temperature replica exchange
  or neMD/MC with non-thermostatted switches to be implemented.
  For now, stochRescalePeriod CANNOT be reset, because this would
  require more significant changes, much like langevin and tCouple
  have now. I can't envision any use case for this, so I did not
  consider it worth the time.

o The stochastic component now has a more rigorous backend inside
  the Random class. The new routine sum_of_squared_gaussians(),
  should be pretty transparent. Really this is a highly-specialized
  case of a Gamma distribution, but again, I don't know of any
  other uses for such a routine, so I stuck with the specialized
  routine that assume integer argument (numDegFreedom - 1 in
  practical use). I didn't take the time to profile _this_ routine,
  but I translated the whole thing from Tcl and tested them for
  proper statistics and timings. It looks like the approximation
  for large arguments saves ~33% compute time per RNG.

Change-Id: Ic613a22ddc88c50f832fd3e889dc926ba7aece82

src/Controller.C
src/Random.h
src/Sequencer.C
src/SimParameters.C
src/SimParameters.h

index ab8faf3..8ef5d55 100644 (file)
@@ -1325,9 +1325,17 @@ void Controller::tcoupleVelocities(int step)
   }
 }
 
-/*
-* Generate random numbers and broadcast the velocity rescaling factor.
-*/
+/** Generate and broadcast the scale factor for stochastic velocity rescaling.
+ *
+ *  Stochastic velocity rescaling couples the system to a heat bath by globally
+ *  scaling the velocites by a single factor. This factor is chosen based on
+ *  the instantaneous and bath temperatures, a user-defined time scale, and a
+ *  stochastic component linked to the number of degrees of freedom in the
+ *  system. All of this information is combined here and sent to the Sequencer
+ *  for the actual rescaling.
+ *  
+ *  \param step the current timestep
+ */
 void Controller::stochRescaleVelocities(int step)
 {
   if ( simParams->stochRescaleOn )
@@ -1341,12 +1349,13 @@ void Controller::stochRescaleVelocities(int step)
       if ( temperature > 0. ) 
       {
         BigReal R1 = random->gaussian();
-        BigReal gammaShape = 0.5*(numDegFreedom - 1);
+        // 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 R2sum = 2*(gammaShape + sqrt(gammaShape)*random->gaussian()) + R1*R1;
+        BigReal R2sum = random->sum_of_squared_gaussians(numDegFreedom-1);
         BigReal tempfactor = stochRescaleTemp/(temperature*numDegFreedom);
 
         coefficient = sqrt(stochRescaleTimefactor + (1 - stochRescaleTimefactor)*tempfactor*R2sum
index b5745f8..3fce49a 100644 (file)
@@ -135,6 +135,75 @@ public:
     }
   }
 
+ /** \brief Return the sum of i.i.d. squared Gaussians.
+  *
+  *  The sum of k squared standard normal random variables is equivalent to
+  *  drawing a value from a chi-squared distribution with the shape given as
+  *  the number of Gaussians. That is,
+  *
+  *  X ~ sum_n=1^k N(0, 1)^2 ~ chi^2(k).
+  *
+  *  This is in turn a special case of the Gamma distribution with shape k/2
+  *  and scale 2 (we could also use rate = 1/scale)
+  *
+  *  X ~ chi^2(k) ~ Gamma(k/2, 2) ~ 2*Gamma(k/2, 1).
+  *
+  *  The second relation follows from the scaling property of the Gamma
+  *  distribution. Furthermore, when a Gamma distribution has unit scale and a
+  *  large shape, it can be well-approximated by a normal distribution with
+  *  mean and variance equal to the shape. Thus,
+  *
+  *  X ~ 2*Gamma(k/2, 1) ~= 2*N(k/2, sqrt(k/2)).
+  *
+  *  A quick numerical test shows that the mean integrated square error for
+  *  this approximation is <10^-5 for shape >30 and <10^-6 for shape >100.
+  *  We'll be conservative and use the latter cutoff.
+  *
+  *  We thus have three cases for k Gaussians:
+  *
+  *   0 < k <=   2 - just brute force generate and sum the Gaussians
+  *   2 < k <= 200 - use a (slightly modified) version of the Gamma
+  *                  distribution algorithm from Marsaglia and Tsang that is
+  *                  implemented in the GNU Science Library (GSL)
+  *   else         - use a single Gaussian distribution
+  *
+  *   The brute force method is almost certainly the slowest method, even for
+  *   k = 3. The rigorous method takes about 150% as long as the approximate
+  *   method (but this is invariant to the value of k).
+  *
+  *  \param num_gaussians A positive integer number of Gaussians
+  *  \return a random variable equal to the sum of num_gaussians squared
+  *    standard normal variables
+  */
+  BigReal sum_of_squared_gaussians(int num_gaussians) {
+    BigReal z, u, v;
+
+    if (num_gaussians <= 2) {
+        v = 0.0;
+        for(int i=0; i<num_gaussians; ++i) {
+          z = gaussian();
+          v += z*z;
+        }
+        return v;
+    } else if (2 < num_gaussians && num_gaussians <= 200) {
+      const BigReal d = 0.5*num_gaussians - 1./3;
+      const BigReal c = 1. / (3*sqrt(d));
+      const BigReal zmin = -1. / c;
+      do {
+        do {
+          z = gaussian();
+        }
+        while(z <= zmin);
+        u = uniform();
+        v = (1 + c*z); v = v*v*v; // v = (1 + c*z)^3
+      }
+      while(log(u) >= (0.5*z*z + d*(1 - v + log(v))));
+      return 2*d*v;
+    } else {
+      return num_gaussians + sqrt(2*num_gaussians)*gaussian();
+    }
+  }
+
   // return a vector of gaussian random numbers
   Vector gaussian_vector(void) {
     return Vector( gaussian(), gaussian(), gaussian() );
index b0264f5..a37b609 100644 (file)
@@ -1673,9 +1673,11 @@ void Sequencer::tcoupleVelocities(BigReal dt_fs, int step)
   }
 }
 
-/*
-* Rescale velocities with the scale factor sent from the corresponding Controller method.
-*/
+/** Rescale velocities with the scale factor sent from the Controller.
+ *
+ *  \param dt_fs The integration increment, in fs (not currently used)
+ *  \param step The current timestep
+ */
 void Sequencer::stochRescaleVelocities(BigReal dt_fs, int step)
 {
   if ( simParams->stochRescaleOn )
index 506135d..8e08dfa 100644 (file)
@@ -204,6 +204,8 @@ void SimParameters::scriptSet(const char *param, const char *value) {
   SCRIPT_PARSE_FLOAT("langevinTemp",langevinTemp)
   SCRIPT_PARSE_BOOL("langevinBAOAB",langevin_useBAOAB) // [!!] Use the BAOAB integrator or not
   SCRIPT_PARSE_FLOAT("loweAndersenTemp",loweAndersenTemp) // BEGIN LA, END LA
+  SCRIPT_PARSE_BOOL("stochRescale",stochRescaleOn)
+  SCRIPT_PARSE_FLOAT("stochRescaleTemp",stochRescaleTemp)
   SCRIPT_PARSE_FLOAT("initialTemp",initialTemp)
   SCRIPT_PARSE_BOOL("useGroupPressure",useGroupPressure)
   SCRIPT_PARSE_BOOL("useFlexibleCell",useFlexibleCell)
@@ -1220,7 +1222,7 @@ void SimParameters::config_parser_methods(ParseOptions &opts) {
       "Time scale for stochastic velocity rescaling (ps)",
        &stochRescalePeriod);
    opts.range("stochRescalePeriod", POSITIVE);
-   opts.require("stochRescale", "stochRescaleFreq",
+   opts.optional("stochRescale", "stochRescaleFreq",
        "Number of steps between stochastic rescalings",
         &stochRescaleFreq);
    opts.range("stochRescaleFreq", POSITIVE);
@@ -3305,6 +3307,15 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
    }
    // END LA
 
+   // BKR - stochastic velocity rescaling
+   if (stochRescaleOn) {
+     if (langevinOn || loweAndersenOn || tCoupleOn ||
+         opts.defined("rescaleFreq") || opts.defined("reassignFreq"))
+       NAMD_die("Stochastic velocity rescaling is incompatible with other temperature control methods");
+     // This is largely the same default used in GROMACS.
+     if (!opts.defined("stochRescaleFreq")) stochRescaleFreq = stepsPerCycle;
+   }
+
    if (opts.defined("rescaleFreq"))
    {
   if (!opts.defined("rescaleTemp"))
@@ -3424,6 +3435,7 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
      if             (rescaleFreq > 0)  alchTemp = rescaleTemp;
      else if (reassignFreq > 0)        alchTemp = reassignTemp;
      else if (langevinOn)      alchTemp = langevinTemp;
+     else if (stochRescaleOn)  alchTemp = stochRescaleTemp;
      else if (tCoupleOn)       alchTemp = tCoupleTemp;
      else NAMD_die("Alchemical FEP can be performed only in constant temperature simulations\n");
 
@@ -5573,6 +5585,18 @@ if ( openatomOn )
       iout << endi;
    }
 
+   if (stochRescaleOn)
+   {
+      iout << iINFO << "STOCHASTIC RESCALING ACTIVE\n";
+      iout << iINFO << "STOCHASTIC RESCALING TEMPERATURE "
+           << stochRescaleTemp << " K\n";
+      iout << iINFO << "STOCHASTIC RESCALING PERIOD " 
+           << stochRescalePeriod << " PS\n";
+      iout << iINFO << "STOCHASTIC RESCALING WILL OCCUR EVERY "
+           << stochRescaleFreq << " STEPS\n";
+      iout << endi;
+   }
+
    if (minimizeOn)
    {
       iout << iINFO << "OLD STYLE MINIMIZATION ACTIVE\n";
index 1f7700d..c1036c9 100644 (file)
@@ -575,7 +575,8 @@ public:
   BigReal stochRescalePeriod;
   /**< Timescale (ps) for stochastic velocity rescaling. */
 
-        int stochRescaleFreq;           //  How frequently stochastic velocity rescaling occurs
+  int stochRescaleFreq;
+  /**< How frequently (time steps) stochastic velocity rescaling occurs. */
 
        int rescaleFreq;                //  Velocity rescale frequency
        BigReal rescaleTemp;            //  Temperature to rescale to