Skip to content

GitLab

Commit d9c27051 authored by peastman's avatar peastman
Browse files

Barostats use checkpointed random number generator

parent e2629476
Hide whitespace changes
Inline Side-by-side
Showing with 20 additions and 47 deletions
openmmapi/include/openmm/internal/MonteCarloAnisotropicBarostatImpl.h
View file @ d9c27051
......@@ -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
View file @ d9c27051
......@@ -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
View file @ d9c27051
......@@ -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
View file @ d9c27051
......@@ -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
View file @ d9c27051
......@@ -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
View file @ d9c27051
......@@ -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);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment