Created reference implementation of AndersenThermostat

openmmapi/include/VerletIntegrator.h

@@ -38,7 +38,7 @@
 namespace OpenMM {
 
 /**
- * This is an Integrator which simulates a System using the velocity Verlet algorithm.
+ * This is an Integrator which simulates a System using the leap-frog Verlet algorithm.
  */
 
 class VerletIntegrator : public Integrator {

openmmapi/src/AndersenThermostatImpl.cpp

@@ -32,12 +32,20 @@
 #include "internal/AndersenThermostatImpl.h"
 #include "internal/OpenMMContextImpl.h"
 #include "Integrator.h"
+#include "System.h"
 #include "kernels.h"
+#include <vector>
 
 using namespace OpenMM;
+using std::vector;
 
 AndersenThermostatImpl::AndersenThermostatImpl(AndersenThermostat& owner, OpenMMContextImpl& context) : owner(owner) {
     kernel = context.getPlatform().createKernel(ApplyAndersenThermostatKernel::Name());
+    const System& system = context.getSystem();
+    vector<double> masses(system.getNumAtoms());
+    for (int i = 0; i < system.getNumAtoms(); ++i)
+        masses[i] = system.getAtomMass(i);
+    dynamic_cast<ApplyAndersenThermostatKernel&>(kernel.getImpl()).initialize(masses);
 }
 
 void AndersenThermostatImpl::updateContextState(OpenMMContextImpl& context) {

openmmapi/src/OpenMMContextImpl.cpp

@@ -36,9 +36,11 @@
 #include "kernels.h"
 #include "internal/ForceImpl.h"
 #include "internal/OpenMMContextImpl.h"
+#include <map>
 #include <vector>
 
 using namespace OpenMM;
+using std::map;
 using std::vector;
 using std::string;
 
@@ -48,6 +50,8 @@ OpenMMContextImpl::OpenMMContextImpl(OpenMMContext& owner, System& system, Integ
     kernelNames.push_back(CalcKineticEnergyKernel::Name());
     for (int i = 0; i < system.getNumForces(); ++i) {
         forceImpls.push_back(system.getForce(i).createImpl(*this));
+        map<string, double> forceParameters = forceImpls[forceImpls.size()-1]->getDefaultParameters();
+        parameters.insert(forceParameters.begin(), forceParameters.end());
         vector<string> forceKernels = forceImpls[forceImpls.size()-1]->getKernelNames();
         kernelNames.insert(kernelNames.begin(), forceKernels.begin(), forceKernels.end());
     }

platforms/reference/src/ReferenceKernels.cpp

@@ -32,6 +32,7 @@
 #include "ReferenceKernels.h"
 #include "ReferenceFloatStreamImpl.h"
 #include "gbsa/CpuObc.h"
+#include "SimTKReference/ReferenceAndersenThermostat.h"
 #include "SimTKReference/ReferenceAngleBondIxn.h"
 #include "SimTKReference/ReferenceBondForce.h"
 #include "SimTKReference/ReferenceBrownianDynamics.h"
@@ -401,12 +402,23 @@ void ReferenceIntegrateBrownianStepKernel::execute(Stream& positions, Stream& ve
     dynamics->update(positions.getSize(), posData, velData, forceData, masses);
 }
 
+ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
+    if (thermostat)
+        delete thermostat;
+    if (masses)
+        delete[] masses;
+}
+
 void ReferenceApplyAndersenThermostatKernel::initialize(const vector<double>& masses) {
-    
+    this->masses = new RealOpenMM[masses.size()];
+    for (int i = 0; i < masses.size(); ++i)
+        this->masses[i] = masses[i];
+    thermostat = new ReferenceAndersenThermostat();
 }
 
 void ReferenceApplyAndersenThermostatKernel::execute(Stream& velocities, double temperature, double collisionFrequency, double stepSize) {
-    
+    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) velocities.getImpl()).getData();
+    thermostat->applyThermostat(velocities.getSize(), velData, masses, temperature, collisionFrequency, stepSize);
 }
 
 void ReferenceCalcKineticEnergyKernel::initialize(const vector<double>& masses) {

platforms/reference/src/ReferenceKernels.h

@@ -36,6 +36,7 @@
 #include "SimTKUtilities/SimTKOpenMMRealType.h"
 
 class CpuObc;
+class ReferenceAndersenThermostat;
 class ReferenceBrownianDynamics;
 class ReferenceStochasticDynamics;
 class ReferenceShakeAlgorithm;
@@ -254,8 +255,9 @@ private:
  */
 class ReferenceApplyAndersenThermostatKernel : public ApplyAndersenThermostatKernel {
 public:
-    ReferenceApplyAndersenThermostatKernel(std::string name, const Platform& platform) : ApplyAndersenThermostatKernel(name, platform) {
+    ReferenceApplyAndersenThermostatKernel(std::string name, const Platform& platform) : ApplyAndersenThermostatKernel(name, platform), thermostat(0) {
     }
+    ~ReferenceApplyAndersenThermostatKernel();
     /**
      * Initialize the kernel, setting up the values of unchanging parameters.
      * 
@@ -271,6 +273,9 @@ public:
      * @param stepSize           the integration step size
      */
     void execute(Stream& velocities, double temperature, double collisionFrequency, double stepSize);
+private:
+    ReferenceAndersenThermostat* thermostat;
+    RealOpenMM* masses;
 };
 
 /**

platforms/reference/src/SimTKReference/ReferenceAndersenThermostat.cpp

+
+/* Portions copyright (c) 2008 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include <string.h>
+#include <sstream>
+
+#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
+#include "ReferenceAndersenThermostat.h"
+
+      /**---------------------------------------------------------------------------------------
+      
+         Constructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceAndersenThermostat::ReferenceAndersenThermostat( ) {
+       }
+
+      /**---------------------------------------------------------------------------------------
+      
+         Destructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceAndersenThermostat::~ReferenceAndersenThermostat( ) {
+       }
+
+      /**---------------------------------------------------------------------------------------
+      
+         Apply the thermostat at the start of a time step.
+      
+         @param numberOfAtoms      number of atoms
+         @param atomVelocities     atom velocities
+         @param temperature        thermostat temperature in Kelvin
+         @param collisionFrequency collision frequency for each atom in fs^-1
+         @param stepSize           integration step size in fs
+                  
+         --------------------------------------------------------------------------------------- */
+          
+      void ReferenceAndersenThermostat::applyThermostat( int numberOfAtoms, RealOpenMM** atomVelocities, RealOpenMM* atomMasses,
+              RealOpenMM temperature, RealOpenMM collisionFrequency, RealOpenMM stepSize ) const {
+          
+          const RealOpenMM collisionProbability = collisionFrequency*stepSize;
+          for (int i = 0; i < numberOfAtoms; ++i) {
+              if (SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber() < collisionProbability) {
+                  
+                  // A collision occurred, so set the velocity to a new value chosen from a Boltzmann distribution.
+                  
+                  const RealOpenMM velocityScale = sqrt(BOLTZ*temperature/atomMasses[i]);
+                  atomVelocities[i][0] = velocityScale*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+                  atomVelocities[i][1] = velocityScale*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+                  atomVelocities[i][2] = velocityScale*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+              }
+          }
+          
+      }
+

platforms/reference/src/SimTKReference/ReferenceAndersenThermostat.h

+
+/* Portions copyright (c) 2008 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef __ReferenceAndersenThermostat_H__
+#define __ReferenceAndersenThermostat_H__
+
+#include "../SimTKUtilities/SimTKOpenMMCommon.h"
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceAndersenThermostat {
+
+   private:
+       
+   public:
+
+      /**---------------------------------------------------------------------------------------
+      
+         Constructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceAndersenThermostat( );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Destructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceAndersenThermostat( );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Apply the thermostat at the start of a time step.
+      
+         @param numberOfAtoms      number of atoms
+         @param atomVelocities     atom velocities
+         @param atomMasses         atom masses
+         @param temperature        thermostat temperature in Kelvin
+         @param collisionFrequency collision frequency for each atom in fs^-1
+         @param stepSize           integration step size in fs
+                  
+         --------------------------------------------------------------------------------------- */
+          
+      void applyThermostat( int numberOfAtoms, RealOpenMM** atomVelocities, RealOpenMM* atomMasses,
+              RealOpenMM temperature, RealOpenMM collisionFrequency, RealOpenMM stepSize ) const;
+      
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // __ReferenceAndersenThermostat_H__

platforms/reference/src/SimTKReference/ReferenceBrownianDynamics.cpp

@@ -200,7 +200,7 @@ int ReferenceBrownianDynamics::update( int numberOfAtoms, RealOpenMM** atomCoord
    const RealOpenMM forceScale = getDeltaT()/getFriction();
    for (int i = 0; i < numberOfAtoms; ++i) {
        for (int j = 0; j < 3; ++j) {
-           xPrime[i][j] = atomCoordinates[i][j] + forceScale*forces[i][j] + noiseAmplitude*getNormallyDistributedRandomNumber();
+           xPrime[i][j] = atomCoordinates[i][j] + forceScale*forces[i][j] + noiseAmplitude*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
        }
    }
    ReferenceShakeAlgorithm* referenceShakeAlgorithm = getReferenceShakeAlgorithm();

platforms/reference/src/SimTKReference/ReferenceDynamics.cpp

@@ -28,7 +28,6 @@
 #include "../SimTKUtilities/SimTKOpenMMCommon.h"
 #include "../SimTKUtilities/SimTKOpenMMLog.h"
 #include "../SimTKUtilities/SimTKOpenMMUtilities.h"
-#include "../sfmt/SFMT.h"
 #include "ReferenceDynamics.h"
 
 const int ReferenceDynamics::DefaultReturn      = 0;
@@ -57,8 +56,6 @@ ReferenceDynamics::ReferenceDynamics( int numberOfAtoms,  RealOpenMM deltaT, Rea
    // ---------------------------------------------------------------------------------------
 
    _timeStep             = 0;
-   _randomNumberSeed     = 0;
-   init_gen_rand(_randomNumberSeed);
 
    _twoDTempArrays       = 0;
    _twoDTempArrays       = NULL;
@@ -384,77 +381,6 @@ RealOpenMM ReferenceDynamics::getTemperature( void ) const {
    return _temperature;
 }
 
-/**---------------------------------------------------------------------------------------
-
-   Get normally distributed random number
-
-   @return random value
-
-   --------------------------------------------------------------------------------------- */
-
-RealOpenMM ReferenceDynamics::getNormallyDistributedRandomNumber( void ){
-    static bool nextValueIsValid = false;
-    static RealOpenMM nextValue = 0;
-    if (nextValueIsValid) {
-        nextValueIsValid = false;
-        return nextValue;
-    }
-    
-    // Use the polar form of the Box-Muller transformation to generate two Gaussian random numbers.
-    
-    RealOpenMM x, y, r2;
-    do {
-        x = 2.0*genrand_real2()-1.0;
-        y = 2.0*genrand_real2()-1.0;
-        r2 = x*x + y*y;
-    } while (r2 >= 1.0 || r2 == 0.0);
-    RealOpenMM multiplier = sqrt((-2.0*log(r2))/r2);
-    nextValue = y*multiplier;
-    nextValueIsValid = true;
-    return x*multiplier;
-}
-
-/**---------------------------------------------------------------------------------------
-
-   Get random number seed
-
-   @return random number seed
-
-   --------------------------------------------------------------------------------------- */
-
-uint32_t ReferenceDynamics::getRandomNumberSeed( void ) const {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName  = "\nReferenceDynamics::getRandomNumberSeed";
-
-   // ---------------------------------------------------------------------------------------
-
-   return _randomNumberSeed;
-}
-
-/**---------------------------------------------------------------------------------------
-
-   Set random number seed
-
-   @param seed    new seed value
-
-   @return ReferenceDynamics::DefaultReturn
-
-   --------------------------------------------------------------------------------------- */
-
-int ReferenceDynamics::setRandomNumberSeed( uint32_t seed ){
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName  = "\nReferenceDynamics::setRandomNumberSeed";
-
-   // ---------------------------------------------------------------------------------------
-
-   _randomNumberSeed = seed;
-   return ReferenceDynamics::DefaultReturn;
-}
-
 /**---------------------------------------------------------------------------------------
 
    Get ReferenceShake
@@ -664,7 +590,7 @@ int ReferenceDynamics::printParameters( std::stringstream& message ) const {
    message << "delta_t=" << getDeltaT()          << " ";
    message << "temperature=" << getTemperature() << " ";
    message << "step=" << getTimeStep()           << " ";
-   message << "seed=" << getRandomNumberSeed()   << " ";
+   message << "seed=" << SimTKOpenMMUtilities::getRandomNumberSeed()   << " ";
    message << std::endl;
 
    return ReferenceDynamics::DefaultReturn;
@@ -759,7 +685,7 @@ int ReferenceDynamics::writeState( int numberOfAtoms, RealOpenMM** atomCoordinat
    scalarR.push_back( getDeltaT() );
    scalarNameR.push_back( "delta_t" );
 
-   scalarR.push_back( getRandomNumberSeed() );
+   scalarR.push_back( SimTKOpenMMUtilities::getRandomNumberSeed() );
    scalarNameR.push_back( "seed" );
 
    scalarR.push_back( getTemperature() );

platforms/reference/src/SimTKReference/ReferenceDynamics.h

@@ -61,7 +61,6 @@ class ReferenceDynamics {
 
       RealOpenMM _deltaT;
       RealOpenMM _temperature;
-      uint32_t _randomNumberSeed;
 
       int _numberOf2DTempArrays;
       RealOpenMM*** _twoDTempArrays;
@@ -259,38 +258,6 @@ class ReferenceDynamics {
       
       virtual int update( int numberOfAtoms, RealOpenMM** atomCoordinates,
                           RealOpenMM** velocities, RealOpenMM** forces, RealOpenMM* masses );
-      
-      /**---------------------------------------------------------------------------------------
-      
-         Get normally distributed random number
-      
-         @return random value
-      
-         --------------------------------------------------------------------------------------- */
-      
-      RealOpenMM getNormallyDistributedRandomNumber( void );
-      
-      /**---------------------------------------------------------------------------------------
-      
-         Get random number seed
-      
-         @return random number seed
-      
-         --------------------------------------------------------------------------------------- */
-      
-      uint32_t getRandomNumberSeed( void ) const;
-      
-      /**---------------------------------------------------------------------------------------
-      
-         Set random number seed
-      
-         @param seed    new seed value
-      
-         @return DefaultReturn
-      
-         --------------------------------------------------------------------------------------- */
-      
-      int setRandomNumberSeed( uint32_t seed );
 
       /**---------------------------------------------------------------------------------------
       

platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.cpp

@@ -284,9 +284,9 @@ int ReferenceStochasticDynamics::updatePart1( int numberOfAtoms, RealOpenMM** at
          oldVelocities[ii][jj] = velocities[ii][jj];
 
          RealOpenMM Vmh        = xVector[ii][jj]*fixedParameters[D]/(tau*fixedParameters[C]) +
-                                 sqrtInvMass*fixedParameters[Yv]*getNormallyDistributedRandomNumber();
+                                 sqrtInvMass*fixedParameters[Yv]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
 
-           vVector[ii][jj]     = sqrtInvMass*fixedParameters[V]*getNormallyDistributedRandomNumber();
+           vVector[ii][jj]     = sqrtInvMass*fixedParameters[V]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
           
            velocities[ii][jj]  = oldVelocities[ii][jj]*fixedParameters[EM] +
                                  inverseMasses[ii]*forces[ii][jj]*tau*( one - fixedParameters[EM]) +
@@ -375,9 +375,9 @@ int ReferenceStochasticDynamics::updatePart2( int numberOfAtoms, RealOpenMM** at
 
            velocities[ii][jj]  = (xPrime[ii][jj] - atomCoordinates[ii][jj])*fix1;
          RealOpenMM Xmh        = vVector[ii][jj]*tau*fixedParameters[D]/(fixedParameters[EM] - one) +
-                                 sqrtInvMass*fixedParameters[Yx]*getNormallyDistributedRandomNumber();
+                                 sqrtInvMass*fixedParameters[Yx]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
 
-           xVector[ii][jj]     = sqrtInvMass*fixedParameters[X]*getNormallyDistributedRandomNumber();
+           xVector[ii][jj]     = sqrtInvMass*fixedParameters[X]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
             
            xPrime[ii][jj]     += xVector[ii][jj] - Xmh;
       }
@@ -480,7 +480,7 @@ int ReferenceStochasticDynamics::update( int numberOfAtoms, RealOpenMM** atomCoo
       for( int ii = 0; ii < numberOfAtoms; ii++ ){
          RealOpenMM sqrtInverseMass = SQRT( inverseMasses[ii] )*fixedParameters[X];
          for( int jj = 0; jj < 3; jj++ ){
-            xVector[ii][jj] = sqrtInverseMass*getNormallyDistributedRandomNumber();
+            xVector[ii][jj] = sqrtInverseMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
          }
       }
 

platforms/reference/src/SimTKUtilities/SimTKOpenMMUtilities.cpp

@@ -26,11 +26,15 @@
 
 #include "SimTKOpenMMUtilities.h"
 #include "SimTKOpenMMLog.h"
+#include "../sfmt/SFMT.h"
 
 // fabs(), ...
 
 #include <math.h>
 
+uint32_t SimTKOpenMMUtilities::_randomNumberSeed = 0;
+bool SimTKOpenMMUtilities::_randomInitialized = false;
+
 /* ---------------------------------------------------------------------------------------
 
    Find distances**2 from a given atom (Simbios)
@@ -1414,3 +1418,94 @@ int SimTKOpenMMUtilities::addTwoDimArray( int dimension1, int dimension2,
       
    return  SimTKOpenMMCommon::DefaultReturn;
 }
+
+/**---------------------------------------------------------------------------------------
+
+   Get normally distributed random number
+
+   @return random value
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM SimTKOpenMMUtilities::getNormallyDistributedRandomNumber( void ) {
+    static bool nextValueIsValid = false;
+    static RealOpenMM nextValue = 0;
+    if (nextValueIsValid) {
+        nextValueIsValid = false;
+        return nextValue;
+    }
+    if (!_randomInitialized) {
+        init_gen_rand(_randomNumberSeed);
+        _randomInitialized = true;
+    }
+    
+    // Use the polar form of the Box-Muller transformation to generate two Gaussian random numbers.
+    
+    RealOpenMM x, y, r2;
+    do {
+        x = 2.0*genrand_real2()-1.0;
+        y = 2.0*genrand_real2()-1.0;
+        r2 = x*x + y*y;
+    } while (r2 >= 1.0 || r2 == 0.0);
+    RealOpenMM multiplier = sqrt((-2.0*log(r2))/r2);
+    nextValue = y*multiplier;
+    nextValueIsValid = true;
+    return x*multiplier;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get uniformly distributed random number in the range [0, 1)
+
+   @return random value
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber( void ) {
+    if (!_randomInitialized) {
+        init_gen_rand(_randomNumberSeed);
+        _randomInitialized = true;
+    }
+    return genrand_real2();
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get random number seed
+
+   @return random number seed
+
+   --------------------------------------------------------------------------------------- */
+
+uint32_t SimTKOpenMMUtilities::getRandomNumberSeed( void ) {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName  = "\nReferenceDynamics::getRandomNumberSeed";
+
+   // ---------------------------------------------------------------------------------------
+
+   return _randomNumberSeed;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set random number seed
+
+   @param seed    new seed value
+
+   @return ReferenceDynamics::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+void SimTKOpenMMUtilities::setRandomNumberSeed( uint32_t seed ) {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName  = "\nReferenceDynamics::setRandomNumberSeed";
+
+   // ---------------------------------------------------------------------------------------
+
+   _randomNumberSeed = seed;
+   init_gen_rand(_randomNumberSeed);
+}

platforms/reference/src/SimTKUtilities/SimTKOpenMMUtilities.h

@@ -52,6 +52,11 @@ bool checkString( T& t, const std::string& s, std::ios_base& (*f)(std::ios_base&
 
 class SimTKOpenMMUtilities {
 
+   private:
+
+       static uint32_t _randomNumberSeed;
+       static bool _randomInitialized;
+
    public:
 
       // file flag enums
@@ -576,6 +581,48 @@ class SimTKOpenMMUtilities {
                                  RealOpenMM** sumArray );
       
       
+      /**---------------------------------------------------------------------------------------
+      
+         Get normally distributed random number
+      
+         @return random value
+      
+         --------------------------------------------------------------------------------------- */
+      
+      static RealOpenMM getNormallyDistributedRandomNumber( void );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Get uniformly distributed random number in the range [0, 1)
+      
+         @return random value
+      
+         --------------------------------------------------------------------------------------- */
+      
+      static RealOpenMM getUniformlyDistributedRandomNumber( void );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Get random number seed
+      
+         @return random number seed
+      
+         --------------------------------------------------------------------------------------- */
+      
+      static uint32_t getRandomNumberSeed( void );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Set random number seed
+      
+         @param seed    new seed value
+      
+         @return DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+      
+      static void setRandomNumberSeed( uint32_t seed );
+      
 };
    
 // ---------------------------------------------------------------------------------------

platforms/reference/tests/TestReferenceAndersenThermostat.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * 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.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the reference implementation of AndersenThermostat.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "AndersenThermostat.h"
+#include "OpenMMContext.h"
+#include "ReferencePlatform.h"
+#include "StandardMMForceField.h"
+#include "System.h"
+#include "VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "../src/sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void testTemperature() {
+    const int numAtoms = 8;
+    const double temp = 100.0;
+    const double collisionFreq = 10.0;
+    ReferencePlatform platform;
+    System system(numAtoms, 0);
+    VerletIntegrator integrator(0.01);
+    StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0);
+    for (int i = 0; i < numAtoms; ++i) {
+        system.setAtomMass(i, 2.0);
+        forceField->setAtomParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
+    }
+    system.addForce(forceField);
+    AndersenThermostat* thermstat = new AndersenThermostat(temp, collisionFreq);
+    system.addForce(thermstat);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(numAtoms);
+    for (int i = 0; i < numAtoms; ++i)
+        positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
+    context.setPositions(positions);
+    
+    // Let it equilibrate.
+    
+    integrator.step(10000);
+    
+    // Now run it for a while and see if the temperature is correct.
+    
+    double ke = 0.0;
+    for (int i = 0; i < 1000; ++i) {
+        State state = context.getState(State::Energy);
+        ke += state.getKineticEnergy();
+        integrator.step(1);
+    }
+    ke /= 1000;
+    double expected = 0.5*numAtoms*3*BOLTZ*temp;
+    ASSERT_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
+}
+
+int main() {
+    try {
+        testTemperature();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}