Use LF-Middle for LangevinIntegrator and VariableLangevinIntegrator (#4440)

* Made LangevinIntegrator identical to LangevinMiddleIntegrator

* Removed unused code

* VariableLangevinIntegrator uses LFMiddle

platforms/reference/include/ReferenceVariableStochasticDynamics.h

 
-/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2024 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -33,7 +33,7 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
 
    private:
 
-      std::vector<OpenMM::Vec3> xPrime;
+      std::vector<OpenMM::Vec3> xPrime, oldx;
       std::vector<double> inverseMasses;
       double friction, _accuracy;
 
@@ -104,39 +104,31 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
                   std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, double maxStepSize, 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 masses              atom masses
          @param inverseMasses       inverse atom masses
-         @param xPrime              xPrime
          @param maxStepSize         maximum time step
-
+      
          --------------------------------------------------------------------------------------- */
-
-      void updatePart1(int numberOfAtoms, std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& velocities,
-                       std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, std::vector<double>& inverseMasses,
-                       std::vector<OpenMM::Vec3>& xPrime, double maxStepSize);
+      
+      void updatePart1(int numberOfAtoms, std::vector<OpenMM::Vec3>& velocities,  std::vector<OpenMM::Vec3>& forces, std::vector<double>& inverseMasses, double maxStepSize);
 
       /**---------------------------------------------------------------------------------------
-
+      
          Second update
-
+      
          @param numberOfAtoms       number of atoms
          @param atomCoordinates     atom coordinates
          @param velocities          velocities
-         @param forces              forces
-         @param masses              atom masses
-
+         @param inverseMasses       inverse atom masses
+         @param xPrime              xPrime
+      
          --------------------------------------------------------------------------------------- */
-
+      
       void updatePart2(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);
       
 };
 

platforms/reference/src/ReferenceKernelFactory.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-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -94,8 +94,6 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceIntegrateVerletStepKernel(name, platform, data);
     if (name == IntegrateNoseHooverStepKernel::Name())
         return new ReferenceIntegrateNoseHooverStepKernel(name, platform, data);
-    if (name == IntegrateLangevinStepKernel::Name())
-        return new ReferenceIntegrateLangevinStepKernel(name, platform, data);
     if (name == IntegrateLangevinMiddleStepKernel::Name())
         return new ReferenceIntegrateLangevinMiddleStepKernel(name, platform, data);
     if (name == IntegrateBrownianStepKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -61,7 +61,6 @@
 #include "ReferenceProperDihedralBond.h"
 #include "ReferenceRbDihedralBond.h"
 #include "ReferenceRMSDForce.h"
-#include "ReferenceStochasticDynamics.h"
 #include "ReferenceTabulatedFunction.h"
 #include "ReferenceVariableStochasticDynamics.h"
 #include "ReferenceVariableVerletDynamics.h"
@@ -2570,51 +2569,6 @@ void ReferenceIntegrateNoseHooverStepKernel::setChainStates(ContextImpl& context
     chainVelocities = velocities;
 }
 
-ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
-    if (dynamics)
-        delete dynamics;
-}
-
-void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
-    int numParticles = system.getNumParticles();
-    masses.resize(numParticles);
-    for (int i = 0; i < numParticles; ++i)
-        masses[i] = system.getParticleMass(i);
-    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
-}
-
-void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
-    double temperature = integrator.getTemperature();
-    double friction = integrator.getFriction();
-    double stepSize = integrator.getStepSize();
-    vector<Vec3>& posData = extractPositions(context);
-    vector<Vec3>& velData = extractVelocities(context);
-    vector<Vec3>& forceData = extractForces(context);
-    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
-        // Recreate the computation objects with the new parameters.
-        
-        if (dynamics)
-            delete dynamics;
-        dynamics = new ReferenceStochasticDynamics(
-                context.getSystem().getNumParticles(), 
-                stepSize, 
-                friction, 
-                temperature);
-        dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
-        dynamics->setVirtualSites(extractVirtualSites(context));
-        prevTemp = temperature;
-        prevFriction = friction;
-        prevStepSize = stepSize;
-    }
-    dynamics->update(context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance());
-    data.time += stepSize;
-    data.stepCount++;
-}
-
-double ReferenceIntegrateLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const LangevinIntegrator& integrator) {
-    return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
-}
-
 ReferenceIntegrateLangevinMiddleStepKernel::~ReferenceIntegrateLangevinMiddleStepKernel() {
     if (dynamics)
         delete dynamics;

platforms/reference/src/ReferencePlatform.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-2023 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -69,7 +69,6 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateNoseHooverStepKernel::Name(), factory);
-    registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinMiddleStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
     registerKernelFactory(IntegrateVariableLangevinStepKernel::Name(), factory);

platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.cpp

-
-/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
- * Contributors: Pande Group
- *
- * Permission is hereby granted, free of charge, to any person obtaining
- * a copy of this software and associated documentation files (the
- * "Software"), to deal in the Software without restriction, including
- * without limitation the rights to use, copy, modify, merge, publish,
- * distribute, sublicense, and/or sell copies of the Software, and to
- * permit persons to whom the Software is furnished to do so, subject
- * to the following conditions:
- *
- * The above copyright notice and this permission notice shall be included
- * in all copies or substantial portions of the Software.
- *
- * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
- * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
- * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
- * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
- * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
- * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
- */
-
-#include <cstring>
-#include <sstream>
-
-#include "SimTKOpenMMUtilities.h"
-#include "ReferenceStochasticDynamics.h"
-#include "ReferenceVirtualSites.h"
-#include "openmm/OpenMMException.h"
-
-#include <cstdio>
-
-using std::vector;
-using namespace OpenMM;
-
-/**---------------------------------------------------------------------------------------
-
-   ReferenceStochasticDynamics constructor
-
-   @param numberOfAtoms  number of atoms
-   @param deltaT         delta t for dynamics
-   @param friction       friction coefficient
-   @param temperature    temperature
-
-   --------------------------------------------------------------------------------------- */
-
-ReferenceStochasticDynamics::ReferenceStochasticDynamics(int numberOfAtoms,
-                                                         double deltaT, double friction,
-                                                         double temperature) : 
-           ReferenceDynamics(numberOfAtoms, deltaT, temperature), friction(friction) {
-   xPrime.resize(numberOfAtoms);
-   inverseMasses.resize(numberOfAtoms);
-}
-
-/**---------------------------------------------------------------------------------------
-
-   ReferenceStochasticDynamics destructor
-
-   --------------------------------------------------------------------------------------- */
-
-ReferenceStochasticDynamics::~ReferenceStochasticDynamics() {
-}
-
-/**---------------------------------------------------------------------------------------
-
-   Get friction coefficient
-
-   --------------------------------------------------------------------------------------- */
-
-double ReferenceStochasticDynamics::getFriction() const {
-   return friction;
-}
-
-/**---------------------------------------------------------------------------------------
-
-   First SD 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
-
-   --------------------------------------------------------------------------------------- */
-
-void ReferenceStochasticDynamics::updatePart1(int numberOfAtoms, vector<Vec3>& atomCoordinates,
-                                              vector<Vec3>& velocities,
-                                              vector<Vec3>& forces, vector<double>& inverseMasses,
-                                              vector<Vec3>& xPrime) {
-   // perform first update
-
-   double dt = getDeltaT();
-   double friction = getFriction();
-   const double vscale = exp(-dt*friction);
-   const double fscale = (friction == 0 ? dt : (1-vscale)/friction);
-   const double kT = BOLTZ*getTemperature();
-   const double noisescale = sqrt(kT*(1-vscale*vscale));
-
-   for (int ii = 0; ii < numberOfAtoms; ii++) {
-       if (inverseMasses[ii] != 0.0) {
-           double sqrtInvMass = sqrt(inverseMasses[ii]);
-           for (int jj = 0; jj < 3; jj++) {
-               velocities[ii][jj]  = vscale*velocities[ii][jj] + fscale*inverseMasses[ii]*forces[ii][jj] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
-           }
-       }
-   }
-}
-
-/**---------------------------------------------------------------------------------------
-
-   Second update; based on code in update.c do_update_sd() w/ bFirstHalf = false in Gromacs 3.1.4
-
-   @param numberOfAtoms       number of atoms
-   @param atomCoordinates     atom coordinates
-   @param velocities          velocities
-   @param forces              forces
-   @param masses              atom masses
-
-   --------------------------------------------------------------------------------------- */
-
-void ReferenceStochasticDynamics::updatePart2(int numberOfAtoms, vector<Vec3>& atomCoordinates,
-                                              vector<Vec3>& velocities,
-                                              vector<Vec3>& forces, vector<double>& inverseMasses,
-                                              vector<Vec3>& xPrime) {
-   // perform second update
-
-   for (int ii = 0; ii < numberOfAtoms; ii++) {
-       if (inverseMasses[ii] != 0.0)
-            xPrime[ii] = atomCoordinates[ii]+velocities[ii]*getDeltaT();
-   }
-}
-
-void ReferenceStochasticDynamics::updatePart3(int numberOfAtoms, vector<Vec3>& atomCoordinates,
-                                              vector<Vec3>& velocities, vector<double>& inverseMasses,
-                                              vector<Vec3>& xPrime) {
-   double invStepSize = 1.0/getDeltaT();
-   for (int i = 0; i < numberOfAtoms; ++i)
-       if (inverseMasses[i] != 0) {
-            velocities[i] = (xPrime[i]-atomCoordinates[i])*invStepSize;
-            atomCoordinates[i] = xPrime[i];
-       }
-}
-
-/**---------------------------------------------------------------------------------------
-
-   Update -- driver routine for performing stochastic dynamics update of coordinates
-   and velocities
-
-   @param system              the System to be integrated
-   @param atomCoordinates     atom coordinates
-   @param velocities          velocities
-   @param forces              forces
-   @param masses              atom masses
-
-   --------------------------------------------------------------------------------------- */
-
-void ReferenceStochasticDynamics::update(const OpenMM::System& system, vector<Vec3>& atomCoordinates,
-                                          vector<Vec3>& velocities, vector<Vec3>& forces, vector<double>& masses, double tolerance) {
-   // first-time-through initialization
-
-   int numberOfAtoms = system.getNumParticles();
-   if (getTimeStep() == 0) {
-      // invert masses
-
-      for (int ii = 0; ii < numberOfAtoms; ii++) {
-         if (masses[ii] == 0.0)
-             inverseMasses[ii] = 0.0;
-         else
-             inverseMasses[ii] = 1.0/masses[ii];
-      }
-   }
-
-   // 1st update
-
-   updatePart1(numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses, xPrime);
-
-   // 2nd update
-
-   updatePart2(numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses, xPrime);
-
-   ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
-   if (referenceConstraintAlgorithm)
-      referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses, tolerance);
-
-   // copy xPrime -> atomCoordinates
-
-   updatePart3(numberOfAtoms, atomCoordinates, velocities, inverseMasses, xPrime);
-
-   getVirtualSites().computePositions(system, atomCoordinates);
-   incrementTimeStep();
-}

