Converted from BAOAB to leapfrog

olla/include/openmm/kernels.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-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -1150,12 +1150,8 @@ public:
      * 
      * @param context    the context in which to execute this kernel
      * @param integrator the BAOABLangevinIntegrator this kernel is being used for
-     * @param forcesAreValid if the context has been modified since the last time step, this will be
-     *                       false to show that cached forces are invalid and must be recalculated.
-     *                       On exit, this should specify whether the cached forces are valid at the
-     *                       end of the step.
      */
-    virtual void execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid) = 0;
+    virtual void execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator) = 0;
     /**
      * Compute the kinetic energy.
      * 

openmmapi/include/openmm/BAOABLangevinIntegrator.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-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -128,10 +128,6 @@ protected:
      * cleanup.  It will also get called again if the application calls reinitialize() on the Context.
      */
     void cleanup();
-    /**
-     * When the user modifies the state, we need to mark that the forces need to be recalculated.
-     */
-    void stateChanged(State::DataType changed);
     /**
      * Get the names of all Kernels used by this Integrator.
      */
@@ -147,7 +143,6 @@ protected:
 private:
     double temperature, friction;
     int randomNumberSeed;
-    bool forcesAreValid;
     Kernel kernel;
 };
 

openmmapi/src/BAOABLangevinIntegrator.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-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -46,7 +46,6 @@ BAOABLangevinIntegrator::BAOABLangevinIntegrator(double temperature, double fric
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
     setRandomNumberSeed(0);
-    forcesAreValid = false;
 }
 
 void BAOABLangevinIntegrator::initialize(ContextImpl& contextRef) {
@@ -62,10 +61,6 @@ void BAOABLangevinIntegrator::cleanup() {
     kernel = Kernel();
 }
 
-void BAOABLangevinIntegrator::stateChanged(State::DataType changed) {
-    forcesAreValid = false;
-}
-
 vector<string> BAOABLangevinIntegrator::getKernelNames() {
     std::vector<std::string> names;
     names.push_back(IntegrateBAOABStepKernel::Name());
@@ -73,7 +68,6 @@ vector<string> BAOABLangevinIntegrator::getKernelNames() {
 }
 
 double BAOABLangevinIntegrator::computeKineticEnergy() {
-    forcesAreValid = false;
     return kernel.getAs<IntegrateBAOABStepKernel>().computeKineticEnergy(*context, *this);
 }
 
@@ -85,8 +79,8 @@ void BAOABLangevinIntegrator::step(int steps) {
     if (context == NULL)
         throw OpenMMException("This Integrator is not bound to a context!");  
     for (int i = 0; i < steps; ++i) {
-        if (context->updateContextState())
-            forcesAreValid = false;
-        kernel.getAs<IntegrateBAOABStepKernel>().execute(*context, *this, forcesAreValid);
+        context->updateContextState();
+        context->calcForcesAndEnergy(true, false);
+        kernel.getAs<IntegrateBAOABStepKernel>().execute(*context, *this);
     }
 }

platforms/common/include/openmm/common/CommonKernels.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-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -924,12 +924,8 @@ public:
      * 
      * @param context    the context in which to execute this kernel
      * @param integrator the BAOABLangevinIntegrator this kernel is being used for
-     * @param forcesAreValid if the context has been modified since the last time step, this will be
-     *                       false to show that cached forces are invalid and must be recalculated.
-     *                       On exit, this should specify whether the cached forces are valid at the
-     *                       end of the step.
      */
-    void execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid);
+    void execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator);
     /**
      * Compute the kinetic energy.
      * 
@@ -942,7 +938,7 @@ private:
     double prevTemp, prevFriction, prevStepSize;
     bool hasInitializedKernels;
     ComputeArray params, oldDelta;
-    ComputeKernel kernel1, kernel2, kernel3, kernel4;
+    ComputeKernel kernel1, kernel2, kernel3;
 };
 
 /**

platforms/common/src/CommonKernels.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-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -5555,7 +5555,6 @@ void CommonIntegrateBAOABStepKernel::initialize(const System& system, const BAOA
     kernel1 = program->createKernel("integrateBAOABPart1");
     kernel2 = program->createKernel("integrateBAOABPart2");
     kernel3 = program->createKernel("integrateBAOABPart3");
-    kernel4 = program->createKernel("integrateBAOABPart4");
     if (cc.getUseDoublePrecision() || cc.getUseMixedPrecision()) {
         params.initialize<double>(cc, 2, "baoabParams");
         oldDelta.initialize<mm_double4>(cc, cc.getPaddedNumAtoms(), "oldDelta");
@@ -5567,7 +5566,7 @@ void CommonIntegrateBAOABStepKernel::initialize(const System& system, const BAOA
     prevStepSize = -1.0;
 }
 
-void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid) {
+void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator) {
     IntegrationUtilities& integration = cc.getIntegrationUtilities();
     int numAtoms = cc.getNumAtoms();
     int paddedNumAtoms = cc.getPaddedNumAtoms();
@@ -5577,11 +5576,8 @@ void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLa
         kernel1->addArg(paddedNumAtoms);
         kernel1->addArg(cc.getVelm());
         kernel1->addArg(cc.getLongForceBuffer());
-        kernel1->addArg(integration.getPosDelta());
-        kernel1->addArg(oldDelta);
         kernel1->addArg(integration.getStepSize());
         kernel2->addArg(numAtoms);
-        kernel2->addArg(cc.getPosq());
         kernel2->addArg(cc.getVelm());
         kernel2->addArg(integration.getPosDelta());
         kernel2->addArg(oldDelta);
@@ -5589,8 +5585,6 @@ void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLa
         kernel2->addArg(integration.getStepSize());
         kernel2->addArg(integration.getRandom());
         kernel2->addArg(); // Random index will be set just before it is executed.
-        if (cc.getUseMixedPrecision())
-            kernel2->addArg(cc.getPosqCorrection());
         kernel3->addArg(numAtoms);
         kernel3->addArg(cc.getPosq());
         kernel3->addArg(cc.getVelm());
@@ -5599,15 +5593,6 @@ void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLa
         kernel3->addArg(integration.getStepSize());
         if (cc.getUseMixedPrecision())
             kernel3->addArg(cc.getPosqCorrection());
-        kernel4->addArg(numAtoms);
-        kernel4->addArg(paddedNumAtoms);
-        kernel4->addArg(cc.getVelm());
-        kernel4->addArg(cc.getLongForceBuffer());
-        kernel4->addArg(integration.getStepSize());
-    }
-    if (!forcesAreValid) {
-        context.calcForcesAndEnergy(true, false);
-        forcesAreValid = true;
     }
     double temperature = integrator.getTemperature();
     double friction = integrator.getFriction();
@@ -5630,15 +5615,12 @@ void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLa
 
     // Perform the integration.
 
-    kernel2->setArg(8, integration.prepareRandomNumbers(cc.getPaddedNumAtoms()));
+    kernel2->setArg(7, integration.prepareRandomNumbers(cc.getPaddedNumAtoms()));
     kernel1->execute(numAtoms);
-    integration.applyConstraints(integrator.getConstraintTolerance());
+    integration.applyVelocityConstraints(integrator.getConstraintTolerance());
     kernel2->execute(numAtoms);
     integration.applyConstraints(integrator.getConstraintTolerance());
     kernel3->execute(numAtoms);
-    context.calcForcesAndEnergy(true, false);
-    kernel4->execute(numAtoms);
-    integration.applyVelocityConstraints(integrator.getConstraintTolerance());
     integration.computeVirtualSites();
 
     // Update the time and step count.
@@ -5646,8 +5628,6 @@ void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLa
     cc.setTime(cc.getTime()+stepSize);
     cc.setStepCount(cc.getStepCount()+1);
     cc.reorderAtoms();
-    if (cc.getAtomsWereReordered())
-        forcesAreValid = false;
     
     // Reduce UI lag.
     

platforms/common/src/kernels/baoab.cc

@@ -4,10 +4,9 @@ enum {VelScale, NoiseScale};
  * Perform the first part of BAOAB integration: velocity half step, then position half step.
  */
 
