Improved algorithm for VariableVerletIntegrator

openmmapi/include/openmm/VariableVerletIntegrator.h

@@ -63,10 +63,9 @@ public:
     /**
      * Create a VariableVerletIntegrator.
      *
-     * @param stepSize the initial step size to use (in picoseconds)
      * @param tol      the error tolerance
      */
-    VariableVerletIntegrator(double stepSize, double errorTol);
+    VariableVerletIntegrator(double errorTol);
     /**
      * Get the error tolerance.
      */

openmmapi/src/VariableVerletIntegrator.cpp

@@ -39,8 +39,7 @@ using namespace OpenMM;
 using std::string;
 using std::vector;
 
-VariableVerletIntegrator::VariableVerletIntegrator(double stepSize, double errorTol) : errorTol(errorTol) {
-    setStepSize(stepSize);
+VariableVerletIntegrator::VariableVerletIntegrator(double errorTol) : errorTol(errorTol) {
     setConstraintTolerance(1e-4);
 }
 

platforms/reference/src/ReferenceKernels.cpp

@@ -814,7 +814,7 @@ void ReferenceIntegrateVariableVerletStepKernel::execute(OpenMMContextImpl& cont
             delete dynamics;
             delete constraints;
         }
-        dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize), errorTol);
+        dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), (RealOpenMM) errorTol);
         vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
         findAnglesForCCMA(context.getSystem(), angles);
         constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
@@ -823,7 +823,7 @@ void ReferenceIntegrateVariableVerletStepKernel::execute(OpenMMContextImpl& cont
         prevErrorTol = errorTol;
     }
     dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
