Merge pull request #2400 from peastman/random

Barostats use checkpointed random number generator

openmmapi/include/openmm/internal/MonteCarloAnisotropicBarostatImpl.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010 Stanford University and the Authors.           *
+ * Portions copyright (c) 2010-2019 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,7 +35,6 @@
 #include "ForceImpl.h"
 #include "openmm/MonteCarloAnisotropicBarostat.h"
 #include "openmm/Kernel.h"
-#include "sfmt/SFMT.h"
 #include <string>
 
 namespace OpenMM {
@@ -62,7 +61,6 @@ private:
     const MonteCarloAnisotropicBarostat& owner;
     int step, numAttempted[3], numAccepted[3];
     double volumeScale[3];
-    OpenMM_SFMT::SFMT random;
     Kernel kernel;
 };
 

openmmapi/include/openmm/internal/MonteCarloBarostatImpl.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010 Stanford University and the Authors.           *
+ * Portions copyright (c) 2010-2019 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,7 +35,6 @@
 #include "ForceImpl.h"
 #include "openmm/MonteCarloBarostat.h"
 #include "openmm/Kernel.h"
-#include "sfmt/SFMT.h"
 #include <string>
 
 namespace OpenMM {
@@ -62,7 +61,6 @@ private:
     const MonteCarloBarostat& owner;
     int step, numAttempted, numAccepted;
     double volumeScale;
-    OpenMM_SFMT::SFMT random;
     Kernel kernel;
 };
 

openmmapi/include/openmm/internal/MonteCarloMembraneBarostatImpl.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010-2014 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2019 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,7 +35,6 @@
 #include "ForceImpl.h"
 #include "openmm/MonteCarloMembraneBarostat.h"
 #include "openmm/Kernel.h"
-#include "sfmt/SFMT.h"
 #include <string>
 
 namespace OpenMM {
@@ -62,7 +61,6 @@ private:
     const MonteCarloMembraneBarostat& owner;
     int step, numAttempted[3], numAccepted[3];
     double volumeScale[3];
-    OpenMM_SFMT::SFMT random;
     Kernel kernel;
 };
 

openmmapi/src/MonteCarloAnisotropicBarostatImpl.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) 2010-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2019 Stanford University and the Authors.      *
  * Authors: Peter Eastman, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,6 +35,7 @@
 #include "openmm/Context.h"
 #include "openmm/kernels.h"
 #include "openmm/OpenMMException.h"
+#include "SimTKOpenMMUtilities.h"
 #include <cmath>
 #include <vector>
 #include <algorithm>
@@ -43,11 +44,6 @@ using namespace OpenMM;
 using namespace OpenMM_SFMT;
 using std::vector;
 
-const float BOLTZMANN = 1.380658e-23f; // (J/K)
-const float AVOGADRO = 6.0221367e23f;
-const float RGAS = BOLTZMANN*AVOGADRO; // (J/(mol K))
-const float BOLTZ = RGAS/1000;         // (kJ/(mol K))
-
 MonteCarloAnisotropicBarostatImpl::MonteCarloAnisotropicBarostatImpl(const MonteCarloAnisotropicBarostat& owner) : owner(owner), step(0) {
 }
 
@@ -64,10 +60,7 @@ void MonteCarloAnisotropicBarostatImpl::initialize(ContextImpl& context) {
         numAttempted[i] = 0;
         numAccepted[i] = 0;
     }
