Merge pull request #765 from swails/stochastic-forces

Make it so that random forces and integrators generate unique seeds each time a Context is made

docs-source/usersguide/theory.rst

@@ -1379,3 +1379,67 @@ specific types of rules.  They are:
 .. math::
    \mathbf{r}=\mathbf{o}+p_1\mathbf{\hat{x}}+p_2\mathbf{\hat{y}}+p_3\mathbf{\hat{z}}
 ..
+
+Random Numbers with Stochastic Integrators and Forces
+*****************************************************
+
+OpenMM includes many stochastic integrators and forces that make extensive use
+of random numbers. It is impossible to generate truly random numbers on a
+computer like you would with a dice roll or coin flip in real life---instead
+programs rely on pseudo-random number generators (PRNGs) that take some sort of
+initial "seed" value and steps through a sequence of seemingly random numbers.
+
+The exact implementation of the PRNGs is not important (in fact, each platform
+may have its own PRNG whose performance is optimized for that hardware).  What
+*is* important, however, is that the PRNG must generate a uniform distribution
+of random numbers between 0 and 1. Random numbers drawn from this distribution
+can be manipulated to yield random integers in a desired range or even a random
+number from a different type of probability distribution function (e.g., a
+normal distribution).
+
+What this means is that the random numbers used by integrators and forces within
+OpenMM cannot have any discernible pattern to them.  Patterns can be induced in
+PRNGs in two principal ways:
+
+1. The PRNG uses a bad algorithm with a short period.
+2. Two PRNGs are started using the same seed
+
+All PRNG algorithms in common use are periodic---that is their sequence of
+random numbers repeats after a given *period*, defined by the number of "unique"
+random numbers in the repeating pattern.  As long as this period is longer than
+the total number of random numbers your application requires (preferably by
+orders of magnitude), the first problem described above is avoided. All PRNGs
+employed by OpenMM have periods far longer than any current simulation can cycle
+through.
+
+Point two is far more common in biomolecular simulation, and can result in very
+strange artifacts that may be difficult to detect. For example, with Langevin
+dynamics, two simulations that use the same sequence of random numbers appear to
+synchronize in their global movements.\ :cite:`Uberuaga2004`\
+:cite:`Sindhikara2009` It is therefore very important that the stochastic forces
+and integrators in OpenMM generate unique sequences of pseudo-random numbers not
+only within a single simulation, but between two different simulations of the
+same system as well (including any restarts of previous simulations).
+
+Every stochastic force and integrator that does (or could) make use of random
+numbers has two instance methods attached to it: :meth:`getRandomNumberSeed()`
+and :meth:`setRandomNumberSeed(int seed)`. If you set a unique random seed for
+two different simulations (or different forces/integrators if applicable),
+OpenMM guarantees that the generated sequences of random numbers will be
+different (by contrast, no guarantee is made that the same seed will result in
+identical random number sequences).
+
+Since breaking simulations up into pieces and/or running multiple replicates of
+a system to obtain more complete statistics is common practice, a new strategy
+has been employed for OpenMM versions 6.3 and later with the aim of trying to
+ensure that each simulation will be started with a unique random seed. A random
+seed value of 0 (the default) will cause a unique random seed to be generated
+when a new :class:`Context` is instantiated.
+
+Prior to the introduction of this feature, deserializing a serialized
+:class:`System` XML file would result in each stochastic force or integrator
+being assigned the same random seed as the original instance that was
+serialized. If you use a :class:`System` XML file generated by a version of
+OpenMM older than 6.3 to start a new simulation, you should manually set the
+random number seed of each stochastic force or integrator to 0 (or another
+unique value).

openmmapi/include/openmm/AndersenThermostat.h

@@ -35,6 +35,7 @@
 #include "Force.h"
 #include <string>
 #include "internal/windowsExport.h"
