Merge pull request #2384 from peastman/maxstep

Added max step size option to variable step size integrators

openmmapi/include/openmm/VariableLangevinIntegrator.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) 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.
+     */
+    void 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

@@ -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.
+     */
+    void 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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -205,6 +205,7 @@ UNITS = {
 ("*", "getSoluteDielectric") : (None, ()),
 ("*", "getSolventDielectric") : (None, ()),
 ("*", "getStepSize") : ("unit.picosecond", ()),
+("*", "getMaximumStepSize") : ("unit.picosecond", ()),
 ("*", "getSystem") : (None, ()),
 ("*", "getTabulatedFunction") : (None, ()),
 ("*", "getUseDispersionCorrection") : (None, ()),