-    data.time += dynamics->getLastStepSize();
+    data.time += dynamics->getDeltaT();
 }
 
 ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {

platforms/reference/src/SimTKReference/ReferenceVariableVerletDynamics.cpp

@@ -40,9 +40,8 @@
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceVariableVerletDynamics::ReferenceVariableVerletDynamics( int numberOfAtoms,
-                                                          RealOpenMM deltaT, RealOpenMM accuracy ) :
-           ReferenceDynamics( numberOfAtoms, deltaT, 0.0 ), _accuracy(accuracy) {
+ReferenceVariableVerletDynamics::ReferenceVariableVerletDynamics( int numberOfAtoms, RealOpenMM accuracy ) :
+           ReferenceDynamics( numberOfAtoms, 0.0f, 0.0f ), _accuracy(accuracy) {
 
    // ---------------------------------------------------------------------------------------
 
@@ -96,18 +95,6 @@ void ReferenceVariableVerletDynamics::setAccuracy( RealOpenMM accuracy ) {
     _accuracy = accuracy;
 }
 
-/**---------------------------------------------------------------------------------------
-
-   Get the actual size of the last step that was taken
-
-   @return step size
-
-   --------------------------------------------------------------------------------------- */
-
-RealOpenMM ReferenceVariableVerletDynamics::getLastStepSize( void ) const {
-    return _lastStepSize;
-}
-
 /**---------------------------------------------------------------------------------------
 
    Print parameters
@@ -195,58 +182,34 @@ int ReferenceVariableVerletDynamics::update( int numberOfAtoms, RealOpenMM** ato
        }
     }
 
-    // Try different step sizes until the accuracy is acceptable.
-
-    bool success = false;
-    RealOpenMM maxStepSize = 5.0f*getDeltaT();
-    while (!success) {
-        // Perform the integration and estimate the error.
-
-        _lastStepSize = getDeltaT();
     RealOpenMM error = zero;
     for (int i = 0; i < numberOfAtoms; ++i) {
         for (int j = 0; j < 3; ++j) {
-                RealOpenMM xref = atomCoordinates[i][j] + velocities[i][j]*getDeltaT();
-                RealOpenMM vPrime = velocities[i][j] + inverseMasses[i]*forces[i][j]*getDeltaT();
-                xPrime[i][j] = atomCoordinates[i][j] + vPrime*getDeltaT();
-                RealOpenMM xerror = xPrime[i][j]-xref;
+            RealOpenMM xerror = inverseMasses[i]*forces[i][j];
             error += xerror*xerror;
         }
     }
     error = SQRT(error/(numberOfAtoms*3));
-
-        // Select a new step size.
-
-        const RealOpenMM Safety = 0.9f, MinShrink = 0.1f;
-        const RealOpenMM HysteresisLow = 0.9f, HysteresisHigh = 1.0f, ErrorOrder = 2.0f;
-
-        RealOpenMM newStepSize = Safety*getDeltaT()*POW(getAccuracy()/error, 1.0f/ErrorOrder);
-        if (newStepSize > getDeltaT()) {
-            if (newStepSize < HysteresisHigh*getDeltaT())
-                newStepSize = getDeltaT();
+    RealOpenMM newStepSize = SQRT(getAccuracy()/error);
+    if (getDeltaT() > 0.0f)
+        newStepSize = std::min(newStepSize, getDeltaT()*2.0f); // For safety, limit how quickly dt can increase.
+    if (newStepSize > getDeltaT() && newStepSize < 1.2f*getDeltaT())
+        newStepSize = getDeltaT(); // Keeping dt constant between steps improves the behavior of the integrator.
+    RealOpenMM vstep = 0.5f*(newStepSize+getDeltaT()); // The time interval by which to advance the velocities
+    setDeltaT(newStepSize);
+    for (int i = 0; i < numberOfAtoms; ++i) {
+        for (int j = 0; j < 3; ++j) {
+            RealOpenMM vPrime = velocities[i][j] + inverseMasses[i]*forces[i][j]*vstep;
+            xPrime[i][j] = atomCoordinates[i][j] + vPrime*getDeltaT();
+        }
     }
-        if (newStepSize < getDeltaT() && error <= getAccuracy())
-            newStepSize = getDeltaT();
-        newStepSize = std::min(newStepSize, maxStepSize);
-        newStepSize = std::max(newStepSize, MinShrink*getDeltaT());
-        if (newStepSize < getDeltaT())
-            newStepSize = std::min(newStepSize, HysteresisLow*getDeltaT());
-        success = (newStepSize >= getDeltaT());
-        if (success) {
     ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
     if (referenceConstraintAlgorithm)
-                success = (referenceConstraintAlgorithm->apply(numberOfAtoms, atomCoordinates, xPrime, inverseMasses) != ErrorReturn);
-            if (!success) {
-                newStepSize *= 0.5f;
-                maxStepSize = newStepSize;
-            }
-        }
-        setDeltaT(newStepSize);
-   }
+        referenceConstraintAlgorithm->apply(numberOfAtoms, atomCoordinates, xPrime, inverseMasses);
 
    // Update the positions and velocities.
 
-   RealOpenMM velocityScale = static_cast<RealOpenMM>( 1.0/_lastStepSize );
+   RealOpenMM velocityScale = one/getDeltaT();
    for (int i = 0; i < numberOfAtoms; ++i) {
        for (int j = 0; j < 3; ++j) {
            velocities[i][j] = velocityScale*(xPrime[i][j] - atomCoordinates[i][j]);

platforms/reference/src/SimTKReference/ReferenceVariableVerletDynamics.h

@@ -35,7 +35,7 @@ class ReferenceVariableVerletDynamics : public ReferenceDynamics {
 
       enum TwoDArrayIndicies { xPrime2D, Max2DArrays };
       enum OneDArrayIndicies { InverseMasses, Max1DArrays };
-      RealOpenMM _accuracy, _lastStepSize;
+      RealOpenMM _accuracy;
 
    public:
 
@@ -44,12 +44,11 @@ class ReferenceVariableVerletDynamics : public ReferenceDynamics {
          Constructor
 
          @param numberOfAtoms  number of atoms
-         @param deltaT         initial delta t for dynamics
          @param accuracy       required accuracy
 
          --------------------------------------------------------------------------------------- */
 
-       ReferenceVariableVerletDynamics( int numberOfAtoms, RealOpenMM deltaT, RealOpenMM accuracy );
+       ReferenceVariableVerletDynamics( int numberOfAtoms, RealOpenMM accuracy );
 
       /**---------------------------------------------------------------------------------------
 
@@ -77,16 +76,6 @@ class ReferenceVariableVerletDynamics : public ReferenceDynamics {
 
       void setAccuracy( RealOpenMM accuracy );
 
-      /**---------------------------------------------------------------------------------------
-
-         Get the actual size of the last step that was taken
-
-         @return step size
-
-         --------------------------------------------------------------------------------------- */
-
-      RealOpenMM getLastStepSize( void ) const;
-
       /**---------------------------------------------------------------------------------------
 
          Print parameters

platforms/reference/tests/TestReferenceVariableVerletIntegrator.cpp

@@ -55,7 +55,7 @@ void testSingleBond() {
     System system;
     system.addParticle(2.0);
     system.addParticle(2.0);
-    VariableVerletIntegrator integrator(0.01, 1e-6);
+    VariableVerletIntegrator integrator(1e-6);
     HarmonicBondForce* forceField = new HarmonicBondForce();
     forceField->addBond(0, 1, 1.5, 1);
     system.addForce(forceField);
@@ -91,7 +91,7 @@ void testConstraints() {
     const double temp = 500.0;
     ReferencePlatform platform;
     System system;
-    VariableVerletIntegrator integrator(0.002, 1e-5);
+    VariableVerletIntegrator integrator(1e-5);
     integrator.setConstraintTolerance(1e-5);
     NonbondedForce* forceField = new NonbondedForce();
     for (int i = 0; i < numParticles; ++i) {
@@ -141,7 +141,7 @@ void testConstrainedClusters() {
     const double temp = 500.0;
     ReferencePlatform platform;
     System system;
-    VariableVerletIntegrator integrator(0.002, 1e-5);
+    VariableVerletIntegrator integrator(1e-5);
     integrator.setConstraintTolerance(1e-5);
     NonbondedForce* forceField = new NonbondedForce();
     for (int i = 0; i < numParticles; ++i) {