Implemented VariableVerletIntegrator::stepTo()

olla/include/openmm/kernels.h

@@ -441,8 +441,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
      */
-    virtual void execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator) = 0;
+    virtual void execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) = 0;
 };
 
 /**

openmmapi/include/openmm/Integrator.h

@@ -56,14 +56,14 @@ public:
     }
     /**
      * Get the size of each time step, in picoseconds.  If this integrator uses variable time steps,
-     * this is the initial step size that will be attempted for the next step.
+     * the size of the most recent step is returned.
      */
     double getStepSize() const {
         return stepSize;
     }
     /**
      * Set the size of each time step, in picoseconds.  If this integrator uses variable time steps,
-     * this is the initial step size that will be attempted for the next step.
+     * the effect of calling this method is undefined, and it may simply be ignored.
      */
     void setStepSize(double size) {
         stepSize = size;

openmmapi/include/openmm/VariableVerletIntegrator.h

@@ -84,6 +84,15 @@ public:
      * @param steps   the number of time steps to take
      */
     void step(int steps);
+    /**
+     * Advance a simulation through time by taking a series of steps until a specified time is
+     * reached.  When this method returns, the simulation time will exactly equal the time which
+     * was specified.  If you call this method and specify a time that is earlier than the
+     * current time, it will return without doing anything.
+     *
+     * @param time   the time to which the simulation should be advanced
+     */
+    void stepTo(double time);
 protected:
     /**
      * This will be called by the OpenMMContext when it is created.  It informs the Integrator

openmmapi/src/VariableVerletIntegrator.cpp

@@ -33,6 +33,7 @@
 #include "openmm/OpenMMContext.h"
 #include "openmm/internal/OpenMMContextImpl.h"
 #include "openmm/kernels.h"
+#include <limits>
 #include <string>
 
 using namespace OpenMM;
@@ -59,7 +60,14 @@ void VariableVerletIntegrator::step(int steps) {
     for (int i = 0; i < steps; ++i) {
         context->updateContextState();
         context->calcForces();
-        dynamic_cast<IntegrateVariableVerletStepKernel&>(kernel.getImpl()).execute(*context, *this);
+        dynamic_cast<IntegrateVariableVerletStepKernel&>(kernel.getImpl()).execute(*context, *this, std::numeric_limits<double>::infinity());
     }
 }
 
+void VariableVerletIntegrator::stepTo(double time) {
+    while (time > context->getTime()) {
+        context->updateContextState();
+        context->calcForces();
+        dynamic_cast<IntegrateVariableVerletStepKernel&>(kernel.getImpl()).execute(*context, *this, time);
+    }
+}

platforms/cuda/src/CudaKernels.cpp

@@ -562,7 +562,7 @@ void CudaIntegrateVariableVerletStepKernel::initialize(const System& system, con
     prevErrorTol = -1.0;
 }
 
-void CudaIntegrateVariableVerletStepKernel::execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator) {
+void CudaIntegrateVariableVerletStepKernel::execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
     _gpuContext* gpu = data.gpu;
     double errorTol = integrator.getErrorTolerance();
     if (errorTol != prevErrorTol) {
@@ -572,7 +572,8 @@ void CudaIntegrateVariableVerletStepKernel::execute(OpenMMContextImpl& context,
         gpuSetConstants(gpu);
         prevErrorTol = errorTol;
     }
-    kSelectVerletStepSize(gpu);
+    float maxStepSize = maxTime-data.time;
+    kSelectVerletStepSize(gpu, maxStepSize);
     kVerletUpdatePart1(gpu);
     kApplyFirstShake(gpu);
     kApplyFirstSettle(gpu);
@@ -583,6 +584,8 @@ void CudaIntegrateVariableVerletStepKernel::execute(OpenMMContextImpl& context,
     kVerletUpdatePart2(gpu);
     gpu->psStepSize->Download();
     data.time += (*gpu->psStepSize)[0].y;
+    if ((*gpu->psStepSize)[0].y == maxStepSize)
+        data.time = maxTime; // Avoid round-off error
     data.stepCount++;
 }
 

platforms/cuda/src/CudaKernels.h

@@ -393,8 +393,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
      */
-    void execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator);
+    void execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
 private:
     CudaPlatform::PlatformData& data;
     double prevErrorTol;

platforms/cuda/src/kernels/cudaKernels.h

@@ -50,7 +50,7 @@ extern void kLangevinUpdatePart1(gpuContext gpu);
 extern void kLangevinUpdatePart2(gpuContext gpu);
 extern void kVerletUpdatePart1(gpuContext gpu);
 extern void kVerletUpdatePart2(gpuContext gpu);