+#include "openmm/internal/OSRngSeed.h"
 
 namespace OpenMM {
 
@@ -113,6 +114,10 @@ public:
      * the other hand, no guarantees are made about the behavior of simulations that use the same seed.
      * In particular, Platforms are permitted to use non-deterministic algorithms which produce different
      * results on successive runs, even if those runs were initialized identically.
+     *
+     * If seed is set to 0 (which is the default value assigned), a unique seed is chosen when a Context
+     * is created from this Force. This is done to ensure that each Context receives unique random seeds
+     * without you needing to set them explicitly.
      */
     void setRandomNumberSeed(int seed) {
         randomNumberSeed = seed;

openmmapi/include/openmm/BrownianIntegrator.h

@@ -99,6 +99,10 @@ public:
      * the other hand, no guarantees are made about the behavior of simulations that use the same seed.
      * In particular, Platforms are permitted to use non-deterministic algorithms which produce different
      * results on successive runs, even if those runs were initialized identically.
+     *
+     * If seed is set to 0 (which is the default value assigned), a unique seed is chosen when a Context
+     * is created from this Force. This is done to ensure that each Context receives unique random seeds
+     * without you needing to set them explicitly.
      */
     void setRandomNumberSeed(int seed) {
         randomNumberSeed = seed;

openmmapi/include/openmm/CustomIntegrator.h

@@ -444,6 +444,10 @@ public:
      * the other hand, no guarantees are made about the behavior of simulations that use the same seed.
      * In particular, Platforms are permitted to use non-deterministic algorithms which produce different
      * results on successive runs, even if those runs were initialized identically.
+     *
+     * If seed is set to 0 (which is the default value assigned), a unique seed is chosen when a Context
+     * is created from this Force. This is done to ensure that each Context receives unique random seeds
+     * without you needing to set them explicitly.
      */
     void setRandomNumberSeed(int seed) {
         randomNumberSeed = seed;

openmmapi/include/openmm/LangevinIntegrator.h

@@ -99,6 +99,10 @@ public:
      * the other hand, no guarantees are made about the behavior of simulations that use the same seed.
      * In particular, Platforms are permitted to use non-deterministic algorithms which produce different
      * results on successive runs, even if those runs were initialized identically.
+     *
+     * If seed is set to 0 (which is the default value assigned), a unique seed is chosen when a Context
+     * is created from this Force. This is done to ensure that each Context receives unique random seeds
+     * without you needing to set them explicitly.
      */
     void setRandomNumberSeed(int seed) {
         randomNumberSeed = seed;

openmmapi/include/openmm/MonteCarloAnisotropicBarostat.h

@@ -167,6 +167,10 @@ public:
      * the other hand, no guarantees are made about the behavior of simulations that use the same seed.
      * In particular, Platforms are permitted to use non-deterministic algorithms which produce different
      * results on successive runs, even if those runs were initialized identically.
+     *
+     * If seed is set to 0 (which is the default value assigned), a unique seed is chosen when a Context
+     * is created from this Force. This is done to ensure that each Context receives unique random seeds
+     * without you needing to set them explicitly.
      */
     void setRandomNumberSeed(int seed) {
         randomNumberSeed = seed;

openmmapi/include/openmm/MonteCarloBarostat.h

@@ -123,6 +123,10 @@ public:
      * the other hand, no guarantees are made about the behavior of simulations that use the same seed.
      * In particular, Platforms are permitted to use non-deterministic algorithms which produce different
      * results on successive runs, even if those runs were initialized identically.
+     *
+     * If seed is set to 0 (which is the default value assigned), a unique seed is chosen when a Context
+     * is created from this Force. This is done to ensure that each Context receives unique random seeds
+     * without you needing to set them explicitly.
      */
     void setRandomNumberSeed(int seed) {
         randomNumberSeed = seed;

openmmapi/include/openmm/MonteCarloMembraneBarostat.h

@@ -220,6 +220,10 @@ public:
      * the other hand, no guarantees are made about the behavior of simulations that use the same seed.
      * In particular, Platforms are permitted to use non-deterministic algorithms which produce different
      * results on successive runs, even if those runs were initialized identically.
+     *
+     * If seed is set to 0 (which is the default value assigned), a unique seed is chosen when a Context
+     * is created from this Force. This is done to ensure that each Context receives unique random seeds
+     * without you needing to set them explicitly.
      */
     void setRandomNumberSeed(int seed) {
         randomNumberSeed = seed;

openmmapi/include/openmm/VariableLangevinIntegrator.h

@@ -121,6 +121,10 @@ public:
      * the other hand, no guarantees are made about the behavior of simulations that use the same seed.
      * In particular, Platforms are permitted to use non-deterministic algorithms which produce different
      * results on successive runs, even if those runs were initialized identically.
+     *
+     * If seed is set to 0 (which is the default value assigned), a unique seed is chosen when a Context
+     * is created from this Force. This is done to ensure that each Context receives unique random seeds
+     * without you needing to set them explicitly.
      */
     void setRandomNumberSeed(int seed) {
         randomNumberSeed = seed;

openmmapi/src/AndersenThermostat.cpp

@@ -37,7 +37,7 @@ using namespace OpenMM;
 
 AndersenThermostat::AndersenThermostat(double defaultTemperature, double defaultCollisionFrequency) :
         defaultTemp(defaultTemperature), defaultFreq(defaultCollisionFrequency) {
-    setRandomNumberSeed(osrngseed());
+    setRandomNumberSeed(0);
 }
 
 ForceImpl* AndersenThermostat::createImpl() const {

openmmapi/src/BrownianIntegrator.cpp

@@ -33,7 +33,6 @@
 #include "openmm/Context.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
-#include "openmm/internal/OSRngSeed.h"
 #include "openmm/kernels.h"
 #include <string>
 
@@ -46,7 +45,7 @@ BrownianIntegrator::BrownianIntegrator(double temperature, double frictionCoeff,
     setFriction(frictionCoeff);
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
-    setRandomNumberSeed(osrngseed());
+    setRandomNumberSeed(0);
 }
 
 void BrownianIntegrator::initialize(ContextImpl& contextRef) {

openmmapi/src/CustomIntegrator.cpp

@@ -34,7 +34,6 @@
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/ContextImpl.h"
-#include "openmm/internal/OSRngSeed.h"
 #include "openmm/kernels.h"
 #include <set>
 #include <string>
@@ -45,7 +44,7 @@ using namespace std;
 CustomIntegrator::CustomIntegrator(double stepSize) : globalsAreCurrent(true), forcesAreValid(false) {
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
-    setRandomNumberSeed(osrngseed());
+    setRandomNumberSeed(0);
     kineticEnergy = "m*v*v/2";
 }
 

openmmapi/src/LangevinIntegrator.cpp

@@ -33,7 +33,6 @@
 #include "openmm/Context.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
-#include "openmm/internal/OSRngSeed.h"
 #include "openmm/kernels.h"
 #include <string>
 
@@ -46,7 +45,7 @@ LangevinIntegrator::LangevinIntegrator(double temperature, double frictionCoeff,
     setFriction(frictionCoeff);
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
-    setRandomNumberSeed(osrngseed());
+    setRandomNumberSeed(0);
 }
 
 void LangevinIntegrator::initialize(ContextImpl& contextRef) {

openmmapi/src/MonteCarloAnisotropicBarostat.cpp

@@ -31,13 +31,12 @@
 
 #include "openmm/MonteCarloAnisotropicBarostat.h"
 #include "openmm/internal/MonteCarloAnisotropicBarostatImpl.h"
-#include "openmm/internal/OSRngSeed.h"
 
 using namespace OpenMM;
 
 MonteCarloAnisotropicBarostat::MonteCarloAnisotropicBarostat(const Vec3& defaultPressure, double temperature, bool scaleX, bool scaleY, bool scaleZ, int frequency) :
         defaultPressure(defaultPressure), temperature(temperature), scaleX(scaleX), scaleY(scaleY), scaleZ(scaleZ), frequency(frequency) {
-    setRandomNumberSeed(osrngseed());
+    setRandomNumberSeed(0);
 }
 
 ForceImpl* MonteCarloAnisotropicBarostat::createImpl() const {

openmmapi/src/MonteCarloAnisotropicBarostatImpl.cpp

@@ -31,6 +31,7 @@
 
 #include "openmm/internal/MonteCarloAnisotropicBarostatImpl.h"
 #include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/OSRngSeed.h"
 #include "openmm/Context.h"
 #include "openmm/kernels.h"
 #include <cmath>
@@ -60,7 +61,10 @@ void MonteCarloAnisotropicBarostatImpl::initialize(ContextImpl& context) {
         numAttempted[i] = 0;
         numAccepted[i] = 0;
     }
-    init_gen_rand(owner.getRandomNumberSeed(), random);
+    int randSeed = owner.getRandomNumberSeed();
+    // A random seed of 0 means use a unique one
+    if (randSeed == 0) randSeed = osrngseed();
+    init_gen_rand(randSeed, random);
 }
 
 void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context) {

openmmapi/src/MonteCarloBarostat.cpp

@@ -31,13 +31,12 @@
 
 #include "openmm/MonteCarloBarostat.h"
 #include "openmm/internal/MonteCarloBarostatImpl.h"
-#include "openmm/internal/OSRngSeed.h"
 
 using namespace OpenMM;
 
 MonteCarloBarostat::MonteCarloBarostat(double defaultPressure, double temperature, int frequency) :
         defaultPressure(defaultPressure), temperature(temperature), frequency(frequency) {
-    setRandomNumberSeed(osrngseed());
+    setRandomNumberSeed(0);
 }
 
 ForceImpl* MonteCarloBarostat::createImpl() const {

openmmapi/src/MonteCarloBarostatImpl.cpp

@@ -31,6 +31,7 @@
 
 #include "openmm/internal/MonteCarloBarostatImpl.h"
 #include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/OSRngSeed.h"
 #include "openmm/Context.h"
 #include "openmm/kernels.h"
 #include <cmath>
@@ -58,7 +59,10 @@ void MonteCarloBarostatImpl::initialize(ContextImpl& context) {
     volumeScale = 0.01*volume;
     numAttempted = 0;
     numAccepted = 0;
-    init_gen_rand(owner.getRandomNumberSeed(), random);
+    int randSeed = owner.getRandomNumberSeed();
+    // A random seed of 0 means use a unique one
+    if (randSeed == 0) randSeed = osrngseed();
+    init_gen_rand(randSeed, random);
 }
 
 void MonteCarloBarostatImpl::updateContextState(ContextImpl& context) {

openmmapi/src/MonteCarloMembraneBarostat.cpp

@@ -31,14 +31,13 @@
 
 #include "openmm/MonteCarloMembraneBarostat.h"
 #include "openmm/internal/MonteCarloMembraneBarostatImpl.h"
-#include "openmm/internal/OSRngSeed.h"
 
 using namespace OpenMM;
 
 MonteCarloMembraneBarostat::MonteCarloMembraneBarostat(double defaultPressure, double defaultSurfaceTension, double temperature, XYMode xymode, ZMode zmode, int frequency) :
         defaultPressure(defaultPressure), defaultSurfaceTension(defaultSurfaceTension), temperature(temperature),
         xymode(xymode), zmode(zmode), frequency(frequency) {
-    setRandomNumberSeed(osrngseed());
+    setRandomNumberSeed(0);
 }
 
 ForceImpl* MonteCarloMembraneBarostat::createImpl() const {

openmmapi/src/MonteCarloMembraneBarostatImpl.cpp

@@ -31,6 +31,7 @@
 
 #include "openmm/internal/MonteCarloMembraneBarostatImpl.h"
 #include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/OSRngSeed.h"
 #include "openmm/Context.h"
 #include "openmm/kernels.h"
 #include <cmath>
@@ -60,7 +61,10 @@ void MonteCarloMembraneBarostatImpl::initialize(ContextImpl& context) {
         numAttempted[i] = 0;
         numAccepted[i] = 0;
     }
-    init_gen_rand(owner.getRandomNumberSeed(), random);
+    int randSeed = owner.getRandomNumberSeed();
+    // A random seed of 0 means use a unique one
+    if (randSeed == 0) randSeed = osrngseed();
+    init_gen_rand(randSeed, random);
 }
 
 void MonteCarloMembraneBarostatImpl::updateContextState(ContextImpl& context) {