Implemented setting random seed for reference platform, as well as adding it...

Implemented setting random seed for reference platform, as well as adding it as an option to BrownianIntegrator and AndersenThermostat.

openmmapi/include/AndersenThermostat.h

@@ -76,10 +76,28 @@ public:
     double getDefaultCollisionFrequency() {
         return defaultFreq;
     }
+    /**
+     * Get the random number seed.  See setRandomNumberSeed() for details.
+     */
+    int getRandomNumberSeed() const {
+        return randomNumberSeed;
+    }
+    /**
+     * Set the random number seed.  The precise meaning of this parameter is undefined, and is left up
+     * to each Platform to interpret in an appropriate way.  It is guaranteed that if two simulations
+     * are run with different random number seeds, the sequence of collisions will be different.  On
+     * 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.
+     */
+    void setRandomNumberSeed(int seed) {
+        randomNumberSeed = seed;
+    }
 protected:
     ForceImpl* createImpl();
 private:
     double defaultTemp, defaultFreq;
+    int randomNumberSeed;
 };
 
 } // namespace OpenMM

openmmapi/include/BrownianIntegrator.h

@@ -78,6 +78,23 @@ public:
     void setFriction(double coeff) {
         friction = coeff;
     }
+    /**
+     * Get the random number seed.  See setRandomNumberSeed() for details.
+     */
+    int getRandomNumberSeed() const {
+        return randomNumberSeed;
+    }
+    /**
+     * Set the random number seed.  The precise meaning of this parameter is undefined, and is left up
+     * to each Platform to interpret in an appropriate way.  It is guaranteed that if two simulations
+     * are run with different random number seeds, the sequence of random forces will be different.  On
+     * 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.
+     */
+    void setRandomNumberSeed(int seed) {
+        randomNumberSeed = seed;
+    }
     /**
      * Advance a simulation through time by taking a series of time steps.
      * 
@@ -97,6 +114,7 @@ protected:
     std::vector<std::string> getKernelNames();
 private:
     double temperature, friction;
+    int randomNumberSeed;
     OpenMMContextImpl* context;
     Kernel kernel;
 };

openmmapi/include/LangevinIntegrator.h

@@ -79,16 +79,21 @@ public:
         friction = coeff;
     }
     /**
-     * Get the random number seed
+     * Get the random number seed.  See setRandomNumberSeed() for details.
      */