platforms/reference/src/SimTKReference/ReferenceVariableStochasticDynamics.cpp

 
-/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2024 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -53,6 +53,7 @@ ReferenceVariableStochasticDynamics::ReferenceVariableStochasticDynamics(int num
                                                           double accuracy) :
            ReferenceDynamics(numberOfAtoms, 0.0f, temperature), friction(friction), _accuracy(accuracy) {
    xPrime.resize(numberOfAtoms);
+   oldx.resize(numberOfAtoms);
    inverseMasses.resize(numberOfAtoms);
 }
 
@@ -97,38 +98,7 @@ double ReferenceVariableStochasticDynamics::getFriction() const {
    return friction;
 }
 
-/**---------------------------------------------------------------------------------------
-
-   First SD 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 masses              atom masses
-   @param inverseMasses       inverse atom masses
-   @param xPrime              xPrime
-   @param maxStepSize         maximum time step
-
-   --------------------------------------------------------------------------------------- */
-
-void ReferenceVariableStochasticDynamics::updatePart1(int numberOfAtoms, vector<Vec3>& atomCoordinates,
-                                              vector<Vec3>& velocities,
-                                              vector<Vec3>& forces, vector<double>& masses, vector<double>& inverseMasses,
-                                              vector<Vec3>& xPrime, double maxStepSize) {
-   // first-time-through initialization
-
-   if (getTimeStep() == 0) {
-      // invert masses
-
-      for (int ii = 0; ii < numberOfAtoms; ii++) {
-         if (masses[ii] == 0)
-             inverseMasses[ii] = 0;
-         else
-             inverseMasses[ii] = 1/masses[ii];
-      }
-   }
-
+void ReferenceVariableStochasticDynamics::updatePart1(int numberOfAtoms, vector<Vec3>& velocities, vector<Vec3>& forces, vector<double>& inverseMasses, double maxStepSize) {
    // Select the step size to use
     double error = 0;
     for (int i = 0; i < numberOfAtoms; ++i) {
@@ -138,59 +108,42 @@ void ReferenceVariableStochasticDynamics::updatePart1(int numberOfAtoms, vector<
         }
     }
     error = sqrt(error/(numberOfAtoms*3));
-    double newStepSize = sqrt(getAccuracy()/error);
+    double dt = sqrt(getAccuracy()/error);
     if (getDeltaT() > 0.0f)
-        newStepSize = std::min(newStepSize, getDeltaT()*2.0f); // For safety, limit how quickly dt can increase.
-    if (newStepSize > getDeltaT() && newStepSize < 1.2f*getDeltaT())
-        newStepSize = getDeltaT(); // Keeping dt constant between steps improves the behavior of the integrator.
-    if (newStepSize > maxStepSize)
-        newStepSize = maxStepSize;
-    setDeltaT(newStepSize);
+        dt = std::min(dt, getDeltaT()*2.0f); // For safety, limit how quickly dt can increase.
+    if (dt > getDeltaT() && dt < 1.2f*getDeltaT())
+        dt = getDeltaT(); // Keeping dt constant between steps improves the behavior of the integrator.
+    if (dt > maxStepSize)
+        dt = maxStepSize;
+    setDeltaT(dt);
  
     // perform first update
 
-   double dt = getDeltaT();
-   double friction = getFriction();
-   const double vscale = exp(-dt*friction);
-   const double fscale = (friction == 0 ? dt : (1-vscale)/friction);
-   const double kT = BOLTZ*getTemperature();
-   const double noisescale = sqrt(kT*(1-vscale*vscale));
-
-   for (int ii = 0; ii < numberOfAtoms; ii++) {
-       if (masses[ii] != 0) {
-           double sqrtInvMass = sqrt(inverseMasses[ii]);
-           for (int jj = 0; jj < 3; jj++) {
-               velocities[ii][jj]  = vscale*velocities[ii][jj] + fscale*inverseMasses[ii]*forces[ii][jj] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
-           }
-       }
-   }
-
+    for (int i = 0; i < numberOfAtoms; i++)
+        if (inverseMasses[i] != 0.0)
+            velocities[i] += (dt*inverseMasses[i])*forces[i];
 }
 
-/**---------------------------------------------------------------------------------------
-
-   Second update; based on code in update.c do_update_sd() w/ bFirstHalf = false in Gromacs 3.1.4
-
-   @param numberOfAtoms       number of atoms
-   @param atomCoordinates     atom coordinates
-   @param velocities          velocities
-   @param forces              forces
-   @param masses              atom masses
-
-   --------------------------------------------------------------------------------------- */
-
 void ReferenceVariableStochasticDynamics::updatePart2(int numberOfAtoms, vector<Vec3>& atomCoordinates,
-                                              vector<Vec3>& velocities,
-                                              vector<Vec3>& forces, vector<double>& inverseMasses,
-                                              vector<Vec3>& xPrime) {
-   // perform second update
-
-   for (int ii = 0; ii < numberOfAtoms; ii++) {
-       if (inverseMasses[ii] != 0.0)
-           for (int jj = 0; jj < 3; jj++)
-               xPrime[ii][jj] = atomCoordinates[ii][jj]+getDeltaT()*velocities[ii][jj];
-   }
-
+                                         vector<Vec3>& velocities, vector<double>& inverseMasses,
+                                         vector<Vec3>& xPrime) {
+    const double halfdt = 0.5*getDeltaT();
+    const double kT = BOLTZ*getTemperature();
+    const double friction = getFriction();
+    const double vscale = exp(-getDeltaT()*friction);
+    const double noisescale = sqrt(1-vscale*vscale);
+
+    for (int i = 0; i < numberOfAtoms; i++) {
+        if (inverseMasses[i] != 0.0) {
+            xPrime[i] = atomCoordinates[i] + velocities[i]*halfdt;
+            velocities[i] = vscale*velocities[i] + noisescale*sqrt(kT*inverseMasses[i])*Vec3(
+                    SimTKOpenMMUtilities::getNormallyDistributedRandomNumber(),
+                    SimTKOpenMMUtilities::getNormallyDistributedRandomNumber(),
+                    SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
+            xPrime[i] = xPrime[i] + velocities[i]*halfdt;
+            oldx[i] = xPrime[i];
+        }
+    }
 }
 
 /**---------------------------------------------------------------------------------------
@@ -209,30 +162,41 @@ void ReferenceVariableStochasticDynamics::updatePart2(int numberOfAtoms, vector<
 void ReferenceVariableStochasticDynamics::update(const OpenMM::System& system, vector<Vec3>& atomCoordinates,
                                           vector<Vec3>& velocities,
                                           vector<Vec3>& forces, vector<double>& masses, double maxStepSize, double tolerance) {
+    int numberOfAtoms = system.getNumParticles();
+    ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
+    if (getTimeStep() == 0) {
+        // Invert masses
+
+        for (int ii = 0; ii < numberOfAtoms; ii++) {
+            if (masses[ii] == 0.0)
+                inverseMasses[ii] = 0.0;
+            else
+                inverseMasses[ii] = 1/masses[ii];
+        }
+    }
 
-   // 1st update
-
-   int numberOfAtoms = system.getNumParticles();
-   updatePart1(numberOfAtoms, atomCoordinates, velocities, forces, masses, inverseMasses, xPrime, maxStepSize);
+    // 1st update
 
-   // 2nd update
+    updatePart1(numberOfAtoms, velocities, forces, inverseMasses, maxStepSize);
+    if (referenceConstraintAlgorithm)
+        referenceConstraintAlgorithm->applyToVelocities(atomCoordinates, velocities, inverseMasses, tolerance);
 
-   updatePart2(numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses, xPrime);
+    // 2nd update
 
-   ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
-   if (referenceConstraintAlgorithm)
-      referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses, tolerance);
+    updatePart2(numberOfAtoms, atomCoordinates, velocities, inverseMasses, xPrime);
+    if (referenceConstraintAlgorithm)
+        referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses, tolerance);
 
-   // copy xPrime -> atomCoordinates
+    // copy xPrime -> atomCoordinates
 
-   double invStepSize = 1.0/getDeltaT();
-   for (int ii = 0; ii < numberOfAtoms; ii++) {
-       if (masses[ii] != 0.0) {
-           velocities[ii] = (xPrime[ii]-atomCoordinates[ii])*invStepSize;
-           atomCoordinates[ii] = xPrime[ii];
-       }
-   }
+    double invStepSize = 1.0/getDeltaT();
+    for (int i = 0; i < numberOfAtoms; i++) {
+        if (masses[i] != 0.0) {
+            velocities[i] += (xPrime[i]-oldx[i])*invStepSize;
+            atomCoordinates[i] = xPrime[i];
+        }
+    }
 
-   getVirtualSites().computePositions(system, atomCoordinates);
-   incrementTimeStep();
+    getVirtualSites().computePositions(system, atomCoordinates);
+    incrementTimeStep();
 }