Fixed bugs in CustomIntegrator, including incorrect examples of velocity verlet with constraints

openmmapi/include/openmm/CustomIntegrator.h

@@ -144,27 +144,19 @@ namespace OpenMM {
  * time step so that, for example, an AndersenThermostat can randomize velocities
  * or a MonteCarloBarostat can scale particle positions.  You need to add a
  * step to tell the integrator when to do this.  The following example corrects
- * both these problems:
+ * both these problems, using the RATTLE algorithm to apply constraints:
  * 
  * <tt><pre>
  * CustomIntegrator integrator(0.001);
+ * integrator.addPerDofVariable("x1", 0);
+ * integrator.addUpdateContextState();
  * integrator.addComputePerDof("v", "v+0.5*dt*f/m");
  * integrator.addComputePerDof("x", "x+dt*v");
  * integrator.addConstrainPositions();
- * integrator.addUpdateContextState();
- * integrator.addComputePerDof("v", "v+0.5*dt*f/m");
+ * integrator.addComputePerDof("v", "v+0.5*dt*f/m+(x-x1)/dt");
  * integrator.addConstrainVelocities();
  * </pre></tt>
  * 
- * This integrator includes two steps that require forces (the two velocity
- * updates) and three steps that can potentially change particle positions and
- * thus invalidate the forces (the position update, applying position constraints,
- * and allowing Forces to update the context state).  We put all three of these
- * steps together to minimize the number of force computations needed.  If we had
- * put addUpdateContextState() at the beginning of the algorithm instead, that would
- * risk invalidating the forces just before the first velocity update, thus
- * requiring two force evaluations per time step instead of one.
- * 
  * CustomIntegrator can be used to implement multiple time step integrators.  The
  * following example shows an r-RESPA integrator.  It assumes the quickly changing
  * forces are in force group 0 and the slowly changing ones are in force group 1.
@@ -173,7 +165,7 @@ namespace OpenMM {
  * <tt><pre>
  * CustomIntegrator integrator(0.004);
  * integrator.addComputePerDof("v", "v+0.5*dt*f1/m");
- * for (int i = 0; i < 4; i++) {
+ * for (int i = 0; i &lt; 4; i++) {
  *     integrator.addComputePerDof("v", "v+0.5*(dt/4)*f0/m");
  *     integrator.addComputePerDof("x", "x+(dt/4)*v");
  *     integrator.addComputePerDof("v", "v+0.5*(dt/4)*f0/m");

platforms/opencl/src/OpenCLKernels.cpp

@@ -3787,12 +3787,11 @@ void OpenCLIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegr
         // Identify steps that can be merged into a single kernel.
         
         for (int step = 1; step < numSteps; step++) {
-            if ((needsForces[step] || needsEnergy[step]) && (invalidatesForces[step-1] || forceGroup[step] != forceGroup[step-1]))
+            if (needsForces[step] || needsEnergy[step])
                 continue;
             if (stepType[step-1] == CustomIntegrator::ComputeGlobal && stepType[step] == CustomIntegrator::ComputeGlobal)
                 merged[step] = true;
             if (stepType[step-1] == CustomIntegrator::ComputePerDof && stepType[step] == CustomIntegrator::ComputePerDof &&
-                    storePosAsDelta[step-1] == storePosAsDelta[step] && loadPosAsDelta[step-1] == loadPosAsDelta[step] &&
                     !usesVariable(expression[step], "uniform"))
                 merged[step] = true;
         }
@@ -3828,11 +3827,13 @@ void OpenCLIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegr
                     }
                     else if (variable[j] == "v")
                         compute << "velm[index] = " << convert << "velocity);\n";
-                    for (int i = 0; i < (int) perDofValues->getBuffers().size(); i++) {
-                        const OpenCLNonbondedUtilities::ParameterInfo& buffer = perDofValues->getBuffers()[i];
-                        compute << "perDofValues"<<intToString(i+1)<<"[3*index] = perDofx"<<intToString(i+1)<<";\n";
-                        compute << "perDofValues"<<intToString(i+1)<<"[3*index+1] = perDofy"<<intToString(i+1)<<";\n";
-                        compute << "perDofValues"<<intToString(i+1)<<"[3*index+2] = perDofz"<<intToString(i+1)<<";\n";
+                    else {
+                        for (int i = 0; i < (int) perDofValues->getBuffers().size(); i++) {
+                            const OpenCLNonbondedUtilities::ParameterInfo& buffer = perDofValues->getBuffers()[i];
+                            compute << "perDofValues"<<intToString(i+1)<<"[3*index] = perDofx"<<intToString(i+1)<<";\n";
+                            compute << "perDofValues"<<intToString(i+1)<<"[3*index+1] = perDofy"<<intToString(i+1)<<";\n";
+                            compute << "perDofValues"<<intToString(i+1)<<"[3*index+2] = perDofz"<<intToString(i+1)<<";\n";
+                        }
                     }
                     compute << "}\n";
                     numGaussian += numAtoms*usesVariable(expression[j], "gaussian");

platforms/opencl/tests/TestOpenCLCustomIntegrator.cpp

@@ -146,7 +146,7 @@ void testConstraints() {
         if (i == 1)
             initialEnergy = energy;
         else if (i > 1)
-            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
         integrator.step(1);
     }
 }
@@ -159,10 +159,12 @@ void testVelocityConstraints() {
     OpenCLPlatform platform;
     System system;
     CustomIntegrator integrator(0.002);
+    integrator.addPerDofVariable("x1", 0);
     integrator.addComputePerDof("v", "v+0.5*dt*f/m");
     integrator.addComputePerDof("x", "x+dt*v");
+    integrator.addComputePerDof("x1", "x");
     integrator.addConstrainPositions();
-    integrator.addComputePerDof("v", "v+0.5*dt*f/m");
+    integrator.addComputePerDof("v", "v+0.5*dt*f/m+(x-x1)/dt");
     integrator.addConstrainVelocities();
     integrator.setConstraintTolerance(1e-5);
     NonbondedForce* forceField = new NonbondedForce();
@@ -225,7 +227,7 @@ void testVelocityConstraints() {
         if (i == 0)
             initialEnergy = energy;
         else if (i > 0)
-            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
     }
 }
 

platforms/reference/tests/TestReferenceCustomIntegrator.cpp

@@ -146,7 +146,7 @@ void testConstraints() {
         if (i == 1)
             initialEnergy = energy;
         else if (i > 1)
-            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
         integrator.step(1);
     }
 }
@@ -159,10 +159,12 @@ void testVelocityConstraints() {
     ReferencePlatform platform;
     System system;
     CustomIntegrator integrator(0.002);
+    integrator.addPerDofVariable("x1", 0);
     integrator.addComputePerDof("v", "v+0.5*dt*f/m");
     integrator.addComputePerDof("x", "x+dt*v");
+    integrator.addComputePerDof("x1", "x");
     integrator.addConstrainPositions();
-    integrator.addComputePerDof("v", "v+0.5*dt*f/m");
+    integrator.addComputePerDof("v", "v+0.5*dt*f/m+(x-x1)/dt");
     integrator.addConstrainVelocities();
     integrator.setConstraintTolerance(1e-5);
     NonbondedForce* forceField = new NonbondedForce();
@@ -210,7 +212,7 @@ void testVelocityConstraints() {
         if (i == 1)
             initialEnergy = energy;
         else if (i > 1)
-            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
         integrator.step(2);
     }
 }