-    int getRandomNumberSeed(void ) const {
+    int getRandomNumberSeed() const {
         return randomNumberSeed;
     }
     /**
-     * Set the random number seed
+     * Set the random number seed.  The precise meaning of this parameter is undefined, and is left up
+     * to each Platform to interpret in an appropriate way.  It is guaranteed that if two simulations
+     * are run with different random number seeds, the sequence of random forces will be different.  On
+     * 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.
      */
-    void setRandomNumberSeed(int inputRandomNumberSeed) {
-        randomNumberSeed = inputRandomNumberSeed;
+    void setRandomNumberSeed(int seed) {
+        randomNumberSeed = seed;
     }
     /**
      * Advance a simulation through time by taking a series of time steps.

openmmapi/src/AndersenThermostat.cpp

@@ -31,11 +31,13 @@
 
 #include "AndersenThermostat.h"
 #include "internal/AndersenThermostatImpl.h"
+#include <ctime>
 
 using namespace OpenMM;
 
 AndersenThermostat::AndersenThermostat(double defaultTemperature, double defaultCollisionFrequency) :
-    defaultTemp(defaultTemperature), defaultFreq(defaultCollisionFrequency) {
+        defaultTemp(defaultTemperature), defaultFreq(defaultCollisionFrequency) {
+    setRandomNumberSeed((int) time(NULL));
 }
 
 ForceImpl* AndersenThermostat::createImpl() {

openmmapi/src/BrownianIntegrator.cpp

@@ -33,6 +33,7 @@
 #include "OpenMMContext.h"
 #include "internal/OpenMMContextImpl.h"
 #include "kernels.h"
+#include <ctime>
 #include <string>
 
 using namespace OpenMM;
@@ -44,6 +45,7 @@ BrownianIntegrator::BrownianIntegrator(double temperature, double frictionCoeff,
     setFriction(frictionCoeff);
     setStepSize(stepSize);
     setConstraintTolerance(1e-4);
+    setRandomNumberSeed((int) time(NULL));
 }
 
 void BrownianIntegrator::initialize(OpenMMContextImpl& contextRef) {

openmmapi/src/LangevinIntegrator.cpp

@@ -33,6 +33,7 @@
 #include "OpenMMContext.h"
 #include "internal/OpenMMContextImpl.h"
 #include "kernels.h"
+#include <ctime>
 #include <string>
 
 using namespace OpenMM;
@@ -44,7 +45,7 @@ LangevinIntegrator::LangevinIntegrator(double temperature, double frictionCoeff,
     setFriction(frictionCoeff);
     setStepSize(stepSize);
     setConstraintTolerance(1e-4);
-    setRandomNumberSeed(1);
+    setRandomNumberSeed((int) time(NULL));
 }
 
 void LangevinIntegrator::initialize(OpenMMContextImpl& contextRef) {

platforms/cuda/src/CudaKernels.cpp

@@ -432,15 +432,9 @@ CudaIntegrateLangevinStepKernel::~CudaIntegrateLangevinStepKernel() {
 
 void CudaIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
     initializeIntegration(system, data, integrator);
-
-    // if LangevinIntegrator seed does not equal default value or is less than/equal to 0, then
-    // set gpu seed and redo random values
-
-    if( integrator.getRandomNumberSeed() <= 1 ){
-       _gpuContext* gpu = data.gpu;
-       gpu->seed        = static_cast<unsigned long>( integrator.getRandomNumberSeed() );
-       gpuInitializeRandoms( gpu );
-    }
+    _gpuContext* gpu = data.gpu;
+    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
+    gpuInitializeRandoms(gpu);
     prevStepSize = -1.0;
 }
 
@@ -478,6 +472,9 @@ CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
 
 void CudaIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
     initializeIntegration(system, data, integrator);
+    _gpuContext* gpu = data.gpu;
+    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
+    gpuInitializeRandoms(gpu);
     prevStepSize = -1.0;
 }
 
@@ -512,6 +509,9 @@ CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
 }
 
 void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
+    _gpuContext* gpu = data.gpu;
+    gpu->seed = (unsigned long) thermostat.getRandomNumberSeed();
+    gpuInitializeRandoms(gpu);
     prevStepSize = -1.0;
 }
 

platforms/cuda/src/kernels/gpu.cpp

@@ -38,7 +38,6 @@
 #include <iostream>
 #include <fstream>
 #include <sstream>
-#include <ctime>
 #include <cmath>
 #include <map>
 #include <algorithm>
@@ -1199,7 +1198,7 @@ int gpuAllocateInitialBuffers(gpuContext gpu)
         gpu->psAtomIndex->_pSysStream[0][i] = i;
     gpu->psAtomIndex->Upload();
     // Determine randoms
-    gpu->seed                           = (unsigned long)time(NULL) & 0x000fffff;
+    gpu->seed                           = 1;
     gpu->sim.randomFrames               = 95;
     gpu->sim.randomIterations           = gpu->sim.randomFrames;
     gpu->sim.randoms                    = gpu->sim.randomFrames * gpu->sim.paddedNumberOfAtoms - 5 * GRID;
@@ -1215,35 +1214,6 @@ int gpuAllocateInitialBuffers(gpuContext gpu)
     gpu->sim.pRandom2b                  = gpu->psRandom2->_pDevStream[0] + gpu->sim.totalRandoms;
     gpu->sim.pRandomPosition            = gpu->psRandomPosition->_pDevStream[0];
     gpu->sim.pRandomSeed                = gpu->psRandomSeed->_pDevStream[0];
-    for (int i = 0; i < (int) gpu->sim.blocks; i++)
-    {
-        gpu->psRandomPosition->_pSysStream[0][i] = 0;
-    }
-    int seed = gpu->seed | ((gpu->seed ^ 0xffffffff) << 16);
-    srand(seed);
-    for (int i = 0; i < (int) (gpu->sim.blocks * gpu->sim.random_threads_per_block); i++)
-    {
-        gpu->psRandomSeed->_pSysStream[0][i].x = rand();
-        gpu->psRandomSeed->_pSysStream[0][i].y = rand();
-        gpu->psRandomSeed->_pSysStream[0][i].z = rand();
-        gpu->psRandomSeed->_pSysStream[0][i].w = rand();
-    }
-
-    float randomValue = 0.0f;
-    for (int i = 0; i < (int) gpu->sim.totalRandomsTimesTwo; i++)
-    {
-        gpu->psRandom4->_pSysStream[0][i].x         = randomValue;
-        gpu->psRandom4->_pSysStream[0][i].y         = randomValue;
-        gpu->psRandom4->_pSysStream[0][i].z         = randomValue;
-        gpu->psRandom4->_pSysStream[0][i].w         = randomValue;
-        gpu->psRandom2->_pSysStream[0][i].x         = randomValue;
-        gpu->psRandom2->_pSysStream[0][i].y         = randomValue;
-    }
-
-    gpu->psRandomSeed->Upload();
-    gpu->psRandom4->Upload();
-    gpu->psRandom2->Upload();
-    gpu->psRandomPosition->Upload();
 
     // Allocate and clear linear momentum buffer
     gpu->psLinearMomentum = new CUDAStream<float4>(gpu->sim.blocks, 1);

platforms/cuda/src/kernels/kBrownianUpdate.cu

@@ -141,12 +141,11 @@ void kBrownianUpdatePart2(gpuContext gpu)
     LAUNCHERROR("kBrownianUpdatePart2");
     
     // Update randoms if necessary
-    static int iteration = 0;
-    iteration++;
-    if (iteration == gpu->sim.randomIterations)
+    gpu->iterations++;
+    if (gpu->iterations == gpu->sim.randomIterations)
     {
         kGenerateRandoms(gpu);
-        iteration = 0;
+        gpu->iterations = 0;
     }
 }
 

platforms/cuda/src/kernels/kCalculateAndersenThermostat.cu

@@ -94,12 +94,11 @@ void kCalculateAndersenThermostat(gpuContext gpu)
     LAUNCHERROR("kCalculateAndersenThermostat");
     
     // Update randoms if necessary
-    static int iteration = 0;
-    iteration++;
-    if (iteration == gpu->sim.randomIterations)
+    gpu->iterations++;
+    if (gpu->iterations == gpu->sim.randomIterations)
     {
         kGenerateRandoms(gpu);
-        iteration = 0;
+        gpu->iterations = 0;
     }
 }
 

platforms/cuda/src/kernels/kUpdateShakeH.cu

@@ -720,12 +720,11 @@ void kUpdatePart2(gpuContext gpu)
     }
     
     // Update randoms if necessary
-    static int iteration = 0;
-    iteration++;
-    if (iteration == gpu->sim.randomIterations)
+    gpu->iterations++;
+    if (gpu->iterations == gpu->sim.randomIterations)
     {
         kGenerateRandoms(gpu);
-        iteration = 0;
+        gpu->iterations = 0;
     }
 }
 

platforms/cuda/tests/TestCudaAndersenThermostat.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -86,9 +86,71 @@ void testTemperature() {
     ASSERT_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
 }
 
+void testRandomSeed() {
+    const int numParticles = 8;
+    const double temp = 100.0;
+    const double collisionFreq = 10.0;
+    CudaPlatform platform;
+    System system(numParticles, 0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* forceField = new NonbondedForce(numParticles, 0);
+    for (int i = 0; i < numParticles; ++i) {
+        system.setParticleMass(i, 2.0);
+        forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
+    }
+    system.addForce(forceField);
+    AndersenThermostat* thermostat = new AndersenThermostat(temp, collisionFreq);
+    system.addForce(thermostat);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
+        velocities[i] = Vec3(0, 0, 0);
+    }
+
+    // Try twice with the same random seed.
+
+    thermostat->setRandomNumberSeed(5);
+    OpenMMContext context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state1 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state2 = context.getState(State::Positions);
+
+    // Try twice with a different random seed.
+
+    thermostat->setRandomNumberSeed(10);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state3 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state4 = context.getState(State::Positions);
+
+    // Compare the results.
+
+    for (int i = 0; i < numParticles; i++) {
+        for (int j = 0; j < 3; j++) {
+            ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
+            ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
+            ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
+        }
+    }
+}
+
 int main() {
     try {
         testTemperature();
+        testRandomSeed();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/cuda/tests/TestCudaBrownianIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -168,11 +168,71 @@ void testConstraints() {
     }
 }
 
+void testRandomSeed() {
+    const int numParticles = 8;
+    const double temp = 100.0;
+    const double collisionFreq = 10.0;
+    CudaPlatform platform;
+    System system(numParticles, 0);
+    BrownianIntegrator integrator(temp, 2.0, 0.001);
+    NonbondedForce* forceField = new NonbondedForce(numParticles, 0);
+    for (int i = 0; i < numParticles; ++i) {
+        system.setParticleMass(i, 2.0);
+        forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
+    }
+    system.addForce(forceField);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
+        velocities[i] = Vec3(0, 0, 0);
+    }
+
+    // Try twice with the same random seed.
+
+    integrator.setRandomNumberSeed(5);
+    OpenMMContext context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state1 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state2 = context.getState(State::Positions);
+
+    // Try twice with a different random seed.
+
+    integrator.setRandomNumberSeed(10);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state3 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state4 = context.getState(State::Positions);
+
+    // Compare the results.
+
+    for (int i = 0; i < numParticles; i++) {
+        for (int j = 0; j < 3; j++) {
+            ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
+            ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
+            ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
+        }
+    }
+}
+
 int main() {
     try {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testRandomSeed();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/cuda/tests/TestCudaLangevinIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -173,11 +173,71 @@ void testConstraints() {
     }
 }
 
+void testRandomSeed() {
+    const int numParticles = 8;
+    const double temp = 100.0;
+    const double collisionFreq = 10.0;
+    CudaPlatform platform;
+    System system(numParticles, 0);
+    LangevinIntegrator integrator(temp, 2.0, 0.01);
+    NonbondedForce* forceField = new NonbondedForce(numParticles, 0);
+    for (int i = 0; i < numParticles; ++i) {
+        system.setParticleMass(i, 2.0);
+        forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
+    }
+    system.addForce(forceField);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
+        velocities[i] = Vec3(0, 0, 0);
+    }
+
+    // Try twice with the same random seed.
+
+    integrator.setRandomNumberSeed(5);
+    OpenMMContext context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state1 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state2 = context.getState(State::Positions);
+
+    // Try twice with a different random seed.
+
+    integrator.setRandomNumberSeed(10);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state3 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state4 = context.getState(State::Positions);
+
+    // Compare the results.
+
+    for (int i = 0; i < numParticles; i++) {
+        for (int j = 0; j < 3; j++) {
+            ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
+            ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
+            ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
+        }
+    }
+}
+
 int main() {
     try {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testRandomSeed();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/reference/src/ReferenceKernels.cpp

@@ -48,6 +48,7 @@
 #include "System.h"
 #include "internal/OpenMMContextImpl.h"
 #include "Integrator.h"
+#include "SimTKUtilities/SimTKOpenMMUtilities.h"
 #include <cmath>
 #include <limits>
 
@@ -538,6 +539,7 @@ void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, cons
         constraintIndices[i][1] = particle2;
         shakeParameters[i][0] = static_cast<RealOpenMM>(distance);
     }
+    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
 }
 
 void ReferenceIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, const LangevinIntegrator& integrator) {
@@ -599,6 +601,7 @@ void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, cons
         constraintIndices[i][1] = particle2;
         shakeParameters[i][0] = static_cast<RealOpenMM>(distance);
     }
+    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
 }
 
 void ReferenceIntegrateBrownianStepKernel::execute(OpenMMContextImpl& context, const BrownianIntegrator& integrator) {
@@ -643,6 +646,7 @@ void ReferenceApplyAndersenThermostatKernel::initialize(const System& system, co
     for (int i = 0; i < numParticles; ++i)
         masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
     this->thermostat = new ReferenceAndersenThermostat();
+    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
 }
 
 void ReferenceApplyAndersenThermostatKernel::execute(OpenMMContextImpl& context) {

platforms/reference/src/SimTKUtilities/SimTKOpenMMUtilities.cpp

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -34,6 +34,8 @@
 
 uint32_t SimTKOpenMMUtilities::_randomNumberSeed = 0;
 bool SimTKOpenMMUtilities::_randomInitialized = false;
+bool SimTKOpenMMUtilities::nextGaussianIsValid = false;
+RealOpenMM SimTKOpenMMUtilities::nextGaussian = 0;
 
 /* ---------------------------------------------------------------------------------------
 
@@ -1428,15 +1430,14 @@ int SimTKOpenMMUtilities::addTwoDimArray( int dimension1, int dimension2,
    --------------------------------------------------------------------------------------- */
 
 RealOpenMM SimTKOpenMMUtilities::getNormallyDistributedRandomNumber( void ) {
-    static bool nextValueIsValid = false;
-    static RealOpenMM nextValue = 0;
-    if (nextValueIsValid) {
-        nextValueIsValid = false;
-        return nextValue;
+    if (nextGaussianIsValid) {
+        nextGaussianIsValid = false;
+        return nextGaussian;
     }
     if (!_randomInitialized) {
         init_gen_rand(_randomNumberSeed);
         _randomInitialized = true;
+        nextGaussianIsValid = false;
     }
     
     // Use the polar form of the Box-Muller transformation to generate two Gaussian random numbers.
@@ -1448,8 +1449,8 @@ RealOpenMM SimTKOpenMMUtilities::getNormallyDistributedRandomNumber( void ) {
         r2 = x*x + y*y;
     } while (r2 >= 1.0 || r2 == 0.0);
     RealOpenMM multiplier = static_cast<RealOpenMM>( sqrt((-2.0*log(r2))/r2) );
-    nextValue = y*multiplier;
-    nextValueIsValid = true;
+    nextGaussian = y*multiplier;
+    nextGaussianIsValid = true;
     return x*multiplier;
 }
 
@@ -1465,8 +1466,10 @@ RealOpenMM SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber( void ) {
     if (!_randomInitialized) {
         init_gen_rand(_randomNumberSeed);
         _randomInitialized = true;
+        nextGaussianIsValid = false;
     }
-    return static_cast<RealOpenMM>( genrand_real2() );
+    RealOpenMM value = static_cast<RealOpenMM>( genrand_real2() );
+    return value;
 }
 
 /**---------------------------------------------------------------------------------------
@@ -1507,5 +1510,5 @@ void SimTKOpenMMUtilities::setRandomNumberSeed( uint32_t seed ) {
    // ---------------------------------------------------------------------------------------
 
    _randomNumberSeed = seed;
-   init_gen_rand(_randomNumberSeed);
+   _randomInitialized = false;
 }

platforms/reference/src/SimTKUtilities/SimTKOpenMMUtilities.h

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -58,6 +58,8 @@ class OPENMM_EXPORT SimTKOpenMMUtilities {
 
        static uint32_t _randomNumberSeed;
        static bool _randomInitialized;
+       static bool nextGaussianIsValid;
+       static RealOpenMM nextGaussian;
 
    public:
 

platforms/reference/tests/TestReferenceAndersenThermostat.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -86,9 +86,71 @@ void testTemperature() {
     ASSERT_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
 }
 
+void testRandomSeed() {
+    const int numParticles = 8;
+    const double temp = 100.0;
+    const double collisionFreq = 10.0;
+    ReferencePlatform platform;
+    System system(numParticles, 0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* forceField = new NonbondedForce(numParticles, 0);
+    for (int i = 0; i < numParticles; ++i) {
+        system.setParticleMass(i, 2.0);
+        forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
+    }
+    system.addForce(forceField);
+    AndersenThermostat* thermostat = new AndersenThermostat(temp, collisionFreq);
+    system.addForce(thermostat);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
+        velocities[i] = Vec3(0, 0, 0);
+    }
+
+    // Try twice with the same random seed.
+
+    thermostat->setRandomNumberSeed(5);
+    OpenMMContext context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state1 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state2 = context.getState(State::Positions);
+
+    // Try twice with a different random seed.
+
+    thermostat->setRandomNumberSeed(10);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state3 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state4 = context.getState(State::Positions);
+
+    // Compare the results.
+
+    for (int i = 0; i < numParticles; i++) {
+        for (int j = 0; j < 3; j++) {
+            ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
+            ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
+            ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
+        }
+    }
+}
+
 int main() {
     try {
         testTemperature();
+        testRandomSeed();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/reference/tests/TestReferenceBrownianIntegrator.cpp

@@ -162,11 +162,71 @@ void testConstraints() {
     }
 }
 
+void testRandomSeed() {
+    const int numParticles = 8;
+    const double temp = 100.0;
+    const double collisionFreq = 10.0;
+    ReferencePlatform platform;
+    System system(numParticles, 0);
+    BrownianIntegrator integrator(temp, 2.0, 0.001);
+    NonbondedForce* forceField = new NonbondedForce(numParticles, 0);
+    for (int i = 0; i < numParticles; ++i) {
+        system.setParticleMass(i, 2.0);
+        forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
+    }
+    system.addForce(forceField);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
+        velocities[i] = Vec3(0, 0, 0);
+    }
+
+    // Try twice with the same random seed.
+
+    integrator.setRandomNumberSeed(5);
+    OpenMMContext context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state1 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state2 = context.getState(State::Positions);
+
+    // Try twice with a different random seed.
+
+    integrator.setRandomNumberSeed(10);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state3 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state4 = context.getState(State::Positions);
+
+    // Compare the results.
+
+    for (int i = 0; i < numParticles; i++) {
+        for (int j = 0; j < 3; j++) {
+            ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
+            ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
+            ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
+        }
+    }
+}
+
 int main() {
     try {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testRandomSeed();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/reference/tests/TestReferenceLangevinIntegrator.cpp

@@ -169,11 +169,71 @@ void testConstraints() {
     }
 }
 
+void testRandomSeed() {
+    const int numParticles = 8;
+    const double temp = 100.0;
+    const double collisionFreq = 10.0;
+    ReferencePlatform platform;
+    System system(numParticles, 0);
+    LangevinIntegrator integrator(temp, 2.0, 0.01);
+    NonbondedForce* forceField = new NonbondedForce(numParticles, 0);
+    for (int i = 0; i < numParticles; ++i) {
+        system.setParticleMass(i, 2.0);
+        forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
+    }
+    system.addForce(forceField);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
+        velocities[i] = Vec3(0, 0, 0);
+    }
+
+    // Try twice with the same random seed.
+
+    integrator.setRandomNumberSeed(5);
+    OpenMMContext context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state1 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state2 = context.getState(State::Positions);
+
+    // Try twice with a different random seed.
+
+    integrator.setRandomNumberSeed(10);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state3 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state4 = context.getState(State::Positions);
+
+    // Compare the results.
+
+    for (int i = 0; i < numParticles; i++) {
+        for (int j = 0; j < 3; j++) {
+            ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
+            ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
+            ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
+        }
+    }
+}
+
 int main() {
     try {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testRandomSeed();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;