-extern void kSelectVerletStepSize(gpuContext gpu);
+extern void kSelectVerletStepSize(gpuContext gpu, float maxTimeStep);
 extern void kBrownianUpdatePart1(gpuContext gpu);
 extern void kBrownianUpdatePart2(gpuContext gpu);
 

platforms/cuda/src/kernels/kVerletUpdate.cu

@@ -89,7 +89,7 @@ void kVerletUpdatePart2(gpuContext gpu)
     }
 }
 
-__global__ void kSelectVerletStepSize_kernel()
+__global__ void kSelectVerletStepSize_kernel(float maxStepSize)
 {
     // Calculate the error.
 
@@ -121,14 +121,16 @@ __global__ void kSelectVerletStepSize_kernel()
             newStepSize = min(newStepSize, oldStepSize*2.0f); // For safety, limit how quickly dt can increase.
         if (newStepSize > oldStepSize && newStepSize < 1.1f*oldStepSize)
             newStepSize = oldStepSize; // Keeping dt constant between steps improves the behavior of the integrator.
+        if (newStepSize > maxStepSize)
+            newStepSize = maxStepSize;
         cSim.pStepSize[0].y = newStepSize;
     }
 }
 
-void kSelectVerletStepSize(gpuContext gpu)
+void kSelectVerletStepSize(gpuContext gpu, float maxTimeStep)
 {
 //    printf("kSelectVerletStepSize\n");
-    kSelectVerletStepSize_kernel<<<1, gpu->sim.update_threads_per_block, sizeof(float)*gpu->sim.update_threads_per_block>>>();
+    kSelectVerletStepSize_kernel<<<1, gpu->sim.update_threads_per_block, sizeof(float)*gpu->sim.update_threads_per_block>>>(maxTimeStep);
     LAUNCHERROR("kSelectVerletStepSize");
 }
 

platforms/cuda/tests/TestCudaVariableVerletIntegrator.cpp

@@ -136,6 +136,14 @@ void testConstraints() {
             ASSERT_EQUAL_TOL(initialEnergy, energy, 0.1);
         integrator.step(1);
     }
+    double finalTime = context.getState(State::Positions).getTime();
+    ASSERT(finalTime > 0.1);
+
+    // Now try the stepTo() method.
+
+    finalTime += 0.5;
+    integrator.stepTo(finalTime);
+    ASSERT_EQUAL(finalTime, context.getState(State::Positions).getTime());
 }
 
 void testConstrainedClusters() {
@@ -198,6 +206,7 @@ void testConstrainedClusters() {
             ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
         integrator.step(1);
     }
+    ASSERT(context.getState(State::Positions).getTime() > 0.1);
 }
 
 int main() {

platforms/reference/include/ReferencePlatform.h

@@ -62,8 +62,9 @@ private:
 
 class ReferencePlatform::PlatformData {
 public:
-    PlatformData() : time(0.0) {
+    PlatformData() : time(0.0), stepCount(0) {
     }
+    int stepCount;
     double time;
 };
 } // namespace OpenMM

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -71,6 +71,6 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
     if (name == CalcKineticEnergyKernel::Name())
         return new ReferenceCalcKineticEnergyKernel(name, platform);
     if (name == RemoveCMMotionKernel::Name())
-        return new ReferenceRemoveCMMotionKernel(name, platform);
+        return new ReferenceRemoveCMMotionKernel(name, platform, data);
     throw OpenMMException( (std::string("Tried to create kernel with illegal kernel name '") + name + "'").c_str() );
 }

platforms/reference/src/ReferenceKernels.cpp

@@ -655,6 +655,7 @@ void ReferenceIntegrateVerletStepKernel::execute(OpenMMContextImpl& context, con
     }
     dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
     data.time += stepSize;
+    data.stepCount++;
 }
 
 ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
@@ -719,6 +720,7 @@ void ReferenceIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, c
     }
     dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
     data.time += stepSize;
+    data.stepCount++;
 }
 
 ReferenceIntegrateBrownianStepKernel::~ReferenceIntegrateBrownianStepKernel() {
@@ -782,6 +784,7 @@ void ReferenceIntegrateBrownianStepKernel::execute(OpenMMContextImpl& context, c
     }
     dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
     data.time += stepSize;
+    data.stepCount++;
 }
 
 ReferenceIntegrateVariableVerletStepKernel::~ReferenceIntegrateVariableVerletStepKernel() {
@@ -815,7 +818,7 @@ void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system
     }
 }
 
