Variable time step integrators correctly report the size of the last step...

Variable time step integrators correctly report the size of the last step taken.  Also improved accuracy of integration unit tests.

olla/include/openmm/kernels.h

@@ -780,8 +780,9 @@ public:
      * @param context    the context in which to execute this kernel
      * @param integrator the LangevinIntegrator this kernel is being used for
      * @param maxTime    the maximum time beyond which the simulation should not be advanced
+     * @return the size of the step that was taken
      */
-    virtual void execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) = 0;
+    virtual double execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) = 0;
 };
 
 /**
@@ -807,8 +808,9 @@ public:
      * @param context    the context in which to execute this kernel
      * @param integrator the VerletIntegrator this kernel is being used for
      * @param maxTime    the maximum time beyond which the simulation should not be advanced
+     * @return the size of the step that was taken
      */
-    virtual void execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) = 0;
+    virtual double execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) = 0;
 };
 
 /**

openmmapi/src/VariableLangevinIntegrator.cpp

@@ -69,7 +69,7 @@ void VariableLangevinIntegrator::step(int steps) {
     for (int i = 0; i < steps; ++i) {
         context->updateContextState();
         context->calcForcesAndEnergy(true, false);
-        dynamic_cast<IntegrateVariableLangevinStepKernel&>(kernel.getImpl()).execute(*context, *this, std::numeric_limits<double>::infinity());
+        setStepSize(dynamic_cast<IntegrateVariableLangevinStepKernel&>(kernel.getImpl()).execute(*context, *this, std::numeric_limits<double>::infinity()));
     }
 }
 
@@ -77,6 +77,6 @@ void VariableLangevinIntegrator::stepTo(double time) {
     while (time > context->getTime()) {
         context->updateContextState();
         context->calcForcesAndEnergy(true, false);
-        dynamic_cast<IntegrateVariableLangevinStepKernel&>(kernel.getImpl()).execute(*context, *this, time);
+        setStepSize(dynamic_cast<IntegrateVariableLangevinStepKernel&>(kernel.getImpl()).execute(*context, *this, time));
     }
 }

openmmapi/src/VariableVerletIntegrator.cpp

@@ -64,7 +64,7 @@ void VariableVerletIntegrator::step(int steps) {
     for (int i = 0; i < steps; ++i) {
         context->updateContextState();
         context->calcForcesAndEnergy(true, false);
-        dynamic_cast<IntegrateVariableVerletStepKernel&>(kernel.getImpl()).execute(*context, *this, std::numeric_limits<double>::infinity());
+        setStepSize(dynamic_cast<IntegrateVariableVerletStepKernel&>(kernel.getImpl()).execute(*context, *this, std::numeric_limits<double>::infinity()));
     }
 }
 
@@ -72,6 +72,6 @@ void VariableVerletIntegrator::stepTo(double time) {
     while (time > context->getTime()) {
         context->updateContextState();
         context->calcForcesAndEnergy(true, false);
-        dynamic_cast<IntegrateVariableVerletStepKernel&>(kernel.getImpl()).execute(*context, *this, time);
+        setStepSize(dynamic_cast<IntegrateVariableVerletStepKernel&>(kernel.getImpl()).execute(*context, *this, time));
     }
 }

platforms/cuda/src/CudaKernels.cpp

@@ -1343,7 +1343,7 @@ void CudaIntegrateVariableVerletStepKernel::initialize(const System& system, con
     prevErrorTol = -1.0;
 }
 
-void CudaIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
+double CudaIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
     _gpuContext* gpu = data.gpu;
     double errorTol = integrator.getErrorTolerance();
     if (errorTol != prevErrorTol) {
@@ -1368,6 +1368,7 @@ void CudaIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const
     if ((*gpu->psStepSize)[0].y == maxStepSize)
         data.time = maxTime; // Avoid round-off error
     data.stepCount++;
+    return (*gpu->psStepSize)[0].y;
 }
 
 CudaIntegrateVariableLangevinStepKernel::~CudaIntegrateVariableLangevinStepKernel() {
@@ -1383,7 +1384,7 @@ void CudaIntegrateVariableLangevinStepKernel::initialize(const System& system, c
     prevErrorTol = -1.0;
 }
 
-void CudaIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
+double CudaIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
     _gpuContext* gpu = data.gpu;
     double temperature = integrator.getTemperature();
     double friction = integrator.getFriction();
@@ -1415,6 +1416,7 @@ void CudaIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, cons
     if ((*gpu->psStepSize)[0].y == maxStepSize)
         data.time = maxTime; // Avoid round-off error
     data.stepCount++;
+    return (*gpu->psStepSize)[0].y;
 }
 
 CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {

platforms/cuda/src/CudaKernels.h

@@ -734,8 +734,9 @@ public:
      * @param context    the context in which to execute this kernel
      * @param integrator the VerletIntegrator this kernel is being used for
      * @param maxTime    the maximum time beyond which the simulation should not be advanced
+     * @return the size of the step that was taken
      */