-KERNEL void integrateBAOABPart1(int numAtoms, int paddedNumAtoms, GLOBAL mixed4* RESTRICT velm, GLOBAL const mm_long* RESTRICT force, GLOBAL mixed4* RESTRICT posDelta,
-        GLOBAL mixed4* RESTRICT oldDelta, GLOBAL const mixed2* RESTRICT dt) {
-    mixed halfdt = 0.5*dt[0].y;
-    mixed fscale = halfdt/(mixed) 0x100000000;
+KERNEL void integrateBAOABPart1(int numAtoms, int paddedNumAtoms, GLOBAL mixed4* RESTRICT velm, GLOBAL const mm_long* RESTRICT force,
+        GLOBAL const mixed2* RESTRICT dt) {
+    mixed fscale = dt[0].y/(mixed) 0x100000000;
     for (int index = GLOBAL_ID; index < numAtoms; index += GLOBAL_SIZE) {
         mixed4 velocity = velm[index];
         if (velocity.w != 0.0) {
@@ -15,9 +14,6 @@ KERNEL void integrateBAOABPart1(int numAtoms, int paddedNumAtoms, GLOBAL mixed4*
             velocity.y += fscale*velocity.w*force[index+paddedNumAtoms];
             velocity.z += fscale*velocity.w*force[index+paddedNumAtoms*2];
             velm[index] = velocity;
-            mixed4 delta = make_mixed4(halfdt*velocity.x, halfdt*velocity.y, halfdt*velocity.z, 0);
-            posDelta[index] = delta;
-            oldDelta[index] = delta;
         }
     }
 }
@@ -27,47 +23,24 @@ KERNEL void integrateBAOABPart1(int numAtoms, int paddedNumAtoms, GLOBAL mixed4*
  * then position half step.
  */
 
-KERNEL void integrateBAOABPart2(int numAtoms, GLOBAL real4* RESTRICT posq, GLOBAL mixed4* RESTRICT velm, GLOBAL mixed4* RESTRICT posDelta,
+KERNEL void integrateBAOABPart2(int numAtoms, GLOBAL mixed4* RESTRICT velm, GLOBAL mixed4* RESTRICT posDelta,
         GLOBAL mixed4* RESTRICT oldDelta, GLOBAL const mixed* RESTRICT paramBuffer, GLOBAL const mixed2* RESTRICT dt, GLOBAL const float4* RESTRICT random, unsigned int randomIndex
-#ifdef USE_MIXED_PRECISION
-        , GLOBAL real4* RESTRICT posqCorrection
-#endif
         ) {
     mixed vscale = paramBuffer[VelScale];
     mixed noisescale = paramBuffer[NoiseScale];
     mixed halfdt = 0.5*dt[0].y;
-    mixed invHalfdt = 1/halfdt;
     int index = GLOBAL_ID;
     randomIndex += index;
     while (index < numAtoms) {
         mixed4 velocity = velm[index];
         if (velocity.w != 0.0) {
-            mixed4 delta = posDelta[index];
+            mixed4 delta = make_mixed4(halfdt*velocity.x, halfdt*velocity.y, halfdt*velocity.z, 0);
             mixed sqrtInvMass = SQRT(velocity.w);
-            velocity.x += (delta.x-oldDelta[index].x)*invHalfdt;
-            velocity.y += (delta.y-oldDelta[index].y)*invHalfdt;
-            velocity.z += (delta.z-oldDelta[index].z)*invHalfdt;
             velocity.x = vscale*velocity.x + noisescale*sqrtInvMass*random[randomIndex].x;
             velocity.y = vscale*velocity.y + noisescale*sqrtInvMass*random[randomIndex].y;
             velocity.z = vscale*velocity.z + noisescale*sqrtInvMass*random[randomIndex].z;
             velm[index] = velocity;
-#ifdef USE_MIXED_PRECISION
-            real4 pos1 = posq[index];
-            real4 pos2 = posqCorrection[index];
-            mixed4 pos = make_mixed4(pos1.x+(mixed)pos2.x, pos1.y+(mixed)pos2.y, pos1.z+(mixed)pos2.z, pos1.w);
-#else
-            real4 pos = posq[index];
-#endif
-            pos.x += delta.x;
-            pos.y += delta.y;
-            pos.z += delta.z;
-#ifdef USE_MIXED_PRECISION
-            posq[index] = make_real4((real) pos.x, (real) pos.y, (real) pos.z, (real) pos.w);
-            posqCorrection[index] = make_real4(pos.x-(real) pos.x, pos.y-(real) pos.y, pos.z-(real) pos.z, 0);
-#else
-            posq[index] = pos;
-#endif
-            delta = make_mixed4(halfdt*velocity.x, halfdt*velocity.y, halfdt*velocity.z, 0);
+            delta += make_mixed4(halfdt*velocity.x, halfdt*velocity.y, halfdt*velocity.z, 0);
             posDelta[index] = delta;
             oldDelta[index] = delta;
         }
@@ -87,15 +60,14 @@ KERNEL void integrateBAOABPart3(int numAtoms, GLOBAL real4* RESTRICT posq, GLOBA
         , GLOBAL real4* RESTRICT posqCorrection
 #endif
         ) {
-    mixed halfdt = 0.5*dt[0].y;
-    mixed invHalfdt = 1/halfdt;
+    mixed invDt = 1/dt[0].y;
     for (int index = GLOBAL_ID; index < numAtoms; index += GLOBAL_SIZE) {
         mixed4 velocity = velm[index];
         if (velocity.w != 0.0) {
             mixed4 delta = posDelta[index];
-            velocity.x += (delta.x-oldDelta[index].x)*invHalfdt;
-            velocity.y += (delta.y-oldDelta[index].y)*invHalfdt;
-            velocity.z += (delta.z-oldDelta[index].z)*invHalfdt;
+            velocity.x += (delta.x-oldDelta[index].x)*invDt;
+            velocity.y += (delta.y-oldDelta[index].y)*invDt;
+            velocity.z += (delta.z-oldDelta[index].z)*invDt;
             velm[index] = velocity;
 #ifdef USE_MIXED_PRECISION
             real4 pos1 = posq[index];
@@ -116,22 +88,3 @@ KERNEL void integrateBAOABPart3(int numAtoms, GLOBAL real4* RESTRICT posq, GLOBA
         }
     }
 }
-
-/**
- * Perform the fourth part of BAOAB integration: velocity half step.
- */
-
-KERNEL void integrateBAOABPart4(int numAtoms, int paddedNumAtoms, GLOBAL mixed4* RESTRICT velm,
-        GLOBAL const mm_long* RESTRICT force, GLOBAL const mixed2* RESTRICT dt) {
-    mixed halfdt = 0.5*dt[0].y;
-    mixed fscale = halfdt/(mixed) 0x100000000;
-    for (int index = GLOBAL_ID; index < numAtoms; index += GLOBAL_SIZE) {
-        mixed4 velocity = velm[index];
-        if (velocity.w != 0.0) {
-            velocity.x += fscale*velocity.w*force[index];
-            velocity.y += fscale*velocity.w*force[index+paddedNumAtoms];
-            velocity.z += fscale*velocity.w*force[index+paddedNumAtoms*2];
-            velm[index] = velocity;
-        }
-    }
-}

platforms/cpu/include/CpuBAOABDynamics.h

 
-/* Portions copyright (c) 2013-2019 Stanford University and Simbios.
+/* Portions copyright (c) 2013-2020 Stanford University and Simbios.
  * Authors: Peter Eastman
  * Contributors: 
  *
@@ -56,14 +56,11 @@ public:
      * First update step.
      * 
      * @param numberOfAtoms       number of atoms
-     * @param atomCoordinates     atom coordinates
      * @param velocities          velocities
      * @param forces              forces
      * @param inverseMasses       inverse atom masses
-     * @param xPrime              xPrime
      */
-    void updatePart1(int numberOfAtoms, std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& velocities,
-                     std::vector<OpenMM::Vec3>& forces, std::vector<double>& inverseMasses, std::vector<OpenMM::Vec3>& xPrime);
+    void updatePart1(int numberOfAtoms, std::vector<OpenMM::Vec3>& velocities,  std::vector<OpenMM::Vec3>& forces, std::vector<double>& inverseMasses);
       
     /**
      * Second update step.
@@ -94,7 +91,6 @@ private:
     void threadUpdate1(int threadIndex);
     void threadUpdate2(int threadIndex);
     void threadUpdate3(int threadIndex);
-    void threadUpdate4(int threadIndex);
     OpenMM::ThreadPool& threads;
     OpenMM::CpuRandom& random;
     std::vector<OpenMM_SFMT::SFMT> threadRandom;

platforms/cpu/include/CpuKernels.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) 2013-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -558,12 +558,8 @@ public:
      * 
      * @param context    the context in which to execute this kernel
      * @param integrator the BAOABLangevinIntegrator this kernel is being used for
-     * @param forcesAreValid if the context has been modified since the last time step, this will be
-     *                       false to show that cached forces are invalid and must be recalculated.
-     *                       On exit, this should specify whether the cached forces are valid at the
-     *                       end of the step.
      */
-    void execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid);
+    void execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator);
     /**
      * Compute the kinetic energy.
      * 

platforms/cpu/src/CpuBAOABDynamics.cpp

-
-/* Portions copyright (c) 2006-2019 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2020 Stanford University and Simbios.
  * Authors: Peter Eastman
  * Contributors: 
  *
@@ -36,16 +35,13 @@ CpuBAOABDynamics::CpuBAOABDynamics(int numberOfAtoms, double deltaT, double fric
 CpuBAOABDynamics::~CpuBAOABDynamics() {
 }
 
-void CpuBAOABDynamics::updatePart1(int numberOfAtoms, vector<Vec3>& atomCoordinates, vector<Vec3>& velocities,
-                                   vector<Vec3>& forces, vector<double>& inverseMasses, vector<Vec3>& xPrime) {
+void CpuBAOABDynamics::updatePart1(int numberOfAtoms, vector<Vec3>& velocities, vector<Vec3>& forces, vector<double>& inverseMasses) {
     // Record the parameters for the threads.
     
     this->numberOfAtoms = numberOfAtoms;
-    this->atomCoordinates = &atomCoordinates[0];
     this->velocities = &velocities[0];
     this->forces = &forces[0];
     this->inverseMasses = &inverseMasses[0];
-    this->xPrime = &xPrime[0];
     
     // Signal the threads to start running and wait for them to finish.
     
@@ -83,23 +79,15 @@ void CpuBAOABDynamics::updatePart3(ContextImpl& context, int numberOfAtoms, vect
     
     threads.execute([&] (ThreadPool& threads, int threadIndex) { threadUpdate3(threadIndex); });
     threads.waitForThreads();
-    context.calcForcesAndEnergy(true, false);
-    threads.execute([&] (ThreadPool& threads, int threadIndex) { threadUpdate4(threadIndex); });
-    threads.waitForThreads();
 }
 
 void CpuBAOABDynamics::threadUpdate1(int threadIndex) {
-    const double halfdt = 0.5*getDeltaT();
     int start = threadIndex*numberOfAtoms/threads.getNumThreads();
     int end = (threadIndex+1)*numberOfAtoms/threads.getNumThreads();
 
-    for (int i = start; i < end; i++) {
-        if (inverseMasses[i] != 0.0) {
-            velocities[i] += (halfdt*inverseMasses[i])*forces[i];
-            xPrime[i] = atomCoordinates[i] + velocities[i]*halfdt;
-            oldx[i] = xPrime[i];
-        }
-    }
+    for (int i = start; i < end; i++)
+        if (inverseMasses[i] != 0.0)
+            velocities[i] += (getDeltaT()*inverseMasses[i])*forces[i];
 }
 
 void CpuBAOABDynamics::threadUpdate2(int threadIndex) {
@@ -113,34 +101,22 @@ void CpuBAOABDynamics::threadUpdate2(int threadIndex) {
 
     for (int i = start; i < end; i++) {
         if (inverseMasses[i] != 0.0) {
-            velocities[i] += (xPrime[i]-oldx[i])/halfdt;
+            xPrime[i] = atomCoordinates[i] + velocities[i]*halfdt;
             Vec3 noise(random.getGaussianRandom(threadIndex), random.getGaussianRandom(threadIndex), random.getGaussianRandom(threadIndex));
             velocities[i] = vscale*velocities[i] + noisescale*sqrt(kT*inverseMasses[i])*noise;
-            atomCoordinates[i] = xPrime[i];
-            xPrime[i] = atomCoordinates[i] + velocities[i]*halfdt;
+            xPrime[i] = xPrime[i] + velocities[i]*halfdt;
             oldx[i] = xPrime[i];
         }
     }
 }
 
 void CpuBAOABDynamics::threadUpdate3(int threadIndex) {
-    const double halfdt = 0.5*getDeltaT();
     int start = threadIndex*numberOfAtoms/threads.getNumThreads();
     int end = (threadIndex+1)*numberOfAtoms/threads.getNumThreads();
 
     for (int i = start; i < end; ++i)
         if (inverseMasses[i] != 0.0) {
-            velocities[i] += (xPrime[i]-oldx[i])/halfdt;
+            velocities[i] += (xPrime[i]-oldx[i])/getDeltaT();
             atomCoordinates[i] = xPrime[i];
         }
 }
-
-void CpuBAOABDynamics::threadUpdate4(int threadIndex) {
-    const double halfdt = 0.5*getDeltaT();
-    int start = threadIndex*numberOfAtoms/threads.getNumThreads();
-    int end = (threadIndex+1)*numberOfAtoms/threads.getNumThreads();
-
-    for (int i = start; i < end; ++i)
-        if (inverseMasses[i] != 0.0)
-            velocities[i] += (halfdt*inverseMasses[i])*forces[i];
-}

platforms/cpu/src/CpuKernels.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) 2013-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -1366,7 +1366,7 @@ void CpuIntegrateBAOABStepKernel::initialize(const System& system, const BAOABLa
     data.random.initialize(integrator.getRandomNumberSeed(), data.threads.getNumThreads());
 }
 
-void CpuIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid) {
+void CpuIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator) {
     double temperature = integrator.getTemperature();
     double friction = integrator.getFriction();
     double stepSize = integrator.getStepSize();
@@ -1383,7 +1383,7 @@ void CpuIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLange
         prevFriction = friction;
         prevStepSize = stepSize;
     }
-    dynamics->update(context, posData, velData, masses, forcesAreValid, integrator.getConstraintTolerance());
+    dynamics->update(context, posData, velData, masses, integrator.getConstraintTolerance());
     ReferencePlatform::PlatformData* refData = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
     refData->time += stepSize;
     refData->stepCount++;

platforms/reference/include/ReferenceBAOABDynamics.h

 
-/* Portions copyright (c) 2006-2019 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2020 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -79,29 +79,25 @@ class OPENMM_EXPORT ReferenceBAOABDynamics : public ReferenceDynamics {
          @param atomCoordinates     atom coordinates
          @param velocities          velocities
          @param masses              atom masses
-         @param forcesAreValid      whether the current forces are valid or need to be recomputed
          @param tolerance           the constraint tolerance
       
          --------------------------------------------------------------------------------------- */
      
       void update(OpenMM::ContextImpl& context, std::vector<OpenMM::Vec3>& atomCoordinates,
-                  std::vector<OpenMM::Vec3>& velocities, std::vector<double>& masses, bool& forcesAreValid, double tolerance);
+                  std::vector<OpenMM::Vec3>& velocities, std::vector<double>& masses, double tolerance);
      
       /**---------------------------------------------------------------------------------------
       
          First update; based on code in update.c do_update_sd() Gromacs 3.1.4
       
          @param numberOfAtoms       number of atoms
-         @param atomCoordinates     atom coordinates
          @param velocities          velocities
          @param forces              forces
          @param inverseMasses       inverse atom masses
-         @param xPrime              xPrime
       
          --------------------------------------------------------------------------------------- */
       
-      virtual void updatePart1(int numberOfAtoms, std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& velocities,
-                       std::vector<OpenMM::Vec3>& forces, std::vector<double>& inverseMasses, std::vector<OpenMM::Vec3>& xPrime);
+      virtual void updatePart1(int numberOfAtoms, std::vector<OpenMM::Vec3>& velocities,  std::vector<OpenMM::Vec3>& forces, std::vector<double>& inverseMasses);
       
       /**---------------------------------------------------------------------------------------
       
@@ -126,14 +122,13 @@ class OPENMM_EXPORT ReferenceBAOABDynamics : public ReferenceDynamics {
          @param numberOfAtoms       number of atoms
          @param atomCoordinates     atom coordinates
          @param velocities          velocities
-         @param forces              forces
          @param inverseMasses       inverse atom masses
          @param xPrime              xPrime
       
          --------------------------------------------------------------------------------------- */
       
       virtual void updatePart3(OpenMM::ContextImpl& context, int numberOfAtoms, std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& velocities,
-                       std::vector<OpenMM::Vec3>& forces, std::vector<double>& inverseMasses, std::vector<OpenMM::Vec3>& xPrime);
+                               std::vector<double>& inverseMasses, std::vector<OpenMM::Vec3>& xPrime);
 };
 
 } // namespace OpenMM

platforms/reference/include/ReferenceKernels.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-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -1233,12 +1233,8 @@ public:
      * 
      * @param context    the context in which to execute this kernel
      * @param integrator the BAOABLangevinIntegrator this kernel is being used for
-     * @param forcesAreValid if the context has been modified since the last time step, this will be
-     *                       false to show that cached forces are invalid and must be recalculated.
-     *                       On exit, this should specify whether the cached forces are valid at the
-     *                       end of the step.
      */
-    void execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid);
+    void execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator);
     /**
      * Compute the kinetic energy.
      * 

platforms/reference/src/ReferenceKernels.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-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -2238,7 +2238,7 @@ void ReferenceIntegrateBAOABStepKernel::initialize(const System& system, const B
     SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
 }
 
-void ReferenceIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid) {
+void ReferenceIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator) {
     double temperature = integrator.getTemperature();
     double friction = integrator.getFriction();
     double stepSize = integrator.getStepSize();
@@ -2259,7 +2259,7 @@ void ReferenceIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOA
         prevFriction = friction;
         prevStepSize = stepSize;
     }
-    dynamics->update(context, posData, velData, masses, forcesAreValid, integrator.getConstraintTolerance());
+    dynamics->update(context, posData, velData, masses, integrator.getConstraintTolerance());
     data.time += stepSize;
     data.stepCount++;
 }

platforms/reference/src/SimTKReference/ReferenceBAOABDynamics.cpp

 
-/* Portions copyright (c) 2006-2019 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2020 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -51,17 +51,10 @@ double ReferenceBAOABDynamics::getFriction() const {
    return friction;
 }
 
-void ReferenceBAOABDynamics::updatePart1(int numberOfAtoms, vector<Vec3>& atomCoordinates,
-                                              vector<Vec3>& velocities, vector<Vec3>& forces, vector<double>& inverseMasses,
-                                              vector<Vec3>& xPrime) {
-    const double halfdt = 0.5*getDeltaT();
-    for (int i = 0; i < numberOfAtoms; i++) {
-        if (inverseMasses[i] != 0.0) {
-            velocities[i] += (halfdt*inverseMasses[i])*forces[i];
-            xPrime[i] = atomCoordinates[i] + velocities[i]*halfdt;
-            oldx[i] = xPrime[i];
-        }
-    }
+void ReferenceBAOABDynamics::updatePart1(int numberOfAtoms, vector<Vec3>& velocities, vector<Vec3>& forces, vector<double>& inverseMasses) {
+    for (int i = 0; i < numberOfAtoms; i++)
+        if (inverseMasses[i] != 0.0)
+            velocities[i] += (getDeltaT()*inverseMasses[i])*forces[i];
 }
 
 void ReferenceBAOABDynamics::updatePart2(int numberOfAtoms, vector<Vec3>& atomCoordinates,
@@ -75,37 +68,29 @@ void ReferenceBAOABDynamics::updatePart2(int numberOfAtoms, vector<Vec3>& atomCo
 
     for (int i = 0; i < numberOfAtoms; i++) {
         if (inverseMasses[i] != 0.0) {
-            velocities[i] += (xPrime[i]-oldx[i])/halfdt;
+            xPrime[i] = atomCoordinates[i] + velocities[i]*halfdt;
             velocities[i] = vscale*velocities[i] + noisescale*sqrt(kT*inverseMasses[i])*Vec3(
                     SimTKOpenMMUtilities::getNormallyDistributedRandomNumber(),
                     SimTKOpenMMUtilities::getNormallyDistributedRandomNumber(),
                     SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
-            atomCoordinates[i] = xPrime[i];
-            xPrime[i] = atomCoordinates[i] + velocities[i]*halfdt;
+            xPrime[i] = xPrime[i] + velocities[i]*halfdt;
             oldx[i] = xPrime[i];
         }
     }
 }
 
 void ReferenceBAOABDynamics::updatePart3(OpenMM::ContextImpl& context, int numberOfAtoms, vector<Vec3>& atomCoordinates,
-                                              vector<Vec3>& velocities, vector<Vec3>& forces, vector<double>& inverseMasses,
-                                              vector<Vec3>& xPrime) {
-    const double halfdt = 0.5*getDeltaT();
+                                         vector<Vec3>& velocities, vector<double>& inverseMasses, vector<Vec3>& xPrime) {
     for (int i = 0; i < numberOfAtoms; i++) {
         if (inverseMasses[i] != 0.0) {
-            velocities[i] += (xPrime[i]-oldx[i])/halfdt;
+            velocities[i] += (xPrime[i]-oldx[i])/getDeltaT();
             atomCoordinates[i] = xPrime[i];
         }
     }
-    context.calcForcesAndEnergy(true, false);
-    for (int i = 0; i < numberOfAtoms; i++) {
-        if (inverseMasses[i] != 0.0)
-            velocities[i] += (halfdt*inverseMasses[i])*forces[i];
-    }
 }
 
 void ReferenceBAOABDynamics::update(ContextImpl& context, vector<Vec3>& atomCoordinates,
-                                          vector<Vec3>& velocities, vector<double>& masses, bool& forcesAreValid, double tolerance) {
+                                    vector<Vec3>& velocities, vector<double>& masses, double tolerance) {
     int numberOfAtoms = context.getSystem().getNumParticles();
     ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
     if (getTimeStep() == 0) {
@@ -118,18 +103,14 @@ void ReferenceBAOABDynamics::update(ContextImpl& context, vector<Vec3>& atomCoor
                 inverseMasses[ii] = 1.0/masses[ii];
         }
     }
-    if (!forcesAreValid) {
-        context.calcForcesAndEnergy(true, false);
-        forcesAreValid = true;
-    }
     ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
     vector<Vec3>& forces = *data->forces;
 
     // 1st update
 
-    updatePart1(numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses, xPrime);
+    updatePart1(numberOfAtoms, velocities, forces, inverseMasses);
     if (referenceConstraintAlgorithm)
-        referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses, tolerance);
+        referenceConstraintAlgorithm->applyToVelocities(atomCoordinates, velocities, inverseMasses, tolerance);
 
     // 2nd update
 
@@ -139,9 +120,7 @@ void ReferenceBAOABDynamics::update(ContextImpl& context, vector<Vec3>& atomCoor
 
     // 3rd update
 
-    updatePart3(context, numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses, xPrime);
-    if (referenceConstraintAlgorithm)
-        referenceConstraintAlgorithm->applyToVelocities(atomCoordinates, velocities, inverseMasses, tolerance);
+    updatePart3(context, numberOfAtoms, atomCoordinates, velocities, inverseMasses, xPrime);
 
     ReferenceVirtualSites::computePositions(context.getSystem(), atomCoordinates);
     incrementTimeStep();