-void ReferenceIntegrateVariableVerletStepKernel::execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator) {
+void ReferenceIntegrateVariableVerletStepKernel::execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
     double errorTol = integrator.getErrorTolerance();
     RealOpenMM** posData = ((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData();
     RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
@@ -834,8 +837,12 @@ void ReferenceIntegrateVariableVerletStepKernel::execute(OpenMMContextImpl& cont
         dynamics->setReferenceConstraintAlgorithm(constraints);
         prevErrorTol = errorTol;
     }
-    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
+    RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
+    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, maxStepSize);
     data.time += dynamics->getDeltaT();
+    if (dynamics->getDeltaT() == maxStepSize)
+        data.time = maxTime; // Avoid round-off error
+    data.stepCount++;
 }
 
 ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
@@ -888,8 +895,7 @@ void ReferenceRemoveCMMotionKernel::initialize(const System& system, const CMMot
 }
 
 void ReferenceRemoveCMMotionKernel::execute(OpenMMContextImpl& context) {
-    int step = (int)std::floor(context.getTime()/context.getIntegrator().getStepSize());
-    if (step%frequency != 0)
+    if (data.stepCount%frequency != 0)
         return;
     RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
     

platforms/reference/src/ReferenceKernels.h

@@ -463,8 +463,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
      */
-    void execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator);
+    void execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
 private:
     ReferencePlatform::PlatformData& data;
     ReferenceVariableVerletDynamics* dynamics;
@@ -530,7 +531,7 @@ private:
  */
 class ReferenceRemoveCMMotionKernel : public RemoveCMMotionKernel {
 public:
-    ReferenceRemoveCMMotionKernel(std::string name, const Platform& platform) : RemoveCMMotionKernel(name, platform) {
+    ReferenceRemoveCMMotionKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : RemoveCMMotionKernel(name, platform), data(data) {
     }
     /**
      * Initialize the kernel, setting up the particle masses.
@@ -546,6 +547,7 @@ public:
      */
     void execute(OpenMMContextImpl& context);
 private:
+    ReferencePlatform::PlatformData& data;
     std::vector<double> masses;
     int frequency;
 };

platforms/reference/src/SimTKReference/ReferenceVariableVerletDynamics.cpp

@@ -131,6 +131,7 @@ int ReferenceVariableVerletDynamics::printParameters( std::stringstream& message
    @param velocities          velocities
    @param forces              forces
    @param masses              atom masses
+   @param maxStepSize         maximum time step
 
    @return ReferenceDynamics::DefaultReturn
 
@@ -138,7 +139,7 @@ int ReferenceVariableVerletDynamics::printParameters( std::stringstream& message
 
 int ReferenceVariableVerletDynamics::update( int numberOfAtoms, RealOpenMM** atomCoordinates,
                                           RealOpenMM** velocities,
-                                          RealOpenMM** forces, RealOpenMM* masses ){
+                                          RealOpenMM** forces, RealOpenMM* masses, RealOpenMM maxStepSize ){
 
     // ---------------------------------------------------------------------------------------
 
@@ -195,6 +196,8 @@ int ReferenceVariableVerletDynamics::update( int numberOfAtoms, RealOpenMM** ato
         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.
+    if (newStepSize > maxStepSize)
+        newStepSize = maxStepSize;
     RealOpenMM vstep = 0.5f*(newStepSize+getDeltaT()); // The time interval by which to advance the velocities
     setDeltaT(newStepSize);
     for (int i = 0; i < numberOfAtoms; ++i) {

platforms/reference/src/SimTKReference/ReferenceVariableVerletDynamics.h

@@ -97,13 +97,14 @@ class ReferenceVariableVerletDynamics : public ReferenceDynamics {
          @param velocities          velocities
          @param forces              forces
          @param masses              atom masses
+         @param maxStepSize         maximum time step
 
          @return ReferenceDynamics::DefaultReturn
 
          --------------------------------------------------------------------------------------- */
 
       int update( int numberOfAtoms, RealOpenMM** atomCoordinates,
-                  RealOpenMM** velocities, RealOpenMM** forces, RealOpenMM* masses );
+                  RealOpenMM** velocities, RealOpenMM** forces, RealOpenMM* masses, RealOpenMM maxStepSize );
 
 };
 

platforms/reference/tests/TestReferenceVariableVerletIntegrator.cpp

@@ -133,7 +133,14 @@ void testConstraints() {
             ASSERT_EQUAL_TOL(initialEnergy, energy, 0.1);
         integrator.step(1);
     }
-    ASSERT(context.getState(State::Positions).getTime() > 0.1);
+    double finalTime = context.getState(State::Positions).getTime();
+    ASSERT(finalTime > 0.1);
+
+    // Now try the stepTo() method.
+
+    finalTime += 0.5;
+    integrator.stepTo(finalTime);
+    ASSERT_EQUAL(finalTime, context.getState(State::Positions).getTime());
 }
 
 void testConstrainedClusters() {