-    int randSeed = owner.getRandomNumberSeed();
-    // A random seed of 0 means use a unique one
-    if (randSeed == 0) randSeed = osrngseed();
-    init_gen_rand(randSeed, random);
+    SimTKOpenMMUtilities::setRandomNumberSeed(owner.getRandomNumberSeed());
 }
 
 void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context) {
@@ -85,7 +78,7 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
     // Choose which axis to modify at random.
     int axis;
     while (true) {
-        double rnd = genrand_real2(random)*3.0;
+        double rnd = SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber()*3.0;
         if (rnd < 1.0) {
             if (owner.getScaleX()) {
                 axis = 0;
@@ -110,7 +103,7 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
     Vec3 box[3];
     context.getPeriodicBoxVectors(box[0], box[1], box[2]);
     double volume = box[0][0]*box[1][1]*box[2][2];
-    double deltaVolume = volumeScale[axis]*2*(genrand_real2(random)-0.5);
+    double deltaVolume = volumeScale[axis]*2*(SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber()-0.5);
     double newVolume = volume+deltaVolume;
     Vec3 lengthScale(1.0, 1.0, 1.0);
     lengthScale[axis] = newVolume/volume;
@@ -124,7 +117,7 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
     double finalEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
     double kT = BOLTZ*context.getParameter(MonteCarloAnisotropicBarostat::Temperature());
     double w = finalEnergy-initialEnergy + pressure*deltaVolume - context.getMolecules().size()*kT*std::log(newVolume/volume);
-    if (w > 0 && genrand_real2(random) > std::exp(-w/kT)) {
+    if (w > 0 && SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber() > std::exp(-w/kT)) {
         // Reject the step.
         
         kernel.getAs<ApplyMonteCarloBarostatKernel>().restoreCoordinates(context);

openmmapi/src/MonteCarloBarostatImpl.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) 2010-2017 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2019 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,6 +35,7 @@
 #include "openmm/Context.h"
 #include "openmm/kernels.h"
 #include "openmm/OpenMMException.h"
+#include "SimTKOpenMMUtilities.h"
 #include <cmath>
 #include <vector>
 #include <algorithm>
@@ -43,11 +44,6 @@ using namespace OpenMM;
 using namespace OpenMM_SFMT;
 using std::vector;
 
-const float BOLTZMANN = 1.380658e-23f; // (J/K)
-const float AVOGADRO = 6.0221367e23f;
-const float RGAS = BOLTZMANN*AVOGADRO; // (J/(mol K))
-const float BOLTZ = RGAS/1000;         // (kJ/(mol K))
-
 MonteCarloBarostatImpl::MonteCarloBarostatImpl(const MonteCarloBarostat& owner) : owner(owner), step(0) {
 }
 
@@ -62,10 +58,7 @@ void MonteCarloBarostatImpl::initialize(ContextImpl& context) {
     volumeScale = 0.01*volume;
     numAttempted = 0;
     numAccepted = 0;
-    int randSeed = owner.getRandomNumberSeed();
-    // A random seed of 0 means use a unique one
-    if (randSeed == 0) randSeed = osrngseed();
-    init_gen_rand(randSeed, random);
+    SimTKOpenMMUtilities::setRandomNumberSeed(owner.getRandomNumberSeed());
 }
 
 void MonteCarloBarostatImpl::updateContextState(ContextImpl& context, bool& forcesInvalid) {
@@ -82,7 +75,7 @@ void MonteCarloBarostatImpl::updateContextState(ContextImpl& context, bool& forc
     Vec3 box[3];
     context.getPeriodicBoxVectors(box[0], box[1], box[2]);
     double volume = box[0][0]*box[1][1]*box[2][2];
-    double deltaVolume = volumeScale*2*(genrand_real2(random)-0.5);
+    double deltaVolume = volumeScale*2*(SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber()-0.5);
     double newVolume = volume+deltaVolume;
     double lengthScale = std::pow(newVolume/volume, 1.0/3.0);
     kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale, lengthScale, lengthScale);
@@ -94,7 +87,7 @@ void MonteCarloBarostatImpl::updateContextState(ContextImpl& context, bool& forc
     double pressure = context.getParameter(MonteCarloBarostat::Pressure())*(AVOGADRO*1e-25);
     double kT = BOLTZ*context.getParameter(MonteCarloBarostat::Temperature());
     double w = finalEnergy-initialEnergy + pressure*deltaVolume - context.getMolecules().size()*kT*std::log(newVolume/volume);
-    if (w > 0 && genrand_real2(random) > std::exp(-w/kT)) {
+    if (w > 0 && SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber() > std::exp(-w/kT)) {
         // Reject the step.
 
         kernel.getAs<ApplyMonteCarloBarostatKernel>().restoreCoordinates(context);

openmmapi/src/MonteCarloMembraneBarostatImpl.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) 2010-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2019 Stanford University and the Authors.      *
  * Authors: Peter Eastman, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,6 +35,7 @@
 #include "openmm/Context.h"
 #include "openmm/kernels.h"
 #include "openmm/OpenMMException.h"
+#include "SimTKOpenMMUtilities.h"
 #include <cmath>
 #include <vector>
 #include <algorithm>
@@ -43,11 +44,6 @@ using namespace OpenMM;
 using namespace OpenMM_SFMT;
 using std::vector;
 
-const float BOLTZMANN = 1.380658e-23f; // (J/K)
-const float AVOGADRO = 6.0221367e23f;
-const float RGAS = BOLTZMANN*AVOGADRO; // (J/(mol K))
-const float BOLTZ = RGAS/1000;         // (kJ/(mol K))
-
 MonteCarloMembraneBarostatImpl::MonteCarloMembraneBarostatImpl(const MonteCarloMembraneBarostat& owner) : owner(owner), step(0) {
 }
 
@@ -64,10 +60,7 @@ void MonteCarloMembraneBarostatImpl::initialize(ContextImpl& context) {
         numAttempted[i] = 0;
         numAccepted[i] = 0;
     }
-    int randSeed = owner.getRandomNumberSeed();
-    // A random seed of 0 means use a unique one
-    if (randSeed == 0) randSeed = osrngseed();
-    init_gen_rand(randSeed, random);
+    SimTKOpenMMUtilities::setRandomNumberSeed(owner.getRandomNumberSeed());
 }
 
 void MonteCarloMembraneBarostatImpl::updateContextState(ContextImpl& context) {
@@ -84,7 +77,7 @@ void MonteCarloMembraneBarostatImpl::updateContextState(ContextImpl& context) {
     // Choose which axis to modify at random.
     int axis;
     while (true) {
-        double rnd = genrand_real2(random)*3.0;
+        double rnd = SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber()*3.0;
         if (rnd < 1.0) {
             axis = 0;
             break;
@@ -102,7 +95,7 @@ void MonteCarloMembraneBarostatImpl::updateContextState(ContextImpl& context) {
     Vec3 box[3];
     context.getPeriodicBoxVectors(box[0], box[1], box[2]);
     double volume = box[0][0]*box[1][1]*box[2][2];
-    double deltaVolume = volumeScale[axis]*2*(genrand_real2(random)-0.5);
+    double deltaVolume = volumeScale[axis]*2*(SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber()-0.5);
     double newVolume = volume+deltaVolume;
     Vec3 lengthScale(1.0, 1.0, 1.0);
     if ((axis == 0 || axis == 1) && owner.getXYMode() == MonteCarloMembraneBarostat::XYIsotropic)
@@ -125,7 +118,7 @@ void MonteCarloMembraneBarostatImpl::updateContextState(ContextImpl& context) {
     double finalEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
     double kT = BOLTZ*context.getParameter(MonteCarloMembraneBarostat::Temperature());
     double w = finalEnergy-initialEnergy + pressure*deltaVolume - tension*deltaArea - context.getMolecules().size()*kT*std::log(newVolume/volume);
-    if (w > 0 && genrand_real2(random) > std::exp(-w/kT)) {
+    if (w > 0 && SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber() > std::exp(-w/kT)) {
         // Reject the step.
         
         kernel.getAs<ApplyMonteCarloBarostatKernel>().restoreCoordinates(context);