Skip to content

GitLab

Commit 0620a02f authored by peastman's avatar peastman
Browse files

Added max step size option to variable step size integrators

parent 68424035
Hide whitespace changes
Inline Side-by-side
Showing with 126 additions and 32 deletions
openmmapi/include/openmm/VariableLangevinIntegrator.h
View file @ 0620a02f
......@@ -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) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2019 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -39,7 +39,7 @@
namespace OpenMM {
/**
* This is an error contolled, variable time step Integrator that simulates a System using Langevin
* This is an error controlled, variable time step Integrator that simulates a System using Langevin
* dynamics. It compares the result of the Langevin integrator to that of an
* explicit Euler integrator, takes the difference between the two as a measure of the integration
* error in each time step, and continuously adjusts the step size to keep the error below a
......@@ -50,6 +50,10 @@ namespace OpenMM {
* adjustable parameter that affects the step size and integration accuracy. You
* should try different values to find the largest one that produces a trajectory sufficiently
* accurate for your purposes. 0.001 is often a good starting point.
*
* You can optionally set a maximum step size it will ever use. This is useful to prevent it
* from taking excessively large steps in usual situations, such as when the system is right at
* a local energy minimum.
*/
class OPENMM_EXPORT VariableLangevinIntegrator : public Integrator {
......@@ -108,6 +112,20 @@ public:
void setErrorTolerance(double tol) {
errorTol = tol;
}
/**
* Get the maximum step size the integrator will ever use, in ps. If this
* is 0 (the default), no limit will be applied to step sizes.
*/
double getMaximumStepSize() const {
return maxStepSize;
}
/**
* Set the maximum step size the integrator will ever use, in ps. If this
* is 0 (the default), no limit will be applied to step sizes.
*/
double setMaximumStepSize(double size) {
maxStepSize = size;
}
/**
* Get the random number seed. See setRandomNumberSeed() for details.
*/
......@@ -165,7 +183,7 @@ protected:
*/
double computeKineticEnergy();
private:
double temperature, friction, errorTol;
double temperature, friction, errorTol, maxStepSize;
int randomNumberSeed;
Kernel kernel;
};
......
openmmapi/include/openmm/VariableVerletIntegrator.h
View file @ 0620a02f
......@@ -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) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2019 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -39,7 +39,7 @@
namespace OpenMM {
/**
* This is an error contolled, variable time step Integrator that simulates a System using the
* This is an error controlled, variable time step Integrator that simulates a System using the
* leap-frog Verlet algorithm. It compares the result of the Verlet integrator to that of an
* explicit Euler integrator, takes the difference between the two as a measure of the integration
* error in each time step, and continuously adjusts the step size to keep the error below a
......@@ -56,6 +56,10 @@ namespace OpenMM {
* This makes it most appropriate for constant temperate simulations. In constant energy simulations
* where precise energy conservation over long time periods is important, a fixed step size Verlet
* integrator may be more appropriate.
*
* You can optionally set a maximum step size it will ever use. This is useful to prevent it
* from taking excessively large steps in usual situations, such as when the system is right at
* a local energy minimum.
*/
class OPENMM_EXPORT VariableVerletIntegrator : public Integrator {
......@@ -78,6 +82,20 @@ public:
void setErrorTolerance(double tol) {
errorTol = tol;
}
/**
* Get the maximum step size the integrator will ever use, in ps. If this
* is 0 (the default), no limit will be applied to step sizes.
*/
double getMaximumStepSize() const {
return maxStepSize;
}
/**
* Set the maximum step size the integrator will ever use, in ps. If this
* is 0 (the default), no limit will be applied to step sizes.
*/
double setMaximumStepSize(double size) {
maxStepSize = size;
}
/**
* Advance a simulation through time by taking a series of time steps.
*
......@@ -114,7 +132,7 @@ protected:
*/
double computeKineticEnergy();
private:
double errorTol;
double errorTol, maxStepSize;
Kernel kernel;
};
......
openmmapi/src/VariableLangevinIntegrator.cpp
View file @ 0620a02f
......@@ -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-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2019 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -46,6 +46,7 @@ VariableLangevinIntegrator::VariableLangevinIntegrator(double temperature, doubl
setTemperature(temperature);
setFriction(frictionCoeff);
setErrorTolerance(errorTol);
setMaximumStepSize(0.0);
setConstraintTolerance(1e-5);
setRandomNumberSeed(0);
setStepSize(0.0);
......
openmmapi/src/VariableVerletIntegrator.cpp
View file @ 0620a02f
......@@ -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-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2019 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -41,7 +41,7 @@ using namespace OpenMM;
using std::string;
using std::vector;
VariableVerletIntegrator::VariableVerletIntegrator(double errorTol) : errorTol(errorTol) {
VariableVerletIntegrator::VariableVerletIntegrator(double errorTol) : errorTol(errorTol), maxStepSize(0.0) {
setConstraintTolerance(1e-5);
setStepSize(0.0);
}
......
platforms/cuda/src/CudaKernels.cpp
View file @ 0620a02f
......@@ -7221,6 +7221,8 @@ double CudaIntegrateVariableVerletStepKernel::execute(ContextImpl& context, cons
// Select the step size to use.
double maxStepSize = maxTime-cu.getTime();
if (integrator.getMaximumStepSize() > 0)
maxStepSize = min(integrator.getMaximumStepSize(), maxStepSize);
float maxStepSizeFloat = (float) maxStepSize;
double tol = integrator.getErrorTolerance();
float tolFloat = (float) tol;
......@@ -7295,6 +7297,8 @@ double CudaIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, co
// Select the step size to use.
double maxStepSize = maxTime-cu.getTime();
if (integrator.getMaximumStepSize() > 0)
maxStepSize = min(integrator.getMaximumStepSize(), maxStepSize);
float maxStepSizeFloat = (float) maxStepSize;
double tol = integrator.getErrorTolerance();
float tolFloat = (float) tol;
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ 0620a02f
......@@ -7594,6 +7594,8 @@ double OpenCLIntegrateVariableVerletStepKernel::execute(ContextImpl& context, co
// Select the step size to use.
double maxStepSize = maxTime-cl.getTime();
if (integrator.getMaximumStepSize() > 0)
maxStepSize = min(integrator.getMaximumStepSize(), maxStepSize);
float maxStepSizeFloat = (float) maxStepSize;
if (useDouble) {
selectSizeKernel.setArg<cl_double>(1, maxStepSize);
......@@ -7691,6 +7693,8 @@ double OpenCLIntegrateVariableLangevinStepKernel::execute(ContextImpl& context,
// Select the step size to use.
double maxStepSize = maxTime-cl.getTime();
if (integrator.getMaximumStepSize() > 0)
maxStepSize = min(integrator.getMaximumStepSize(), maxStepSize);
float maxStepSizeFloat = (float) maxStepSize;
if (useDouble) {
selectSizeKernel.setArg<cl_double>(0, maxStepSize);
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 0620a02f
......@@ -2225,6 +2225,8 @@ double ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& contex
prevErrorTol = errorTol;
}
double maxStepSize = maxTime-data.time;
if (integrator.getMaximumStepSize() > 0)
maxStepSize = min(integrator.getMaximumStepSize(), maxStepSize);
dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize, integrator.getConstraintTolerance());
data.time += dynamics->getDeltaT();
if (dynamics->getDeltaT() == maxStepSize)
......@@ -2264,6 +2266,8 @@ double ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context,
prevErrorTol = errorTol;
}
double maxStepSize = maxTime-data.time;
if (integrator.getMaximumStepSize() > 0)
maxStepSize = min(integrator.getMaximumStepSize(), maxStepSize);
dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize, integrator.getConstraintTolerance());
data.time += dynamics->getDeltaT();
if (dynamics->getDeltaT() == maxStepSize)
......
serialization/src/VariableLangevinIntegratorProxy.cpp
View file @ 0620a02f
......@@ -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 Stanford University and the Authors. *
* Portions copyright (c) 2010-2019 Stanford University and the Authors. *
* Authors: Peter Eastman, Yutong Zhao *
* Contributors: *
* *
......@@ -36,22 +36,23 @@ using namespace std;
using namespace OpenMM;
VariableLangevinIntegratorProxy::VariableLangevinIntegratorProxy() : SerializationProxy("VariableLangevinIntegrator") {
}
void VariableLangevinIntegratorProxy::serialize(const void* object, SerializationNode& node) const {
node.setIntProperty("version", 1);
node.setIntProperty("version", 2);
const VariableLangevinIntegrator& integrator = *reinterpret_cast<const VariableLangevinIntegrator*>(object);
node.setDoubleProperty("stepSize", integrator.getStepSize());
node.setDoubleProperty("constraintTolerance", integrator.getConstraintTolerance());
node.setDoubleProperty("temperature", integrator.getTemperature());
node.setDoubleProperty("friction", integrator.getFriction());
node.setDoubleProperty("errorTol", integrator.getErrorTolerance());
node.setDoubleProperty("maxStepSize", integrator.getMaximumStepSize());
node.setIntProperty("randomSeed", integrator.getRandomNumberSeed());
}
void* VariableLangevinIntegratorProxy::deserialize(const SerializationNode& node) const {
if (node.getIntProperty("version") != 1)
int version = node.getIntProperty("version");
if (version < 1 || version > 2)
throw OpenMMException("Unsupported version number");
VariableLangevinIntegrator *integrator = new VariableLangevinIntegrator(node.getDoubleProperty("temperature"),
node.getDoubleProperty("friction"),
......@@ -59,5 +60,7 @@ void* VariableLangevinIntegratorProxy::deserialize(const SerializationNode& node
integrator->setStepSize(node.getDoubleProperty("stepSize"));
integrator->setConstraintTolerance(node.getDoubleProperty("constraintTolerance"));
integrator->setRandomNumberSeed(node.getIntProperty("randomSeed"));
if (version > 1)
integrator->setMaximumStepSize(node.getDoubleProperty("maxStepSize"));
return integrator;
}
\ No newline at end of file
serialization/src/VariableVerletIntegratorProxy.cpp
View file @ 0620a02f
......@@ -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 Stanford University and the Authors. *
* Portions copyright (c) 2010-2019 Stanford University and the Authors. *
* Authors: Peter Eastman, Yutong Zhao *
* Contributors: *
* *
......@@ -36,22 +36,25 @@ using namespace std;
using namespace OpenMM;
VariableVerletIntegratorProxy::VariableVerletIntegratorProxy() : SerializationProxy("VariableVerletIntegrator") {
}
void VariableVerletIntegratorProxy::serialize(const void* object, SerializationNode& node) const {
node.setIntProperty("version", 1);
node.setIntProperty("version", 2);
const VariableVerletIntegrator& integrator = *reinterpret_cast<const VariableVerletIntegrator*>(object);
node.setDoubleProperty("errorTol", integrator.getErrorTolerance());
node.setDoubleProperty("maxStepSize", integrator.getMaximumStepSize());
node.setDoubleProperty("stepSize", integrator.getStepSize());
node.setDoubleProperty("constraintTolerance", integrator.getConstraintTolerance());
}
void* VariableVerletIntegratorProxy::deserialize(const SerializationNode& node) const {
if (node.getIntProperty("version") != 1)
int version = node.getIntProperty("version");
if (version < 1 || version > 2)
throw OpenMMException("Unsupported version number");
VariableVerletIntegrator *integrator = new VariableVerletIntegrator(node.getDoubleProperty("errorTol"));
integrator->setStepSize(node.getDoubleProperty("stepSize"));
integrator->setConstraintTolerance(node.getDoubleProperty("constraintTolerance"));
if (version > 1)
integrator->setMaximumStepSize(node.getDoubleProperty("maxStepSize"));
return integrator;
}
\ No newline at end of file
serialization/tests/TestSerializeIntegrator.cpp
View file @ 0620a02f
......@@ -89,29 +89,31 @@ void testSerializeBrownianIntegrator() {
}
void testSerializeVariableVerletIntegrator() {
VariableVerletIntegrator *intg = new VariableVerletIntegrator(0.04234);
VariableVerletIntegrator intg(0.04234);
intg.setMaximumStepSize(0.32);
stringstream ss;
XmlSerializer::serialize<Integrator>(intg, "VariableVerletIntegrator", ss);
XmlSerializer::serialize<Integrator>(&intg, "VariableVerletIntegrator", ss);
VariableVerletIntegrator *intg2 = dynamic_cast<VariableVerletIntegrator*>(XmlSerializer::deserialize<Integrator>(ss));
ASSERT_EQUAL(intg->getConstraintTolerance(), intg2->getConstraintTolerance());
ASSERT_EQUAL(intg->getStepSize(), intg2->getStepSize());
ASSERT_EQUAL(intg->getErrorTolerance(), intg2->getErrorTolerance());
delete intg;
ASSERT_EQUAL(intg.getConstraintTolerance(), intg2->getConstraintTolerance());
ASSERT_EQUAL(intg.getStepSize(), intg2->getStepSize());
ASSERT_EQUAL(intg.getErrorTolerance(), intg2->getErrorTolerance());
ASSERT_EQUAL(intg.getMaximumStepSize(), intg2->getMaximumStepSize());
delete intg2;
}
void testSerializeVariableLangevinIntegrator() {
VariableLangevinIntegrator *intg = new VariableLangevinIntegrator(243.1, 3.234, 0.0021);
VariableLangevinIntegrator intg(243.1, 3.234, 0.0021);
intg.setMaximumStepSize(0.32);
stringstream ss;
XmlSerializer::serialize<Integrator>(intg, "VariableLangevinIntegrator", ss);
XmlSerializer::serialize<Integrator>(&intg, "VariableLangevinIntegrator", ss);
VariableLangevinIntegrator *intg2 = dynamic_cast<VariableLangevinIntegrator*>(XmlSerializer::deserialize<Integrator>(ss));
ASSERT_EQUAL(intg->getConstraintTolerance(), intg2->getConstraintTolerance());
ASSERT_EQUAL(intg->getStepSize(), intg2->getStepSize());
ASSERT_EQUAL(intg->getErrorTolerance(), intg2->getErrorTolerance());
ASSERT_EQUAL(intg->getFriction(), intg2->getFriction());
ASSERT_EQUAL(intg->getTemperature(), intg2->getTemperature());
ASSERT_EQUAL(intg->getRandomNumberSeed(), intg2->getRandomNumberSeed());
delete intg;
ASSERT_EQUAL(intg.getConstraintTolerance(), intg2->getConstraintTolerance());
ASSERT_EQUAL(intg.getStepSize(), intg2->getStepSize());
ASSERT_EQUAL(intg.getErrorTolerance(), intg2->getErrorTolerance());
ASSERT_EQUAL(intg.getFriction(), intg2->getFriction());
ASSERT_EQUAL(intg.getTemperature(), intg2->getTemperature());
ASSERT_EQUAL(intg.getRandomNumberSeed(), intg2->getRandomNumberSeed());
ASSERT_EQUAL(intg.getMaximumStepSize(), intg2->getMaximumStepSize());
delete intg2;
}
......
tests/TestVariableLangevinIntegrator.h
View file @ 0620a02f
......@@ -310,6 +310,24 @@ void testArgonBox() {
ke /= 1000;
double expected = 1.5 * numParticles * BOLTZ * temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.01);
// Compute the mean step size.
double meanStep = 0;
for (int i = 0; i < 100; i++) {
integrator.step(1);
meanStep += integrator.getStepSize();
}
meanStep /= 100;
double maxStep = meanStep/2;
// Now set a limit on the step size and see if it is obeyed.
integrator.setMaximumStepSize(maxStep);
for (int i = 0; i < 100; i++) {
integrator.step(1);
ASSERT(integrator.getStepSize() <= maxStep);
}
}
void runPlatformTests();
......
tests/TestVariableVerletIntegrator.h
View file @ 0620a02f
......@@ -289,6 +289,24 @@ void testArgonBox() {
double energy = state.getKineticEnergy() + state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
}
// Compute the mean step size.
double meanStep = 0;
for (int i = 0; i < 100; i++) {
integrator.step(1);
meanStep += integrator.getStepSize();
}
meanStep /= 100;
double maxStep = meanStep/2;
// Now set a limit on the step size and see if it is obeyed.
integrator.setMaximumStepSize(maxStep);
for (int i = 0; i < 100; i++) {
integrator.step(1);
ASSERT(integrator.getStepSize() <= maxStep);
}
}
void runPlatformTests();
......
wrappers/python/src/swig_doxygen/swigInputConfig.py
View file @ 0620a02f
......@@ -205,6 +205,7 @@ UNITS = {
("*", "getSoluteDielectric") : (None, ()),
("*", "getSolventDielectric") : (None, ()),
("*", "getStepSize") : ("unit.picosecond", ()),
("*", "getMaximumStepSize") : ("unit.picosecond", ()),
("*", "getSystem") : (None, ()),
("*", "getTabulatedFunction") : (None, ()),
("*", "getUseDispersionCorrection") : (None, ()),
......
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