-    void execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
+    double execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
 private:
     CudaPlatform::PlatformData& data;
     double prevErrorTol;
@@ -762,8 +763,9 @@ public:
      * @param context    the context in which to execute this kernel
      * @param integrator the VariableLangevinIntegrator this kernel is being used for
      * @param maxTime    the maximum time beyond which the simulation should not be advanced
+     * @return the size of the step that was taken
      */
-    void execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
+    double execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
 private:
     CudaPlatform::PlatformData& data;
     double prevTemp, prevFriction, prevErrorTol;

platforms/cuda/tests/TestCudaVariableVerletIntegrator.cpp

@@ -50,6 +50,23 @@ using namespace std;
 
 const double TOL = 1e-5;
 
+/**
+ * Compute the energy of a state, taking into account the half step offset between
+ * positions and velocities.
+ */
+
+static double computeEnergy(const State& state, const System& system, double dt) {
+    const vector<Vec3>& v = state.getVelocities();
+    const vector<Vec3>& f = state.getForces();
+    double energy = 0.0;
+    for (int i = 0; i < system.getNumParticles(); i++) {
+        double m = system.getParticleMass(i);
+        Vec3 vel = v[i]+f[i]*(0.5*dt/m);
+        energy += 0.5*m*vel.dot(vel);
+    }
+    return energy+state.getPotentialEnergy();
+}
+
 void testSingleBond() {
     CudaPlatform platform;
     System system;
@@ -121,7 +138,7 @@ void testConstraints() {
 
     double initialEnergy = 0.0;
     for (int i = 0; i < 1000; ++i) {
-        State state = context.getState(State::Positions | State::Energy);
+        State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
         for (int j = 0; j < numConstraints; ++j) {
             int particle1, particle2;
             double distance;
@@ -131,11 +148,11 @@ void testConstraints() {
             double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
             ASSERT_EQUAL_TOL(distance, dist, 1e-4);
         }
-        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        double energy = computeEnergy(state, system, integrator.getStepSize());
         if (i == 1)
             initialEnergy = energy;
         else if (i > 1)
-            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.1);
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
         integrator.step(1);
     }
     double finalTime = context.getState(State::Positions).getTime();
@@ -193,7 +210,7 @@ void testConstrainedClusters() {
 
     double initialEnergy = 0.0;
     for (int i = 0; i < 1000; ++i) {
-        State state = context.getState(State::Positions | State::Energy);
+        State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
         for (int j = 0; j < system.getNumConstraints(); ++j) {
             int particle1, particle2;
             double distance;
@@ -203,11 +220,11 @@ void testConstrainedClusters() {
             double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
             ASSERT_EQUAL_TOL(distance, dist, 2e-5);
         }
-        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        double energy = computeEnergy(state, system, integrator.getStepSize());
         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);
     }
     ASSERT(context.getState(State::Positions).getTime() > 0.1);

platforms/cuda/tests/TestCudaVerletIntegrator.cpp

@@ -50,6 +50,23 @@ using namespace std;
 
 const double TOL = 1e-5;
 
+/**
+ * Compute the energy of a state, taking into account the half step offset between
+ * positions and velocities.
+ */
+
+static double computeEnergy(const State& state, const System& system, double dt) {
+    const vector<Vec3>& v = state.getVelocities();
+    const vector<Vec3>& f = state.getForces();
+    double energy = 0.0;
+    for (int i = 0; i < system.getNumParticles(); i++) {
+        double m = system.getParticleMass(i);
+        Vec3 vel = v[i]+f[i]*(0.5*dt/m);
+        energy += 0.5*m*vel.dot(vel);
+    }
+    return energy+state.getPotentialEnergy();
+}
+
 void testSingleBond() {
     CudaPlatform platform;
     System system;
@@ -121,7 +138,7 @@ void testConstraints() {
     
     double initialEnergy = 0.0;
     for (int i = 0; i < 1000; ++i) {
-        State state = context.getState(State::Positions | State::Energy);
+        State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
             for (int j = 0; j < numConstraints; ++j) {
                 int particle1, particle2;
                 double distance;
@@ -131,11 +148,11 @@ void testConstraints() {
                 double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
                 ASSERT_EQUAL_TOL(distance, dist, 1e-4);
             }
-        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        double energy = computeEnergy(state, system, integrator.getStepSize());
         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);
     }
 }
@@ -145,7 +162,7 @@ void testConstrainedClusters() {
     const double temp = 500.0;
     CudaPlatform platform;
     System system;
-    VerletIntegrator integrator(0.002);
+    VerletIntegrator integrator(0.001);
     integrator.setConstraintTolerance(1e-5);
     NonbondedForce* forceField = new NonbondedForce();
     for (int i = 0; i < numParticles; ++i) {
@@ -185,7 +202,7 @@ void testConstrainedClusters() {
 
     double initialEnergy = 0.0;
     for (int i = 0; i < 1000; ++i) {
-        State state = context.getState(State::Positions | State::Energy);
+        State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
         for (int j = 0; j < system.getNumConstraints(); ++j) {
             int particle1, particle2;
             double distance;
@@ -195,11 +212,11 @@ void testConstrainedClusters() {
             double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
             ASSERT_EQUAL_TOL(distance, dist, 2e-5);
         }
-        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        double energy = computeEnergy(state, system, integrator.getStepSize());
         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);
     }
 }

platforms/opencl/src/OpenCLKernels.cpp

@@ -3400,7 +3400,7 @@ void OpenCLIntegrateVariableVerletStepKernel::initialize(const System& system, c
     blockSize = std::min(std::min(256, system.getNumParticles()), (int) cl.getDevice().getInfo<CL_DEVICE_MAX_WORK_GROUP_SIZE>());
 }
 
-void OpenCLIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
+double OpenCLIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
     OpenCLIntegrationUtilities& integration = cl.getIntegrationUtilities();
     int numAtoms = cl.getNumAtoms();
     if (!hasInitializedKernels) {
@@ -3452,6 +3452,7 @@ void OpenCLIntegrateVariableVerletStepKernel::execute(ContextImpl& context, cons
         time = maxTime; // Avoid round-off error
     cl.setTime(time);
     cl.setStepCount(cl.getStepCount()+1);
+    return dt;
 }
 
 OpenCLIntegrateVariableLangevinStepKernel::~OpenCLIntegrateVariableLangevinStepKernel() {
@@ -3475,7 +3476,7 @@ void OpenCLIntegrateVariableLangevinStepKernel::initialize(const System& system,
     blockSize = std::min(blockSize, (int) cl.getDevice().getInfo<CL_DEVICE_MAX_WORK_GROUP_SIZE>());
 }
 
-void OpenCLIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
+double OpenCLIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
     OpenCLIntegrationUtilities& integration = cl.getIntegrationUtilities();
     int numAtoms = cl.getNumAtoms();
     if (!hasInitializedKernels) {
@@ -3530,6 +3531,7 @@ void OpenCLIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, co
         time = maxTime; // Avoid round-off error
     cl.setTime(time);
     cl.setStepCount(cl.getStepCount()+1);
+    return dt;
 }
 
 class OpenCLIntegrateCustomStepKernel::ReorderListener : public OpenCLContext::ReorderListener {

platforms/opencl/src/OpenCLKernels.h

@@ -875,8 +875,9 @@ public:
      * @param context    the context in which to execute this kernel
      * @param integrator the VerletIntegrator this kernel is being used for
      * @param maxTime    the maximum time beyond which the simulation should not be advanced
+     * @return the size of the step that was taken
      */
-    void execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
+    double execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
 private:
     OpenCLContext& cl;
     bool hasInitializedKernels;
@@ -906,8 +907,9 @@ public:
      * @param context    the context in which to execute this kernel
      * @param integrator the VariableLangevinIntegrator this kernel is being used for
      * @param maxTime    the maximum time beyond which the simulation should not be advanced
+     * @return the size of the step that was taken
      */
-    void execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
+    double execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
 private:
     OpenCLContext& cl;
     bool hasInitializedKernels;

platforms/opencl/tests/TestOpenCLVariableVerletIntegrator.cpp

@@ -50,6 +50,23 @@ using namespace std;
 
 const double TOL = 1e-5;
 
+/**
+ * Compute the energy of a state, taking into account the half step offset between
+ * positions and velocities.
+ */
+
+static double computeEnergy(const State& state, const System& system, double dt) {
+    const vector<Vec3>& v = state.getVelocities();
+    const vector<Vec3>& f = state.getForces();
+    double energy = 0.0;
+    for (int i = 0; i < system.getNumParticles(); i++) {
+        double m = system.getParticleMass(i);
+        Vec3 vel = v[i]+f[i]*(0.5*dt/m);
+        energy += 0.5*m*vel.dot(vel);
+    }
+    return energy+state.getPotentialEnergy();
+}
+
 void testSingleBond() {
     OpenCLPlatform platform;
     System system;
@@ -121,7 +138,7 @@ void testConstraints() {
 
     double initialEnergy = 0.0;
     for (int i = 0; i < 1000; ++i) {
-        State state = context.getState(State::Positions | State::Energy);
+        State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
             for (int j = 0; j < numConstraints; ++j) {
                 int particle1, particle2;
                 double distance;
@@ -131,11 +148,11 @@ void testConstraints() {
                 double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
                 ASSERT_EQUAL_TOL(distance, dist, 1e-4);
             }
-        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        double energy = computeEnergy(state, system, integrator.getStepSize());
         if (i == 1)
             initialEnergy = energy;
         else if (i > 1)
-            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.1);
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
         integrator.step(1);
     }
     double finalTime = context.getState(State::Positions).getTime();
@@ -193,7 +210,7 @@ void testConstrainedClusters() {
 
     double initialEnergy = 0.0;
     for (int i = 0; i < 1000; ++i) {
-        State state = context.getState(State::Positions | State::Energy);
+        State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
         for (int j = 0; j < system.getNumConstraints(); ++j) {
             int particle1, particle2;
             double distance;
@@ -203,11 +220,11 @@ void testConstrainedClusters() {
             double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
             ASSERT_EQUAL_TOL(distance, dist, 2e-5);
         }
-        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        double energy = computeEnergy(state, system, integrator.getStepSize());
         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);
     }
     ASSERT(context.getState(State::Positions).getTime() > 0.1);

platforms/opencl/tests/TestOpenCLVerletIntegrator.cpp

@@ -50,6 +50,23 @@ using namespace std;
 
 const double TOL = 1e-5;
 
+/**
+ * Compute the energy of a state, taking into account the half step offset between
+ * positions and velocities.
+ */
+
+static double computeEnergy(const State& state, const System& system, double dt) {
+    const vector<Vec3>& v = state.getVelocities();
+    const vector<Vec3>& f = state.getForces();
+    double energy = 0.0;
+    for (int i = 0; i < system.getNumParticles(); i++) {
+        double m = system.getParticleMass(i);
+        Vec3 vel = v[i]+f[i]*(0.5*dt/m);
+        energy += 0.5*m*vel.dot(vel);
+    }
+    return energy+state.getPotentialEnergy();
+}
+
 void testSingleBond() {
     OpenCLPlatform platform;
     System system;
@@ -121,7 +138,7 @@ void testConstraints() {
 
     double initialEnergy = 0.0;
     for (int i = 0; i < 1000; ++i) {
-        State state = context.getState(State::Positions | State::Energy);
+        State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
         for (int j = 0; j < numConstraints; ++j) {
             int particle1, particle2;
             double distance;
@@ -131,11 +148,11 @@ void testConstraints() {
             double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
             ASSERT_EQUAL_TOL(distance, dist, 1e-4);
         }
-        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        double energy = computeEnergy(state, system, integrator.getStepSize());
         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);
     }
 }
@@ -145,7 +162,7 @@ void testConstrainedClusters() {
     const double temp = 500.0;
     OpenCLPlatform platform;
     System system;
-    VerletIntegrator integrator(0.002);
+    VerletIntegrator integrator(0.001);
     integrator.setConstraintTolerance(1e-5);
     NonbondedForce* forceField = new NonbondedForce();
     for (int i = 0; i < numParticles; ++i) {
@@ -185,7 +202,7 @@ void testConstrainedClusters() {
 
     double initialEnergy = 0.0;
     for (int i = 0; i < 1000; ++i) {
-        State state = context.getState(State::Positions | State::Energy);
+        State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
         for (int j = 0; j < system.getNumConstraints(); ++j) {
             int particle1, particle2;
             double distance;
@@ -195,11 +212,11 @@ void testConstrainedClusters() {
             double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
             ASSERT_EQUAL_TOL(distance, dist, 2e-5);
         }
-        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        double energy = computeEnergy(state, system, integrator.getStepSize());
         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);
     }
 }

platforms/reference/src/ReferenceKernels.cpp

@@ -1447,7 +1447,7 @@ void ReferenceIntegrateVariableLangevinStepKernel::initialize(const System& syst
     SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
 }
 
-void ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
+double ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
     double temperature = integrator.getTemperature();
     double friction = integrator.getFriction();
     double errorTol = integrator.getErrorTolerance();
@@ -1478,6 +1478,7 @@ void ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context,
     if (dynamics->getDeltaT() == maxStepSize)
         data.time = maxTime; // Avoid round-off error
     data.stepCount++;
+    return dynamics->getDeltaT();
 }
 
 ReferenceIntegrateVariableVerletStepKernel::~ReferenceIntegrateVariableVerletStepKernel() {
@@ -1509,7 +1510,7 @@ void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system
     }
 }
 
-void ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
+double ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
     double errorTol = integrator.getErrorTolerance();
     vector<RealVec>& posData = extractPositions(context);
     vector<RealVec>& velData = extractVelocities(context);
@@ -1535,6 +1536,7 @@ void ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, c
     if (dynamics->getDeltaT() == maxStepSize)
         data.time = maxTime; // Avoid round-off error
     data.stepCount++;
+    return dynamics->getDeltaT();
 }
 
 ReferenceIntegrateCustomStepKernel::~ReferenceIntegrateCustomStepKernel() {

platforms/reference/src/ReferenceKernels.h

@@ -836,8 +836,9 @@ public:
      * @param context    the context in which to execute this kernel
      * @param integrator the LangevinIntegrator this kernel is being used for
      * @param maxTime    the maximum time beyond which the simulation should not be advanced
+     * @return the size of the step that was taken
      */
-    void execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
+    double execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
 private:
     ReferencePlatform::PlatformData& data;
     ReferenceVariableStochasticDynamics* dynamics;
@@ -871,8 +872,9 @@ public:
      * @param context    the context in which to execute this kernel
      * @param integrator the VerletIntegrator this kernel is being used for
      * @param maxTime    the maximum time beyond which the simulation should not be advanced
+     * @return the size of the step that was taken
      */
-    void execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
+    double execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
 private:
     ReferencePlatform::PlatformData& data;
     ReferenceVariableVerletDynamics* dynamics;

platforms/reference/tests/TestReferenceVariableVerletIntegrator.cpp

@@ -50,6 +50,23 @@ using namespace std;
 
 const double TOL = 1e-5;
 
+/**
+ * Compute the energy of a state, taking into account the half step offset between
+ * positions and velocities.
+ */
+
+static double computeEnergy(const State& state, const System& system, double dt) {
+    const vector<Vec3>& v = state.getVelocities();
+    const vector<Vec3>& f = state.getForces();
+    double energy = 0.0;
+    for (int i = 0; i < system.getNumParticles(); i++) {
+        double m = system.getParticleMass(i);
+        Vec3 vel = v[i]+f[i]*(0.5*dt/m);
+        energy += 0.5*m*vel.dot(vel);
+    }
+    return energy+state.getPotentialEnergy();
+}
+
 void testSingleBond() {
     ReferencePlatform platform;
     System system;
@@ -118,7 +135,7 @@ void testConstraints() {
 
     double initialEnergy = 0.0;
     for (int i = 0; i < 1000; ++i) {
-        State state = context.getState(State::Positions | State::Energy);
+        State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
         for (int j = 0; j < system.getNumConstraints(); ++j) {
             int particle1, particle2;
             double distance;
@@ -128,11 +145,11 @@ void testConstraints() {
             double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
             ASSERT_EQUAL_TOL(distance, dist, 2e-5);
         }
-        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        double energy = computeEnergy(state, system, integrator.getStepSize());
         if (i == 1)
             initialEnergy = energy;
         else if (i > 1)
-            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.1);
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
         integrator.step(1);
     }
     double finalTime = context.getState(State::Positions).getTime();
@@ -190,7 +207,7 @@ void testConstrainedClusters() {
 
     double initialEnergy = 0.0;
     for (int i = 0; i < 1000; ++i) {
-        State state = context.getState(State::Positions | State::Energy);
+        State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
         for (int j = 0; j < system.getNumConstraints(); ++j) {
             int particle1, particle2;
             double distance;
@@ -200,11 +217,11 @@ void testConstrainedClusters() {
             double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
             ASSERT_EQUAL_TOL(distance, dist, 2e-5);
         }
-        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        double energy = computeEnergy(state, system, integrator.getStepSize());
         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);
     }
     ASSERT(context.getState(State::Positions).getTime() > 0.1);

platforms/reference/tests/TestReferenceVerletIntegrator.cpp

@@ -50,6 +50,23 @@ using namespace std;
 
 const double TOL = 1e-5;
 
+/**
+ * Compute the energy of a state, taking into account the half step offset between
+ * positions and velocities.
+ */
+
+static double computeEnergy(const State& state, const System& system, double dt) {
+    const vector<Vec3>& v = state.getVelocities();
+    const vector<Vec3>& f = state.getForces();
+    double energy = 0.0;
+    for (int i = 0; i < system.getNumParticles(); i++) {
+        double m = system.getParticleMass(i);
+        Vec3 vel = v[i]+f[i]*(0.5*dt/m);
+        energy += 0.5*m*vel.dot(vel);
+    }
+    return energy+state.getPotentialEnergy();
+}
+
 void testSingleBond() {
     ReferencePlatform platform;
     System system;
@@ -117,7 +134,7 @@ void testConstraints() {
     
     double initialEnergy = 0.0;
     for (int i = 0; i < 1000; ++i) {
-        State state = context.getState(State::Positions | State::Energy);
+        State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
         for (int j = 0; j < system.getNumConstraints(); ++j) {
             int particle1, particle2;
             double distance;
@@ -127,11 +144,11 @@ void testConstraints() {
             double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
             ASSERT_EQUAL_TOL(distance, dist, 2e-5);
         }
-        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        double energy = computeEnergy(state, system, integrator.getStepSize());
         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);
     }
 }
@@ -141,7 +158,7 @@ void testConstrainedClusters() {
     const double temp = 500.0;
     ReferencePlatform platform;
     System system;
-    VerletIntegrator integrator(0.002);
+    VerletIntegrator integrator(0.001);
     integrator.setConstraintTolerance(1e-5);
     NonbondedForce* forceField = new NonbondedForce();
     for (int i = 0; i < numParticles; ++i) {
@@ -181,7 +198,7 @@ void testConstrainedClusters() {
 
     double initialEnergy = 0.0;
     for (int i = 0; i < 1000; ++i) {
-        State state = context.getState(State::Positions | State::Energy);
+        State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
         for (int j = 0; j < system.getNumConstraints(); ++j) {
             int particle1, particle2;
             double distance;
@@ -191,11 +208,11 @@ void testConstrainedClusters() {
             double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
             ASSERT_EQUAL_TOL(distance, dist, 2e-5);
         }
-        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        double energy = computeEnergy(state, system, integrator.getStepSize());
         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);
     }
 }