Nosé Hoover Middle scheme (#2600)

* Convert Nose-Hoover into LF middle scheme by copying NH kernels

* Rebrand VelocityVerletIntegrator as NoseHooverIntegrator

* Consolidate NH tests

* NoseHooverChainKernel begone

* Make Windows builds happy

* Add missing header for Windows build

* Fix mistake in CommonKernels header

* Add 6th Yoshida-Suzuki and make it the default

olla/include/openmm/kernels.h

@@ -1063,12 +1063,12 @@ public:
 /**
  * This kernel is invoked by NoseHooverIntegrator to take one time step.
  */
-class IntegrateVelocityVerletStepKernel : public KernelImpl {
+class IntegrateNoseHooverStepKernel : public KernelImpl {
 public:
     static std::string Name() {
-        return "IntegrateVelocityVerletStep";
+        return "IntegrateNoseHooverStep";
     }
-    IntegrateVelocityVerletStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    IntegrateNoseHooverStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
     }
     /**
      * Initialize the kernel.
@@ -1093,6 +1093,41 @@ public:
      * @param integrator the NoseHooverIntegrator this kernel is being used for
      */
     virtual double computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator) = 0;
+    /**
+     * Execute the kernel that propagates the Nose Hoover chain and determines the velocity scale factor.
+     *
+     * @param context  the context in which to execute this kernel
+     * @param noseHooverChain the object describing the chain to be propagated.
+     * @param kineticEnergy the {center of mass, relative} kineticEnergies of the particles being thermostated by this chain.
+     * @param timeStep the time step used by the integrator.
+     * @return the velocity scale factor to apply to the particles associated with this heat bath.
+     */
+    virtual std::pair<double, double> propagateChain(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> kineticEnergy, double timeStep) = 0;
+    /**
+     * Execute the kernal that computes the total (kinetic + potential) heat bath energy.
+     *
+     * @param context the context in which to execute this kernel
+     * @param noseHooverChain the chain whose energy is to be determined.
+     * @return the total heat bath energy.
+     */
+    virtual double computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain) = 0;
+    /**
+     * Execute the kernel that computes the kinetic energy for a subset of atoms,
+     * or the relative kinetic energy of Drude particles with respect to their parent atoms
+     *
+     * @param context the context in which to execute this kernel
+     * @param noseHooverChain the chain whose energy is to be determined.
+     * @param downloadValue whether the computed value should be downloaded and returned.
+     */
+    virtual std::pair<double, double> computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue) = 0;
+    /**
+     * Execute the kernel that scales the velocities of particles associated with a nose hoover chain
+     *
+     * @param context the context in which to execute this kernel
+     * @param noseHooverChain the chain whose energy is to be determined.
+     * @param scaleFactor the multiplicative factor by which {absolute, relative} velocities are scaled.
+     */
+    virtual void scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactor) = 0;
 };
 
 /**
@@ -1360,58 +1395,6 @@ public:
     virtual void execute(ContextImpl& context) = 0;
 };
 
-/**
- * This kernel is invoked by NoseHooverChainThermostat at the beginning and end of each time step
- * to update the Nose-Hoover chain.
- */
-class NoseHooverChainKernel : public KernelImpl {
-public:
-    static std::string Name() {
-        return "NoseHooverChain";
-    }
-    NoseHooverChainKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
-    }
-    /**
-     * Initialize the kernel.
-     */
-    virtual void initialize() = 0;
-    /**
-     * Execute the kernel that propagates the Nose Hoover chain and determines the velocity scale factor.
-     * 
-     * @param context  the context in which to execute this kernel
-     * @param noseHooverChain the object describing the chain to be propagated.
-     * @param kineticEnergy the {center of mass, relative} kineticEnergies of the particles being thermostated by this chain.
-     * @param timeStep the time step used by the integrator.
-     * @return the velocity scale factor to apply to the particles associated with this heat bath.
-     */
-    virtual std::pair<double, double> propagateChain(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> kineticEnergy, double timeStep) = 0;
-    /**
-     * Execute the kernal that computes the total (kinetic + potential) heat bath energy.
-     *
-     * @param context the context in which to execute this kernel
-     * @param noseHooverChain the chain whose energy is to be determined.
-     * @return the total heat bath energy.
-     */
-    virtual double computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain) = 0;
-    /**
-     * Execute the kernel that computes the kinetic energy for a subset of atoms,
-     * or the relative kinetic energy of Drude particles with respect to their parent atoms
-     *
-     * @param context the context in which to execute this kernel
-     * @param noseHooverChain the chain whose energy is to be determined.
-     * @param downloadValue whether the computed value should be downloaded and returned.
-     */
-    virtual std::pair<double, double> computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue) = 0;
-    /**
-     * Execute the kernel that scales the velocities of particles associated with a nose hoover chain
-     *
-     * @param context the context in which to execute this kernel
-     * @param noseHooverChain the chain whose energy is to be determined.
-     * @param scaleFactor the multiplicative factor by which {absolute, relative} velocities are scaled.
-     */
-    virtual void scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactor) = 0;
-};
-
 /**
  * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
  */

openmmapi/include/openmm/NoseHooverChain.h

@@ -72,7 +72,7 @@ public:
      *                                   interact with this chain
      * @param chainLength                the length of (number of particles in) this heat bath
      * @param numMTS                     the number of multi time steps used to propagate this chain
-     * @param numYoshidaSuzuki           the number of Yoshida Suzuki steps used to propagate this chain (1, 3, or 5).
+     * @param numYoshidaSuzuki           the number of Yoshida Suzuki steps used to propagate this chain (1, 3, 5, or 7).
      * @param chainID                    the chain id used to distinguish this Nose-Hoover chain from others that may
      *                                          be used to control a different set of particles, e.g. for Drude oscillators
      * @param thermostatedAtoms                 the list of atoms to be handled by this thermostat

openmmapi/include/openmm/NoseHooverIntegrator.h

@@ -45,13 +45,14 @@ namespace OpenMM {
 class System;
 /**
  * This is an Integrator which simulates a System using one or more Nose Hoover chain
- * thermostats, using the velocity Verlet propagation algorithm.
+ * thermostats, using the "middle" leapfrog propagation algorithm described in
+ * J. Phys. Chem. A 2019, 123, 6056-6079.
  */
 
 class OPENMM_EXPORT NoseHooverIntegrator : public Integrator {
 public:
     /**
-     * Create a NoseHooverIntegrator.  This version creates a bare velocity Verlet integrator
+     * Create a NoseHooverIntegrator.  This version creates a bare leapfrog integrator
      * with no thermostats; any thermostats should be added by calling addThermostat.
      * 
      * @param stepSize the step size with which to integrate the system (in picoseconds)
@@ -66,10 +67,10 @@ public:
      * @param chainLength the number of beads in the Nose-Hoover chain.
      * @param numMTS the number of step in the  multiple time step chain propagation algorithm.
      * @param numYoshidaSuzuki the number of terms in the Yoshida-Suzuki multi time step decomposition
-     *        used in the chain propagation algorithm (must be 1, 3, or 5).
+     *        used in the chain propagation algorithm (must be 1, 3, 5, or 7).
      */
     explicit NoseHooverIntegrator(double temperature, double collisionFrequency, double stepSize,
-                                  int chainLength = 3, int numMTS = 3, int numYoshidaSuzuki = 3);
+                                  int chainLength = 3, int numMTS = 3, int numYoshidaSuzuki = 7);
 
     virtual ~NoseHooverIntegrator();
    /**
@@ -86,7 +87,7 @@ public:
      * @param chainLength the number of beads in the Nose-Hoover chain
      * @param numMTS the number of step in the  multiple time step chain propagation algorithm.
      * @param numYoshidaSuzuki the number of terms in the Yoshida-Suzuki multi time step decomposition
-     *        used in the chain propagation algorithm (must be 1, 3, or 5).
+     *        used in the chain propagation algorithm (must be 1, 3, 5, or 7).
      */
      int addThermostat(double temperature, double collisionFrequency,
                        int chainLength, int numMTS, int numYoshidaSuzuki);
@@ -110,13 +111,13 @@ public:
      * @param chainLength the number of beads in the Nose-Hoover chain.
      * @param numMTS the number of step in the  multiple time step chain propagation algorithm.
      * @param numYoshidaSuzuki the number of terms in the Yoshida-Suzuki multi time step decomposition
-     *        used in the chain propagation algorithm (must be 1, 3, or 5).
+     *        used in the chain propagation algorithm (must be 1, 3, 5, or 7).
      */
      int addSubsystemThermostat(const std::vector<int>& thermostatedParticles,
                                 const std::vector< std::pair< int, int> >& thermostatedPairs,
                                 double temperature, double collisionFrequency, double relativeTemperature,
                                 double relativeCollisionFrequency,
-                                int chainLength = 3, int numMTS = 3, int numYoshidaSuzuki = 3);
+                                int chainLength = 3, int numMTS = 3, int numYoshidaSuzuki = 7);
     /**
      * Get the temperature of the i-th chain for controling absolute particle motion (in Kelvin).
      * 
@@ -232,15 +233,6 @@ public:
      */
     const std::vector<std::tuple<int, int, double> > & getAllThermostatedPairs() const { return allPairs; }
 protected:
-   /**
-     * Advance any Nose-Hoover chains associated with this integrator and determine
-     * scale factor for the velocities.
-     * 
-     * @param kineticEnergy  the {absolute, relative} kinetic energies of the system that the chain is thermostating
-     * @param chainID        id of the Nose-Hoover-Chain
-     * @return the scale factor to be applied to the velocities of the particles thermostated by the chain.
-     */
-    std::pair<double, double> propagateChain(std::pair<double, double> kineticEnergy, int chainID=0);
     /**
      * This will be called by the Context when it is created.  It informs the Integrator
      * of what context it will be integrating, and gives it a chance to do any necessary initialization.
@@ -273,7 +265,7 @@ protected:
     std::vector<int> allAtoms;
     std::vector<std::tuple<int, int, double> > allPairs;
     bool forcesAreValid;
-    Kernel vvKernel, nhcKernel;
+    Kernel kernel;
     bool hasSubsystemThermostats_;
     double maxPairDistance_;
 };

openmmapi/src/NoseHooverChain.cpp

@@ -55,7 +55,10 @@ std::vector<double> NoseHooverChain::getYoshidaSuzukiWeights() const {
         case 5:
             return {0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065,
                     0.2967324292201065};
+        case 7:
+            return {0.784513610477560, 0.235573213359357, -1.17767998417887, 1.31518632068391,
+                    -1.17767998417887, 0.235573213359357, 0.784513610477560};
         default:
-            throw OpenMMException("The number of Yoshida-Suzuki weights must be 1,3, or 5.");
+            throw OpenMMException("The number of Yoshida-Suzuki weights must be 1,3,5, or 7.");
     }
 }

openmmapi/src/NoseHooverIntegrator.cpp

@@ -68,11 +68,6 @@ NoseHooverIntegrator::NoseHooverIntegrator(double temperature, double collisionF
 
 NoseHooverIntegrator::~NoseHooverIntegrator() {}
 
-std::pair<double, double> NoseHooverIntegrator::propagateChain(std::pair<double, double> kineticEnergy, int chainID) {
-    return nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, noseHooverChains.at(chainID), kineticEnergy, getStepSize());
-}
-
-
 int NoseHooverIntegrator::addThermostat(double temperature, double collisionFrequency,
                                         int chainLength, int numMTS, int numYoshidaSuzuki) {
     hasSubsystemThermostats_ = false;
@@ -271,10 +266,10 @@ double NoseHooverIntegrator::computeKineticEnergy() {
     double kE = 0.0;
     if(noseHooverChains.size() > 0) {
         for (const auto &nhc: noseHooverChains){
-            kE += nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, true).first;
+            kE += kernel.getAs<IntegrateNoseHooverStepKernel>().computeMaskedKineticEnergy(*context, nhc, true).first;
         }
     } else {
-        kE = vvKernel.getAs<IntegrateVelocityVerletStepKernel>().computeKineticEnergy(*context, *this);
+        kE = kernel.getAs<IntegrateNoseHooverStepKernel>().computeKineticEnergy(*context, *this);
     }
     return kE;
 }
@@ -287,7 +282,7 @@ double NoseHooverIntegrator::computeHeatBathEnergy() {
     double energy = 0;
     for(auto &nhc : noseHooverChains) {
         if (context && (nhc.getNumDegreesOfFreedom() > 0)) {
-            energy += nhcKernel.getAs<NoseHooverChainKernel>().computeHeatBathEnergy(*context, nhc);
+            energy += kernel.getAs<IntegrateNoseHooverStepKernel>().computeHeatBathEnergy(*context, nhc);
         }
     }
     return energy;
@@ -300,10 +295,8 @@ void NoseHooverIntegrator::initialize(ContextImpl& contextRef) {
     context = &contextRef;
     const System& system = context->getSystem();
     owner = &contextRef.getOwner();
-    vvKernel = context->getPlatform().createKernel(IntegrateVelocityVerletStepKernel::Name(), contextRef);
-    vvKernel.getAs<IntegrateVelocityVerletStepKernel>().initialize(contextRef.getSystem(), *this);
-    nhcKernel = context->getPlatform().createKernel(NoseHooverChainKernel::Name(), contextRef);
-    nhcKernel.getAs<NoseHooverChainKernel>().initialize();
+    kernel = context->getPlatform().createKernel(IntegrateNoseHooverStepKernel::Name(), contextRef);
+    kernel.getAs<IntegrateNoseHooverStepKernel>().initialize(contextRef.getSystem(), *this);
     forcesAreValid = false;
 
     // check for drude particles and build the Nose-Hoover Chains
@@ -329,34 +322,22 @@ void NoseHooverIntegrator::initialize(ContextImpl& contextRef) {
 }
 
 void NoseHooverIntegrator::cleanup() {
-    vvKernel = Kernel();
-    nhcKernel = Kernel();
+    kernel = Kernel();
 }
 
 vector<string> NoseHooverIntegrator::getKernelNames() {
     std::vector<std::string> names;
-    names.push_back(NoseHooverChainKernel::Name());
-    names.push_back(IntegrateVelocityVerletStepKernel::Name());
+    names.push_back(IntegrateNoseHooverStepKernel::Name());
     return names;
 }
 
 void NoseHooverIntegrator::step(int steps) {
     if (context == NULL)
         throw OpenMMException("This Integrator is not bound to a context!");
-    std::pair<double, double> scale, kineticEnergy;
     for (int i = 0; i < steps; ++i) {
         if(context->updateContextState())
             forcesAreValid = false;
-        for(auto &nhc : noseHooverChains) {
-            kineticEnergy = nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, false);
-            scale = nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, nhc, kineticEnergy, getStepSize());
-            nhcKernel.getAs<NoseHooverChainKernel>().scaleVelocities(*context, nhc, scale);
-        }
-        vvKernel.getAs<IntegrateVelocityVerletStepKernel>().execute(*context, *this, forcesAreValid);
-        for(auto &nhc : noseHooverChains) {
-            kineticEnergy = nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, false);
-            scale = nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, nhc, kineticEnergy, getStepSize());
-            nhcKernel.getAs<NoseHooverChainKernel>().scaleVelocities(*context, nhc, scale);
-        }
+        context->calcForcesAndEnergy(true, false);
+        kernel.getAs<IntegrateNoseHooverStepKernel>().execute(*context, *this, forcesAreValid);
     }
 }

platforms/common/include/openmm/common/CommonKernels.h

@@ -944,11 +944,13 @@ private:
 /*
  * This kernel is invoked by NoseHooverIntegrator to take one time step.
  */
-class CommonIntegrateVelocityVerletStepKernel : public IntegrateVelocityVerletStepKernel {
+class CommonIntegrateNoseHooverStepKernel : public IntegrateNoseHooverStepKernel {
 public:
-    CommonIntegrateVelocityVerletStepKernel(std::string name, const Platform& platform, ComputeContext& cc) :
-                                  IntegrateVelocityVerletStepKernel(name, platform), cc(cc), hasInitializedKernels(false) { }
-    ~CommonIntegrateVelocityVerletStepKernel() {}
+    CommonIntegrateNoseHooverStepKernel(std::string name, const Platform& platform, ComputeContext& cc) :
+                                  IntegrateNoseHooverStepKernel(name, platform), cc(cc), hasInitializedKernels(false),
+                                  hasInitializedKineticEnergyKernel(false), hasInitializedHeatBathEnergyKernel(false),
+                                  hasInitializedScaleVelocitiesKernel(false), hasInitializedPropagateKernel(false) {}
+    ~CommonIntegrateNoseHooverStepKernel() {}
     /**
      * Initialize the kernel.
      *
@@ -972,39 +974,16 @@ public:
      * @param integrator the NoseHooverIntegrator this kernel is being used for
      */
     double computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator);
-private:
-    ComputeContext& cc;
-    float prevMaxPairDistance;
-    ComputeArray maxPairDistanceBuffer, pairListBuffer, atomListBuffer, pairTemperatureBuffer;
-    ComputeKernel kernel1, kernel2, kernel3, kernelHardWall;
-    bool hasInitializedKernels;
-};
-
-/**
- * This kernel is invoked by NoseHooverChain at the start of each time step to adjust the thermostat
- * and update the associated particle velocities.
- */
-class CommonNoseHooverChainKernel : public NoseHooverChainKernel {
-public:
-    CommonNoseHooverChainKernel(std::string name, const Platform& platform, ComputeContext& cc) :
-                 NoseHooverChainKernel(name, platform), cc(cc), hasInitializedPropagateKernel(false),
-                 hasInitializedKineticEnergyKernel(false), hasInitializedHeatBathEnergyKernel(false),
-                 hasInitializedScaleVelocitiesKernel(false) {}
-    ~CommonNoseHooverChainKernel() {}
-    /**
-     * Initialize the kernel.
-     */
-    void initialize();
     /**
      * Execute the kernel that propagates the Nose Hoover chain and determines the velocity scale factor.
      * 
      * @param context  the context in which to execute this kernel
      * @param noseHooverChain the object describing the chain to be propagated.
-     * @param kineticEnergies the {absolute, relative} kineticEnergy of the particles being thermostated by this chain.
+     * @param kineticEnergy the {center of mass, relative} kineticEnergies of the particles being thermostated by this chain.
      * @param timeStep the time step used by the integrator.
-     * @return the {absolute, relative} velocity scale factor to apply to the particles associated with this heat bath.
+     * @return the velocity scale factor to apply to the particles associated with this heat bath.
      */
-    std::pair<double, double> propagateChain(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> kineticEnergies, double timeStep);
+    std::pair<double, double> propagateChain(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> kineticEnergy, double timeStep);
     /**
      * Execute the kernal that computes the total (kinetic + potential) heat bath energy.
      *
@@ -1012,7 +991,7 @@ public:
      * @param noseHooverChain the chain whose energy is to be determined.
      * @return the total heat bath energy.
      */
-    double computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &nhc);
+    double computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain);
     /**
      * Execute the kernel that computes the kinetic energy for a subset of atoms,
      * or the relative kinetic energy of Drude particles with respect to their parent atoms
@@ -1020,22 +999,28 @@ public:
      * @param context the context in which to execute this kernel
      * @param noseHooverChain the chain whose energy is to be determined.
      * @param downloadValue whether the computed value should be downloaded and returned.
-     *
      */
-    std::pair<double,double> computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue);
-
+    std::pair<double, double> computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue);
     /**
      * Execute the kernel that scales the velocities of particles associated with a nose hoover chain
      *
      * @param context the context in which to execute this kernel
      * @param noseHooverChain the chain whose energy is to be determined.
-     * @param scaleFactors the {absolute, relative} multiplicative factor by which velocities are scaled.
+     * @param scaleFactor the multiplicative factor by which {absolute, relative} velocities are scaled.
      */
-    void scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactors);
-
+    void scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactor);
 private:
-    int sumWorkGroupSize;
     ComputeContext& cc;
+    float prevMaxPairDistance;
+    ComputeArray maxPairDistanceBuffer, pairListBuffer, atomListBuffer, pairTemperatureBuffer, oldDelta;
+    ComputeKernel kernel1, kernel2, kernel3, kernel4, kernelHardWall;
+    bool hasInitializedKernels;
+    ComputeKernel reduceEnergyKernel;
+    ComputeKernel computeHeatBathEnergyKernel;
+    ComputeKernel computeAtomsKineticEnergyKernel;
+    ComputeKernel computePairsKineticEnergyKernel;
+    ComputeKernel scaleAtomsVelocitiesKernel;
+    ComputeKernel scalePairsVelocitiesKernel;
     ComputeArray energyBuffer, scaleFactorBuffer, kineticEnergyBuffer, chainMasses, chainForces, heatBathEnergy;
     std::map<int, ComputeArray> atomlists, pairlists;
     std::map<int, ComputeKernel> propagateKernels;
@@ -1043,12 +1028,6 @@ private:
     bool hasInitializedKineticEnergyKernel;
     bool hasInitializedHeatBathEnergyKernel;
     bool hasInitializedScaleVelocitiesKernel;
-    ComputeKernel reduceEnergyKernel;
-    ComputeKernel computeHeatBathEnergyKernel;
-    ComputeKernel computeAtomsKineticEnergyKernel;
-    ComputeKernel computePairsKineticEnergyKernel;
-    ComputeKernel scaleAtomsVelocitiesKernel;
-    ComputeKernel scalePairsVelocitiesKernel;
 };
 
 /**

platforms/common/src/CommonKernels.cpp

@@ -5640,20 +5640,68 @@ double CommonIntegrateLangevinMiddleStepKernel::computeKineticEnergy(ContextImpl
     return cc.getIntegrationUtilities().computeKineticEnergy(0.0);
 }
 
-void CommonIntegrateVelocityVerletStepKernel::initialize(const System& system, const NoseHooverIntegrator& integrator) {
+void CommonIntegrateNoseHooverStepKernel::initialize(const System& system, const NoseHooverIntegrator& integrator) {
     cc.initializeContexts();
+    bool useDouble = cc.getUseDoublePrecision() || cc.getUseMixedPrecision();
     map<string, string> defines;
     defines["BOLTZ"] = cc.doubleToString(BOLTZ);
-    ComputeProgram program = cc.compileProgram(CommonKernelSources::velocityVerlet, defines);
-    kernel1 = program->createKernel("integrateVelocityVerletPart1");
-    kernel2 = program->createKernel("integrateVelocityVerletPart2");
-    kernel3 = program->createKernel("integrateVelocityVerletPart3");
-    kernelHardWall = program->createKernel("integrateVelocityVerletHardWall");
+    ComputeProgram program = cc.compileProgram(CommonKernelSources::noseHooverIntegrator, defines);
+    kernel1 = program->createKernel("integrateNoseHooverMiddlePart1");
+    kernel2 = program->createKernel("integrateNoseHooverMiddlePart2");
+    kernel3 = program->createKernel("integrateNoseHooverMiddlePart3");
+    kernel4 = program->createKernel("integrateNoseHooverMiddlePart4");
+    if (useDouble) {
+        oldDelta.initialize<mm_double4>(cc, cc.getPaddedNumAtoms(), "oldDelta");
+    } else {
+        oldDelta.initialize<mm_float4>(cc, cc.getPaddedNumAtoms(), "oldDelta");
+    }
+    kernelHardWall = program->createKernel("integrateNoseHooverHardWall");
     prevMaxPairDistance = -1.0f;
     maxPairDistanceBuffer.initialize<float>(cc, 1, "maxPairDistanceBuffer");
+
+    int workGroupSize = std::min(cc.getMaxThreadBlockSize(), 512);
+    defines["WORK_GROUP_SIZE"] = std::to_string(workGroupSize);
+
+    defines["BEGIN_YS_LOOP"] = "const real arr[1] = {1.0};"
+                               "for(int i=0;i<1;++i) {"
+                               "const real ys = arr[i];";
+    defines["END_YS_LOOP"] = "}";
+    program = cc.compileProgram(CommonKernelSources::noseHooverChain, defines);
+    propagateKernels[1] = program->createKernel("propagateNoseHooverChain");
+
+    defines["BEGIN_YS_LOOP"] = "const real arr[3] = {0.828981543588751, -0.657963087177502, 0.828981543588751};"
+                               "for(int i=0;i<3;++i) {"
+                               "const real ys = arr[i];";
+    program = cc.compileProgram(CommonKernelSources::noseHooverChain, defines);
+    propagateKernels[3] = program->createKernel("propagateNoseHooverChain");
+
+    defines["BEGIN_YS_LOOP"] = "const real arr[5] = {0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065, 0.2967324292201065};"
+                               "for(int i=0;i<5;++i) {"
+                               "const real ys = arr[i];";
+    program = cc.compileProgram(CommonKernelSources::noseHooverChain, defines);
+    propagateKernels[5] = program->createKernel("propagateNoseHooverChain");
+
+    defines["BEGIN_YS_LOOP"] = "const real arr[7] = {0.784513610477560, 0.235573213359357, -1.17767998417887, 1.31518632068391,-1.17767998417887, 0.235573213359357, 0.784513610477560};"
+                               "for(int i=0;i<7;++i) {"
+                               "const real ys = arr[i];";
+    program = cc.compileProgram(CommonKernelSources::noseHooverChain, defines);
+    propagateKernels[7] = program->createKernel("propagateNoseHooverChain");
+    program = cc.compileProgram(CommonKernelSources::noseHooverChain, defines);
+    reduceEnergyKernel = program->createKernel("reduceEnergyPair");
+
+    computeHeatBathEnergyKernel = program->createKernel("computeHeatBathEnergy");
+    computeAtomsKineticEnergyKernel = program->createKernel("computeAtomsKineticEnergy");
+    computePairsKineticEnergyKernel = program->createKernel("computePairsKineticEnergy");
+    scaleAtomsVelocitiesKernel = program->createKernel("scaleAtomsVelocities");
+    scalePairsVelocitiesKernel = program->createKernel("scalePairsVelocities");
+    int energyBufferSize = cc.getEnergyBuffer().getSize();
+    if (cc.getUseDoublePrecision() || cc.getUseMixedPrecision())
+        energyBuffer.initialize<mm_double2>(cc, energyBufferSize, "energyBuffer");
+    else
+        energyBuffer.initialize<mm_float2>(cc, energyBufferSize, "energyBuffer");
 }
 
-void CommonIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) {
+void CommonIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) {
     IntegrationUtilities& integration = cc.getIntegrationUtilities();
     int paddedNumAtoms = cc.getPaddedNumAtoms();
     double dt = integrator.getStepSize();
@@ -5700,32 +5748,39 @@ void CommonIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, cons
         pairListBuffer.upload(tmp);
         pairTemperatureBuffer.upload(tmp2);
     }
-
+    int totalAtoms = cc.getNumAtoms();
     if (!hasInitializedKernels) {
         hasInitializedKernels = true;
         kernel1->addArg(numAtoms);
         kernel1->addArg(numPairs);
         kernel1->addArg(paddedNumAtoms);
-        kernel1->addArg(cc.getIntegrationUtilities().getStepSize());
-        kernel1->addArg(cc.getPosq());
         kernel1->addArg(cc.getVelm());
         kernel1->addArg(cc.getLongForceBuffer());
-        kernel1->addArg(integration.getPosDelta());
+        kernel1->addArg(integration.getStepSize());
         kernel1->addArg(numAtoms > 0 ? atomListBuffer : cc.getEnergyBuffer()); // The array is not used if num == 0
         kernel1->addArg(numPairs > 0 ? pairListBuffer : cc.getEnergyBuffer()); // The array is not used if num == 0
-        if (cc.getUseMixedPrecision())
-            kernel1->addArg(cc.getPosqCorrection());
-        kernel2->addArg(numParticles);
-        kernel2->addArg(cc.getIntegrationUtilities().getStepSize());
-        kernel2->addArg(cc.getPosq());
+        kernel2->addArg(totalAtoms);
         kernel2->addArg(cc.getVelm());
         kernel2->addArg(integration.getPosDelta());
+        kernel2->addArg(oldDelta);
+        kernel2->addArg(integration.getStepSize());
+        kernel3->addArg(totalAtoms);
+        kernel3->addArg(cc.getVelm());
+        kernel3->addArg(integration.getPosDelta());
+        kernel3->addArg(oldDelta);
+        kernel3->addArg(integration.getStepSize());
+        kernel4->addArg(totalAtoms);
+        kernel4->addArg(cc.getPosq());
+        kernel4->addArg(cc.getVelm());
+        kernel4->addArg(integration.getPosDelta());
+        kernel4->addArg(oldDelta);
+        kernel4->addArg(integration.getStepSize());
         if (cc.getUseMixedPrecision())
-            kernel2->addArg(cc.getPosqCorrection());
+            kernel4->addArg(cc.getPosqCorrection());
         if (numPairs > 0) {
             kernelHardWall->addArg(numPairs);
             kernelHardWall->addArg(maxPairDistanceBuffer);
-            kernelHardWall->addArg(cc.getIntegrationUtilities().getStepSize());
+            kernelHardWall->addArg(integration.getStepSize());
             kernelHardWall->addArg(cc.getPosq());
             kernelHardWall->addArg(cc.getVelm());
             kernelHardWall->addArg(pairListBuffer);
@@ -5733,78 +5788,51 @@ void CommonIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, cons
             if (cc.getUseMixedPrecision())
                 kernelHardWall->addArg(cc.getPosqCorrection());
         }
-        kernel3->addArg(numAtoms);
-        kernel3->addArg(numPairs);
-        kernel3->addArg(paddedNumAtoms);
-        kernel3->addArg(cc.getIntegrationUtilities().getStepSize());
-        kernel3->addArg(cc.getPosq());
-        kernel3->addArg(cc.getVelm());
-        kernel3->addArg(cc.getLongForceBuffer());
-        kernel3->addArg(integration.getPosDelta());
-        kernel3->addArg(numAtoms > 0 ? atomListBuffer : cc.getEnergyBuffer()); // The array is not used if num == 0
-        kernel3->addArg(numPairs > 0 ? pairListBuffer : cc.getEnergyBuffer()); // The array is not used if num == 0
-        if (cc.getUseMixedPrecision())
-            kernel3->addArg(cc.getPosqCorrection());
     }
 
     /*
-     * Carry out the integration
+     * Carry out the BAOAB integration, also known as the "LF-middle" scheme
+     * c.f. J. Phys. Chem. A 2019, 123, 6056−6079
      */
-
-    // Advance the velocities a half step
+    // B
     kernel1->execute(std::max(numAtoms, numPairs));
-    integration.applyConstraints(integrator.getConstraintTolerance());
-    // Advance particle positions a full step
+    integration.applyVelocityConstraints(integrator.getConstraintTolerance());
+    // A
     kernel2->execute(numParticles);
+    // O
+    int numChains = integrator.getNumThermostats();
+    for(int chain = 0; chain < numChains; ++chain) {
+        const auto &thermostatChain = integrator.getThermostat(chain);
+        auto KEs = computeMaskedKineticEnergy(context, thermostatChain, false);
+        auto scaleFactors = propagateChain(context, thermostatChain, KEs, dt);
+        scaleVelocities(context, thermostatChain, scaleFactors);
+    }
+    // A
+    kernel3->execute(numParticles);
+    integration.applyConstraints(integrator.getConstraintTolerance());
+    // B
+    kernel4->execute(numAtoms);
     // Make sure any Drude-like particles have not wandered too far from home
     if (numPairs > 0) kernelHardWall->execute(numPairs);
     integration.computeVirtualSites();
-    context.calcForcesAndEnergy(true, false);
-    forcesAreValid = true;
-    // Update velocities another half step
-    kernel3->execute(std::max(numAtoms, numPairs));
-    integration.applyVelocityConstraints(integrator.getConstraintTolerance());
 
+    // Update the time and step count.
     cc.setTime(cc.getTime()+dt);
     cc.setStepCount(cc.getStepCount()+1);
     cc.reorderAtoms();
+
+    // Reduce UI lag.
+#ifdef WIN32
+    cc.flushQueue();
+#endif
 }
 
-double CommonIntegrateVelocityVerletStepKernel::computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator) {
+double CommonIntegrateNoseHooverStepKernel::computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator) {
     return cc.getIntegrationUtilities().computeKineticEnergy(0);
 }
 
-void CommonNoseHooverChainKernel::initialize() {
-    bool useDouble = cc.getUseDoublePrecision() || cc.getUseMixedPrecision();
-    map<string, string> defines;
-    int workGroupSize = std::min(cc.getMaxThreadBlockSize(), 512);
-    defines["WORK_GROUP_SIZE"] = std::to_string(workGroupSize);
-    defines["BEGIN_YS_LOOP"] = "const real arr[1] = {1.0}; for(int i=0;i<1;++i) { const real ys = arr[i];";
-    defines["END_YS_LOOP"] = "}";
-    ComputeProgram program = cc.compileProgram(CommonKernelSources::noseHooverChain, defines);
-    propagateKernels[1] = program->createKernel("propagateNoseHooverChain");
-    defines["BEGIN_YS_LOOP"] = "const real arr[3] = {0.828981543588751, -0.657963087177502, 0.828981543588751}; for(int i=0;i<3;++i) { const real ys = arr[i];";
-    program = cc.compileProgram(CommonKernelSources::noseHooverChain, defines);
-    propagateKernels[3] = program->createKernel("propagateNoseHooverChain");
-    defines["BEGIN_YS_LOOP"] = "const real arr[5] = {0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065, 0.2967324292201065}; for(int i=0;i<5;++i) { const real ys = arr[i];";
-    program = cc.compileProgram(CommonKernelSources::noseHooverChain, defines);
-    propagateKernels[5] = program->createKernel("propagateNoseHooverChain");
-    program = cc.compileProgram(CommonKernelSources::noseHooverChain, defines);
-    reduceEnergyKernel = program->createKernel("reduceEnergyPair");
-
-    computeHeatBathEnergyKernel = program->createKernel("computeHeatBathEnergy");
-    computeAtomsKineticEnergyKernel = program->createKernel("computeAtomsKineticEnergy");
-    computePairsKineticEnergyKernel = program->createKernel("computePairsKineticEnergy");
-    scaleAtomsVelocitiesKernel = program->createKernel("scaleAtomsVelocities");
-    scalePairsVelocitiesKernel = program->createKernel("scalePairsVelocities");
-    int energyBufferSize = cc.getEnergyBuffer().getSize();
-    if (cc.getUseDoublePrecision() || cc.getUseMixedPrecision())
-        energyBuffer.initialize<mm_double2>(cc, energyBufferSize, "energyBuffer");
-    else
-        energyBuffer.initialize<mm_float2>(cc, energyBufferSize, "energyBuffer");
-}
 
-std::pair<double, double> CommonNoseHooverChainKernel::propagateChain(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> kineticEnergies, double timeStep) {
+std::pair<double, double> CommonIntegrateNoseHooverStepKernel::propagateChain(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> kineticEnergies, double timeStep) {
     bool useDouble = cc.getUseDoublePrecision() || cc.getUseMixedPrecision();
     int chainID = nhc.getChainID();
     int nAtoms = nhc.getThermostatedAtoms().size();
@@ -5813,8 +5841,8 @@ std::pair<double, double> CommonNoseHooverChainKernel::propagateChain(ContextImp
     int numYS = nhc.getNumYoshidaSuzukiTimeSteps();
     int numMTS = nhc.getNumMultiTimeSteps();
 
-    if (numYS != 1 && numYS != 3 && numYS != 5) {
-        throw OpenMMException("Number of Yoshida Suzuki time steps has to be 1, 3, or 5.");
+    if (numYS != 1 && numYS != 3 && numYS != 5 && numYS != 7) {
+        throw OpenMMException("Number of Yoshida Suzuki time steps has to be 1, 3, 5, or 7.");
     }
 
     auto & chainState = cc.getIntegrationUtilities().getNoseHooverChainState();
@@ -5973,7 +6001,7 @@ std::pair<double, double> CommonNoseHooverChainKernel::propagateChain(ContextImp
     return {0, 0};
 }
 
-double CommonNoseHooverChainKernel::computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &nhc) {
+double CommonIntegrateNoseHooverStepKernel::computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &nhc) {
 
     bool useDouble = cc.getUseDoublePrecision() || cc.getUseMixedPrecision();
 
@@ -6058,7 +6086,7 @@ double CommonNoseHooverChainKernel::computeHeatBathEnergy(ContextImpl& context,
         return *((float*) pinnedBuffer);
 }
 
-std::pair<double, double> CommonNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &nhc, bool downloadValue) {
+std::pair<double, double> CommonIntegrateNoseHooverStepKernel::computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &nhc, bool downloadValue) {
 
     bool useDouble = cc.getUseDoublePrecision() || cc.getUseMixedPrecision();
 
@@ -6153,7 +6181,7 @@ std::pair<double, double> CommonNoseHooverChainKernel::computeMaskedKineticEnerg
     return KEs;
 }
 
-void CommonNoseHooverChainKernel::scaleVelocities(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> scaleFactor) {
+void CommonIntegrateNoseHooverStepKernel::scaleVelocities(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> scaleFactor) {
     // For now we assume that the atoms and pairs info is valid, because compute{Atoms|Pairs}KineticEnergy must have been
     // called before this kernel.  If that ever ceases to be true, some sanity checks are needed here.
 

platforms/common/src/kernels/langevinMiddle.cc

@@ -28,7 +28,7 @@ KERNEL void integrateLangevinMiddlePart2(int numAtoms, GLOBAL mixed4* RESTRICT v
         ) {
     mixed vscale = paramBuffer[VelScale];
     mixed noisescale = paramBuffer[NoiseScale];
-    mixed halfdt = 0.5*dt[0].y;
+    mixed halfdt = 0.5f*dt[0].y;
     int index = GLOBAL_ID;
     randomIndex += index;
     while (index < numAtoms) {

platforms/common/src/kernels/noseHooverChain.cc

+// Propagates a Nose Hoover chain a full timestep
 KERNEL void propagateNoseHooverChain(GLOBAL mixed2* RESTRICT chainData, GLOBAL const mixed2 * RESTRICT energySum, GLOBAL mixed2* RESTRICT scaleFactor,
                                      GLOBAL mixed* RESTRICT chainMasses, GLOBAL mixed* RESTRICT chainForces, int chainType, int chainLength, int numMTS,
                                      int numDOFs, float timeStep, mixed kT, float frequency){
@@ -15,30 +16,29 @@ KERNEL void propagateNoseHooverChain(GLOBAL mixed2* RESTRICT chainData, GLOBAL c
     for (int mts = 0; mts < numMTS; ++mts) {
         BEGIN_YS_LOOP
             mixed wdt = ys * timeOverMTS;
-            chainData[chainLength-1].y += 0.25f * wdt * chainForces[chainLength-1];
+            chainData[chainLength-1].y += 0.5f * wdt * chainForces[chainLength-1];
             for (int bead = chainLength - 2; bead >= 0; --bead) {
-                mixed aa = EXP(-0.125f * wdt * chainData[bead + 1].y);
-                chainData[bead].y = aa * (chainData[bead].y * aa + 0.25f * wdt * chainForces[bead]);
+                mixed aa = EXP(-0.25f * wdt * chainData[bead + 1].y);
+                chainData[bead].y = aa * (chainData[bead].y * aa + 0.5f * wdt * chainForces[bead]);
             }
             // update particle velocities
-            mixed aa = EXP(-0.5f * wdt * chainData[0].y);
+            mixed aa = EXP(-wdt * chainData[0].y);
             scale *= aa;
             // update the thermostat positions
             for (int bead = 0; bead < chainLength; ++bead) {
-                chainData[bead].x += 0.5f * chainData[bead].y * wdt;
+                chainData[bead].x += chainData[bead].y * wdt;
             }
             // update the forces
             chainForces[0] = (scale * scale * KE2 - numDOFs * kT) / chainMasses[0];
             // update thermostat velocities
             for (int bead = 0; bead < chainLength - 1; ++bead) {
-                mixed aa = EXP(-0.125f * wdt * chainData[bead + 1].y);
-                chainData[bead].y = aa * (aa * chainData[bead].y + 0.25f * wdt * chainForces[bead]);
+                mixed aa = EXP(-0.25f * wdt * chainData[bead + 1].y);
+                chainData[bead].y = aa * (aa * chainData[bead].y + 0.5f * wdt * chainForces[bead]);
                 chainForces[bead + 1] = (chainMasses[bead] * chainData[bead].y * chainData[bead].y - kT) / chainMasses[bead + 1];
             }
-            chainData[chainLength-1].y += 0.25f * wdt * chainForces[chainLength-1];
+            chainData[chainLength-1].y += 0.5f * wdt * chainForces[chainLength-1];
         END_YS_LOOP
     } // MTS loop
-
     if (chainType == 0) {
         scaleFactor[0].x = scale;
     } else {

platforms/common/src/kernels/velocityVerlet.cc → platforms/common/src/kernels/noseHooverIntegrator.cc

 /**
- * Perform the first step of Velocity Verlet integration.
+ * Perform the first part of integration: velocity step.
  */
-
-KERNEL void integrateVelocityVerletPart1(int numAtoms, int numPairs, int paddedNumAtoms, GLOBAL const mixed2* RESTRICT dt, GLOBAL const real4* RESTRICT posq,
-                                         GLOBAL mixed4* RESTRICT velm, GLOBAL const mm_long* RESTRICT force, GLOBAL mixed4* RESTRICT posDelta,
-                                         GLOBAL const int* RESTRICT atomList, GLOBAL const int2* RESTRICT pairList
-#ifdef USE_MIXED_PRECISION
-                                        ,GLOBAL const real4* RESTRICT posqCorrection
-#endif
-    ){
-    const mixed2 stepSize = dt[0];
-    const mixed dtPos = stepSize.y;
-    const mixed dtVel = 0.5f*(stepSize.x+stepSize.y);
-    const mixed scale = 0.5f * dtVel/(mixed) 0x100000000;
+KERNEL void integrateNoseHooverMiddlePart1(int numAtoms, int numPairs, int paddedNumAtoms, GLOBAL mixed4* RESTRICT velm, GLOBAL const mm_long* RESTRICT force,
+        GLOBAL const mixed2* RESTRICT dt, GLOBAL const int* RESTRICT atomList, GLOBAL const int2* RESTRICT pairList) {
+    mixed fscale = dt[0].y/(mixed) 0x100000000;
     int index = GLOBAL_ID;
     while (index < numAtoms) {
         int atom = atomList[index];
         mixed4 velocity = velm[atom];
         if (velocity.w != 0.0) {
-#ifdef USE_MIXED_PRECISION
-            real4 pos1 = posq[atom];
-            real4 pos2 = posqCorrection[atom];
-            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[atom];
-#endif
-            velocity.x += scale*force[atom]*velocity.w;
-            velocity.y += scale*force[atom+paddedNumAtoms]*velocity.w;
-            velocity.z += scale*force[atom+paddedNumAtoms*2]*velocity.w;
-            pos.x = velocity.x*dtPos;
-            pos.y = velocity.y*dtPos;
-            pos.z = velocity.z*dtPos;
-            posDelta[atom] = pos;
+            velocity.x += fscale*force[atom]*velocity.w;
+            velocity.y += fscale*force[atom+paddedNumAtoms]*velocity.w;
+            velocity.z += fscale*force[atom+paddedNumAtoms*2]*velocity.w;
             velm[atom] = velocity;
         }
         index += GLOBAL_SIZE;
@@ -58,12 +38,12 @@ KERNEL void integrateVelocityVerletPart1(int numAtoms, int numPairs, int paddedN
         relVel.z= v2.z - v1.z;
 
         mixed3 comFrc;
-        mixed F1x = scale*force[atom1];
-        mixed F1y = scale*force[atom1+paddedNumAtoms];
-        mixed F1z = scale*force[atom1+paddedNumAtoms*2];
-        mixed F2x = scale*force[atom2];
-        mixed F2y = scale*force[atom2+paddedNumAtoms];
-        mixed F2z = scale*force[atom2+paddedNumAtoms*2];
+        mixed F1x = fscale*force[atom1];
+        mixed F1y = fscale*force[atom1+paddedNumAtoms];
+        mixed F1z = fscale*force[atom1+paddedNumAtoms*2];
+        mixed F2x = fscale*force[atom2];
+        mixed F2y = fscale*force[atom2+paddedNumAtoms];
+        mixed F2z = fscale*force[atom2+paddedNumAtoms*2];
         comFrc.x = F1x + F2x;
         comFrc.y = F1y + F2y;
         comFrc.z = F1z + F2z;
@@ -77,35 +57,16 @@ KERNEL void integrateVelocityVerletPart1(int numAtoms, int numPairs, int paddedN
         relVel.x += relFrc.x * invRedMass;
         relVel.y += relFrc.y * invRedMass;
         relVel.z += relFrc.z * invRedMass;
-#ifdef USE_MIXED_PRECISION
-        real4 posv1 = posq[atom1];
-        real4 posv2 = posq[atom2];
-        real4 posc1 = posqCorrection[atom1];
-        real4 posc2 = posqCorrection[atom2];
-        mixed4 pos1 = make_mixed4(posv1.x+(mixed)posc1.x, posv1.y+(mixed)posc1.y, posv1.z+(mixed)posc1.z, posv1.w);
-        mixed4 pos2 = make_mixed4(posv2.x+(mixed)posc2.x, posv2.y+(mixed)posc2.y, posv2.z+(mixed)posc2.z, posv2.w);
-#else
-        real4 pos1 = posq[atom1];
-        real4 pos2 = posq[atom2];
-#endif
         if (v1.w != 0.0f) {
             v1.x = comVel.x - relVel.x*mass2fract;
             v1.y = comVel.y - relVel.y*mass2fract;
             v1.z = comVel.z - relVel.z*mass2fract;
-            pos1.x = v1.x*dtPos;
-            pos1.y = v1.y*dtPos;
-            pos1.z = v1.z*dtPos;
-            posDelta[atom1] = pos1;
             velm[atom1] = v1;
         }
         if (v2.w != 0.0f) {
             v2.x = comVel.x + relVel.x*mass1fract;
             v2.y = comVel.y + relVel.y*mass1fract;
             v2.z = comVel.z + relVel.z*mass1fract;
-            pos2.x = v2.x*dtPos;
-            pos2.y = v2.y*dtPos;
-            pos2.z = v2.z*dtPos;
-            posDelta[atom2] = pos2;
             velm[atom2] = v2;
         }
         index += GLOBAL_SIZE;
@@ -113,22 +74,60 @@ KERNEL void integrateVelocityVerletPart1(int numAtoms, int numPairs, int paddedN
 }
 
 /**
- * Perform the second step of Velocity Verlet integration.
+ * Perform the second part of integration: position half step
+ */
+KERNEL void integrateNoseHooverMiddlePart2(int numAtoms, GLOBAL mixed4* RESTRICT velm, GLOBAL mixed4* RESTRICT posDelta,
+        GLOBAL mixed4* RESTRICT oldDelta, GLOBAL const mixed2* RESTRICT dt) {
+    mixed halfdt = 0.5f*dt[0].y;
+    int index = GLOBAL_ID;
+    while (index < numAtoms) {
+        mixed4 velocity = velm[index];
+        if (velocity.w != 0.0) {
+            mixed4 delta = make_mixed4(halfdt*velocity.x, halfdt*velocity.y, halfdt*velocity.z, 0);
+            posDelta[index] = delta;
+            oldDelta[index] = delta;
+        }
+        index += GLOBAL_SIZE;
+    }
+}
+
+/**
+ * Perform the third part of integration: another position half step
  */
+KERNEL void integrateNoseHooverMiddlePart3(int numAtoms, GLOBAL mixed4* RESTRICT velm, GLOBAL mixed4* RESTRICT posDelta,
+        GLOBAL mixed4* RESTRICT oldDelta, GLOBAL const mixed2* RESTRICT dt) {
+    mixed halfdt = 0.5f*dt[0].y;
+    int index = GLOBAL_ID;
+    while (index < numAtoms) {
+        mixed4 velocity = velm[index];
+        if (velocity.w != 0.0) {
+            mixed4 delta = make_mixed4(halfdt*velocity.x, halfdt*velocity.y, halfdt*velocity.z, 0);
+            posDelta[index] += delta;
+            oldDelta[index] += delta;
+        }
+        index += GLOBAL_SIZE;
+    }
+}
 
-KERNEL void integrateVelocityVerletPart2(int numAtoms, GLOBAL mixed2* RESTRICT dt, GLOBAL real4* RESTRICT posq, GLOBAL mixed4* RESTRICT velm,
-                                         GLOBAL const mixed4* RESTRICT posDelta
+/**
+ * Perform the fourth part of integration: apply constraint forces to velocities, then record
+ * the constrained positions.
+ */
+KERNEL void integrateNoseHooverMiddlePart4(int numAtoms, GLOBAL real4* RESTRICT posq, GLOBAL mixed4* RESTRICT velm,
+         GLOBAL mixed4* RESTRICT posDelta, GLOBAL mixed4* RESTRICT oldDelta, GLOBAL const mixed2* RESTRICT dt
 #ifdef USE_MIXED_PRECISION
-                                        ,GLOBAL real4* RESTRICT posqCorrection
+        , GLOBAL real4* RESTRICT posqCorrection
 #endif
-    ){
-    mixed2 stepSize = dt[0];
-    int index = GLOBAL_ID;
-    if (index == 0)
-        dt[0].x = stepSize.y;
-    while(index < numAtoms) {
+        ) {
+    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)*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];
             real4 pos2 = posqCorrection[index];
@@ -136,7 +135,6 @@ KERNEL void integrateVelocityVerletPart2(int numAtoms, GLOBAL mixed2* RESTRICT d
 #else
             real4 pos = posq[index];
 #endif
-            mixed4 delta = posDelta[index];
             pos.x += delta.x;
             pos.y += delta.y;
             pos.z += delta.z;
@@ -147,120 +145,10 @@ KERNEL void integrateVelocityVerletPart2(int numAtoms, GLOBAL mixed2* RESTRICT d
             posq[index] = pos;
 #endif
         }
-        index += GLOBAL_SIZE;
     }
 }
 
-/**
- * Perform the third step of Velocity Verlet integration.
- */
-
-KERNEL void integrateVelocityVerletPart3(int numAtoms, int numPairs, int paddedNumAtoms, GLOBAL mixed2* RESTRICT dt, GLOBAL real4* RESTRICT posq,
-                                         GLOBAL mixed4* RESTRICT velm, GLOBAL const mm_long* RESTRICT force, GLOBAL const mixed4* RESTRICT posDelta,
-                                         GLOBAL const int* RESTRICT atomList, GLOBAL const int2* RESTRICT pairList
-#ifdef USE_MIXED_PRECISION
-                                        ,GLOBAL const real4* RESTRICT posqCorrection
-#endif
-    ){
-    mixed2 stepSize = dt[0];
-#ifdef SUPPORTS_DOUBLE_PRECISION
-    double oneOverDt = 1.0/stepSize.y;
-#else
-    float oneOverDt = 1.0f/stepSize.y;
-    float correction = (1.0f-oneOverDt*stepSize.y)/stepSize.y;
-#endif
-    const mixed dtVel = 0.5f*(stepSize.x+stepSize.y);
-    const mixed scale = 0.5f*dtVel/(mixed) 0x100000000;
-    int index = GLOBAL_ID;
-    if (index == 0)
-        dt[0].x = stepSize.y;
-    while(index < numAtoms) {
-        int atom = atomList[index];
-        mixed4 velocity = velm[atom];
-        if (velocity.w != 0.0) {
-            mixed4 deltaXconstrained = posDelta[atom];
-            velocity.x += scale*force[atom]*velocity.w + (deltaXconstrained.x - velocity.x*stepSize.y)*oneOverDt;
-            velocity.y += scale*force[atom+paddedNumAtoms]*velocity.w + (deltaXconstrained.y - velocity.y*stepSize.y)*oneOverDt;
-            velocity.z += scale*force[atom+paddedNumAtoms*2]*velocity.w + (deltaXconstrained.z - velocity.z*stepSize.y)*oneOverDt;
-#ifndef SUPPORTS_DOUBLE_PRECISION
-            velocity.x += (deltaXconstrained.x - velocity.x*stepSize.y)*correction;
-            velocity.y += (deltaXconstrained.y - velocity.y*stepSize.y)*correction;
-            velocity.z += (deltaXconstrained.z - velocity.z*stepSize.y)*correction;
-#endif
-            velm[atom] = velocity;
-        }
-        index += GLOBAL_SIZE;
-    }
-    index = GLOBAL_ID;
-    while(index < numPairs) {
-        int atom1 = pairList[index].x;
-        int atom2 = pairList[index].y;
-        mixed4 v1 = velm[atom1];
-        mixed4 v2 = velm[atom2];
-        mixed m1 = v1.w == 0.0f ? 0.0f : 1.0f / v1.w;
-        mixed m2 = v2.w == 0.0f ? 0.0f : 1.0f / v2.w;
-        mixed mass1fract = m1 / (m1 + m2);
-        mixed mass2fract = m2 / (m1 + m2);
-        mixed invRedMass = (m1 * m2 != 0.0f) ? (m1 + m2)/(m1 * m2) : 0.0f;
-        mixed invTotMass = (m1 + m2 != 0.0f) ? 1.0f /(m1 + m2) : 0.0f;
-        mixed3 comVel;
-        comVel.x= v1.x*mass1fract + v2.x*mass2fract;
-        comVel.y= v1.y*mass1fract + v2.y*mass2fract;
-        comVel.z= v1.z*mass1fract + v2.z*mass2fract;
-        mixed3 relVel;
-        relVel.x= v2.x - v1.x;
-        relVel.y= v2.y - v1.y;
-        relVel.z= v2.z - v1.z;
-
-        mixed3 comFrc;
-        mixed F1x = scale*force[atom1];
-        mixed F1y = scale*force[atom1+paddedNumAtoms];
-        mixed F1z = scale*force[atom1+paddedNumAtoms*2];
-        mixed F2x = scale*force[atom2];
-        mixed F2y = scale*force[atom2+paddedNumAtoms];
-        mixed F2z = scale*force[atom2+paddedNumAtoms*2];
-        comFrc.x = F1x + F2x;
-        comFrc.y = F1y + F2y;
-        comFrc.z = F1z + F2z;
-        mixed3 relFrc;
-        relFrc.x = mass1fract*F2x - mass2fract*F1x;
-        relFrc.y = mass1fract*F2y - mass2fract*F1y;
-        relFrc.z = mass1fract*F2z - mass2fract*F1z;
-        comVel.x += comFrc.x * invTotMass;
-        comVel.y += comFrc.y * invTotMass;
-        comVel.z += comFrc.z * invTotMass;
-        relVel.x += relFrc.x * invRedMass;
-        relVel.y += relFrc.y * invRedMass;
-        relVel.z += relFrc.z * invRedMass;
-        if (v1.w != 0.0f) {
-            mixed4 deltaXconstrained = posDelta[atom1];
-            v1.x = comVel.x - relVel.x*mass2fract + (deltaXconstrained.x - v1.x*stepSize.y)*oneOverDt;
-            v1.y = comVel.y - relVel.y*mass2fract + (deltaXconstrained.y - v1.y*stepSize.y)*oneOverDt;
-            v1.z = comVel.z - relVel.z*mass2fract + (deltaXconstrained.z - v1.z*stepSize.y)*oneOverDt;
-#ifndef SUPPORTS_DOUBLE_PRECISION
-            v1.x += (deltaXconstrained.x - v1.x*stepSize.y)*correction;
-            v1.y += (deltaXconstrained.y - v1.y*stepSize.y)*correction;
-            v1.z += (deltaXconstrained.z - v1.z*stepSize.y)*correction;
-#endif
-            velm[atom1] = v1;
-        }
-        if (v2.w != 0.0f) {
-            mixed4 deltaXconstrained = posDelta[atom2];
-            v2.x = comVel.x + relVel.x*mass1fract + (deltaXconstrained.x - v2.x*stepSize.y)*oneOverDt;
-            v2.y = comVel.y + relVel.y*mass1fract + (deltaXconstrained.y - v2.y*stepSize.y)*oneOverDt;
-            v2.z = comVel.z + relVel.z*mass1fract + (deltaXconstrained.z - v2.z*stepSize.y)*oneOverDt;
-#ifndef SUPPORTS_DOUBLE_PRECISION
-            v2.x += (deltaXconstrained.x - v2.x*stepSize.y)*correction;
-            v2.y += (deltaXconstrained.y - v2.y*stepSize.y)*correction;
-            v2.z += (deltaXconstrained.z - v2.z*stepSize.y)*correction;
-#endif
-            velm[atom2] = v2;
-        }
-        index += GLOBAL_SIZE;
-    }
-}
-
-KERNEL void integrateVelocityVerletHardWall(int numPairs, GLOBAL const float* RESTRICT maxPairDistance, 
+KERNEL void integrateNoseHooverHardWall(int numPairs, GLOBAL const float* RESTRICT maxPairDistance, 
                                         GLOBAL mixed2* RESTRICT dt, GLOBAL real4* RESTRICT posq,
                                         GLOBAL mixed4* RESTRICT velm, GLOBAL const int2* RESTRICT pairList,
                                         GLOBAL const float* RESTRICT pairTemperature
@@ -370,4 +258,3 @@ KERNEL void integrateVelocityVerletHardWall(int numPairs, GLOBAL const float* RE
         }
     }
 }
-

platforms/cuda/src/CudaKernelFactory.cpp

@@ -133,10 +133,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CommonIntegrateCustomStepKernel(name, platform, cu);
     if (name == ApplyAndersenThermostatKernel::Name())
         return new CommonApplyAndersenThermostatKernel(name, platform, cu);
-    if (name == NoseHooverChainKernel::Name())
-        return new CommonNoseHooverChainKernel(name, platform, cu);
-    if (name == IntegrateVelocityVerletStepKernel::Name())
-        return new CommonIntegrateVelocityVerletStepKernel(name, platform, cu);
+    if (name == IntegrateNoseHooverStepKernel::Name())
+        return new CommonIntegrateNoseHooverStepKernel(name, platform, cu);
     if (name == ApplyMonteCarloBarostatKernel::Name())
         return new CudaApplyMonteCarloBarostatKernel(name, platform, cu);
     if (name == RemoveCMMotionKernel::Name())

platforms/cuda/src/CudaPlatform.cpp

@@ -96,7 +96,7 @@ CudaPlatform::CudaPlatform() {
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
-    registerKernelFactory(IntegrateVelocityVerletStepKernel::Name(), factory);
+    registerKernelFactory(IntegrateNoseHooverStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinMiddleStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
@@ -104,7 +104,6 @@ CudaPlatform::CudaPlatform() {
     registerKernelFactory(IntegrateVariableLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateCustomStepKernel::Name(), factory);
     registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
-    registerKernelFactory(NoseHooverChainKernel::Name(), factory);
     registerKernelFactory(ApplyMonteCarloBarostatKernel::Name(), factory);
     registerKernelFactory(RemoveCMMotionKernel::Name(), factory);
     platformProperties.push_back(CudaDeviceIndex());

platforms/cuda/tests/TestCudaNoseHooverThermostat.cpp

-/* -------------------------------------------------------------------------- *
- *                                   OpenMM                                   *
- * -------------------------------------------------------------------------- *
- * This is part of the OpenMM molecular simulation toolkit originating from   *
- * Simbios, the NIH National Center for Physics-Based Simulation of           *
- * Biological Structures at Stanford, funded under the NIH Roadmap for        *
- * Medical Research, grant U54 GM072970. See https://simtk.org.               *
- *                                                                            *
- * Portions copyright (c) 2019 Stanford University and the Authors.           *
- * Authors: Andreas Krämer and Andrew C. Simmonett                            *
- * Contributors:                                                              *
- *                                                                            *
- * 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 "CudaTests.h"
-#include "TestNoseHooverThermostat.h"
-
-void runPlatformTests() {
-}

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -131,10 +131,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
         return new CommonIntegrateCustomStepKernel(name, platform, cl);
     if (name == ApplyAndersenThermostatKernel::Name())
         return new CommonApplyAndersenThermostatKernel(name, platform, cl);
-    if (name == NoseHooverChainKernel::Name())
-        return new CommonNoseHooverChainKernel(name, platform, cl);
-    if (name == IntegrateVelocityVerletStepKernel::Name())
-        return new CommonIntegrateVelocityVerletStepKernel(name, platform, cl);
+    if (name == IntegrateNoseHooverStepKernel::Name())
+        return new CommonIntegrateNoseHooverStepKernel(name, platform, cl);
     if (name == ApplyMonteCarloBarostatKernel::Name())
         return new OpenCLApplyMonteCarloBarostatKernel(name, platform, cl);
     if (name == RemoveCMMotionKernel::Name())

platforms/opencl/src/OpenCLPlatform.cpp

@@ -87,7 +87,7 @@ OpenCLPlatform::OpenCLPlatform() {
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
-    registerKernelFactory(IntegrateVelocityVerletStepKernel::Name(), factory);
+    registerKernelFactory(IntegrateNoseHooverStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinMiddleStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
@@ -95,7 +95,6 @@ OpenCLPlatform::OpenCLPlatform() {
     registerKernelFactory(IntegrateVariableLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateCustomStepKernel::Name(), factory);
     registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
-    registerKernelFactory(NoseHooverChainKernel::Name(), factory);
     registerKernelFactory(ApplyMonteCarloBarostatKernel::Name(), factory);
     registerKernelFactory(RemoveCMMotionKernel::Name(), factory);
     platformProperties.push_back(OpenCLDeviceIndex());

platforms/opencl/tests/TestOpenCLNoseHooverThermostat.cpp

-/* -------------------------------------------------------------------------- *
- *                                   OpenMM                                   *
- * -------------------------------------------------------------------------- *
- * This is part of the OpenMM molecular simulation toolkit originating from   *
- * Simbios, the NIH National Center for Physics-Based Simulation of           *
- * Biological Structures at Stanford, funded under the NIH Roadmap for        *
- * Medical Research, grant U54 GM072970. See https://simtk.org.               *
- *                                                                            *
- * Portions copyright (c) 2019 Stanford University and the Authors.           *
- * Authors: Andreas Krämer and Andrew C. Simmonett                            *
- * Contributors:                                                              *
- *                                                                            *
- * 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 "OpenCLTests.h"
-#include "TestNoseHooverThermostat.h"
-
-void runPlatformTests() {
-}

platforms/reference/include/ReferenceKernels.h

@@ -61,7 +61,7 @@ class ReferenceStochasticDynamics;
 class ReferenceConstraintAlgorithm;
 class ReferenceNoseHooverChain;
 class ReferenceMonteCarloBarostat;
-class ReferenceVelocityVerletDynamics;
+class ReferenceNoseHooverDynamics;
 class ReferenceVariableStochasticDynamics;
 class ReferenceVariableVerletDynamics;
 class ReferenceVerletDynamics;
@@ -1139,12 +1139,12 @@ private:
 /**
  * This kernel is invoked by NoseHooverIntegrator to take one time step.
  */
-class ReferenceIntegrateVelocityVerletStepKernel : public IntegrateVelocityVerletStepKernel {
+class ReferenceIntegrateNoseHooverStepKernel : public IntegrateNoseHooverStepKernel {
 public:
-    ReferenceIntegrateVelocityVerletStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVelocityVerletStepKernel(name, platform),
+    ReferenceIntegrateNoseHooverStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateNoseHooverStepKernel(name, platform),
         data(data), dynamics(0) {
     }
-    ~ReferenceIntegrateVelocityVerletStepKernel();
+    ~ReferenceIntegrateNoseHooverStepKernel();
     /**
      * Initialize the kernel.
      * 
@@ -1168,9 +1168,45 @@ public:
      * @param integrator the NoseHooverIntegrator this kernel is being used for
      */
     double computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator);
+    /**
+     * Execute the kernel that propagates the Nose Hoover chain and determines the velocity scale factor.
+     * 
+     * @param context  the context in which to execute this kernel
+     * @param noseHooverChain the object describing the chain to be propagated.
+     * @param kineticEnergy the {center of mass, relative} kineticEnergies of the particles being thermostated by this chain.
+     * @param timeStep the time step used by the integrator.
+     * @return the velocity scale factor to apply to the particles associated with this heat bath.
+     */
+    std::pair<double, double> propagateChain(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> kineticEnergy, double timeStep);
+    /**
+     * Execute the kernal that computes the total (kinetic + potential) heat bath energy.
+     *
+     * @param context the context in which to execute this kernel
+     * @param noseHooverChain the chain whose energy is to be determined.
+     * @return the total heat bath energy.
+     */
+    double computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain);
+    /**
+     * Execute the kernel that computes the kinetic energy for a subset of atoms,
+     * or the relative kinetic energy of Drude particles with respect to their parent atoms
+     *
+     * @param context the context in which to execute this kernel
+     * @param noseHooverChain the chain whose energy is to be determined.
+     * @param downloadValue whether the computed value should be downloaded and returned.
+     */
+    std::pair<double, double> computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue);
+    /**
+     * Execute the kernel that scales the velocities of particles associated with a nose hoover chain
+     *
+     * @param context the context in which to execute this kernel
+     * @param noseHooverChain the chain whose energy is to be determined.
+     * @param scaleFactor the multiplicative factor by which {absolute, relative} velocities are scaled.
+     */
+    void scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactor);
 private:
     ReferencePlatform::PlatformData& data;
-    ReferenceVelocityVerletDynamics* dynamics;
+    ReferenceNoseHooverChain* chainPropagator;
+    ReferenceNoseHooverDynamics* dynamics;
     std::vector<double> masses;
     double prevStepSize;
 };
@@ -1466,62 +1502,6 @@ private:
     std::vector<double> masses;
 };
 
-/**
- * This kernel is invoked by NoseHooverChain at the start of each time step to adjust the thermostat
- * and update the associated particle velocities.
- */
-class ReferenceNoseHooverChainKernel : public NoseHooverChainKernel {
-public:
-    ReferenceNoseHooverChainKernel(std::string name, const Platform& platform) : NoseHooverChainKernel(name, platform), chainPropagator(0) {
-    }
-    ~ReferenceNoseHooverChainKernel();
-    /**
-     * Initialize the kernel.
-     */
-    virtual void initialize();
-    /**
-     * Execute the kernel that propagates the Nose Hoover chain and determines the velocity scale factor.
-     * 
-     * @param context  the context in which to execute this kernel
-     * @param noseHooverChain the object describing the chain to be propagated.
-     * @param kineticEnergy the {absolute, relative} kinetic energies of the particles being thermostated by this chain.
-     * @param timeStep the time step used by the integrator.
-     * @return the velocity scale factors to apply to the {absolute, relative} motion of particles associated with this heat bath.
-     */
-    virtual std::pair<double, double> propagateChain(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> kineticEnergy, double timeStep);
-    /**
-     * Execute the kernal that computes the total (kinetic + potential) heat bath energy.
-     *
-     * @param context the context in which to execute this kernel
-     * @param noseHooverChain the chain whose energy is to be determined.
-     * @return the total heat bath energy.
-     */
-    virtual double computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &nhc);
-    /**
-     * Execute the kernel that computes the kinetic energy for a subset of atoms,
-     * or the relative kinetic energy of Drude particles with respect to their parent atoms
-     *
-     * @param context the context in which to execute this kernel
-     * @param noseHooverChain the chain whose energy is to be determined.
-     * @param downloadValue whether the computed value should be downloaded and returned.
-     *
-     */
-     virtual std::pair<double, double> computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue);
-
-    /**
-     * Execute the kernel that scales the velocities of particles associated with a nose hoover chain
-     *
-     * @param context the context in which to execute this kernel
-     * @param noseHooverChain the chain whose energy is to be determined.
-     * @param scaleFactor the multiplicative factor by which {absolute, relative} velocities are scaled.
-     */
-    virtual void scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactor);
-
-
-private:
-    ReferenceNoseHooverChain* chainPropagator;
-};
-
 /**
  * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
  */

platforms/reference/include/ReferenceVelocityVerletDynamics.h → platforms/reference/include/ReferenceNoseHooverDynamics.h

@@ -22,20 +22,23 @@
  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  */
 
-#ifndef __ReferenceVelocityVerletDynamics_H__
-#define __ReferenceVelocityVerletDynamics_H__
+#ifndef __ReferenceNoseHooverDynamics_H__
+#define __ReferenceNoseHooverDynamics_H__
 
 #include "ReferenceDynamics.h"
+#include <tuple>
 
 namespace OpenMM {
 
 class ContextImpl;
 
-class ReferenceVelocityVerletDynamics : public ReferenceDynamics {
+class ReferenceNoseHooverDynamics : public ReferenceDynamics {
 
    private:
       std::vector<OpenMM::Vec3> xPrime;
+      std::vector<OpenMM::Vec3> oldx;
       std::vector<double> inverseMasses;
+      int numberOfAtoms;
 
    public:
 
@@ -50,7 +53,7 @@ class ReferenceVelocityVerletDynamics : public ReferenceDynamics {
       
          --------------------------------------------------------------------------------------- */
 
-       ReferenceVelocityVerletDynamics(int numberOfAtoms, double deltaT);
+       ReferenceNoseHooverDynamics(int numberOfAtoms, double deltaT);
 
       /**---------------------------------------------------------------------------------------
       
@@ -58,11 +61,11 @@ class ReferenceVelocityVerletDynamics : public ReferenceDynamics {
       
          --------------------------------------------------------------------------------------- */
 
-       ~ReferenceVelocityVerletDynamics();
+       ~ReferenceNoseHooverDynamics();
 
       /**---------------------------------------------------------------------------------------
       
-         Update
+         BAstep = update the velocities and positions to start a BAOAB step using the leapfrog "middle" scheme
       
          @param system              the System to be integrated
          @param atomCoordinates     atom coordinates
@@ -77,7 +80,26 @@ class ReferenceVelocityVerletDynamics : public ReferenceDynamics {
       
          --------------------------------------------------------------------------------------- */
      
-      void update(OpenMM::ContextImpl &context, const OpenMM::System& system, std::vector<OpenMM::Vec3>& atomCoordinates,
+      void BAstep(OpenMM::ContextImpl &context, const OpenMM::System& system, std::vector<OpenMM::Vec3>& atomCoordinates,
+                  std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, double tolerance, bool &forcesAreValid,
+                  const std::vector<int> & allAtoms, const std::vector<std::tuple<int, int, double>> & allPairs, double maxPairDistance);
+      /**---------------------------------------------------------------------------------------
+      
+         ABstep = update the positions and velocities to end a BAOAB step using the leapfrog "middle" scheme
+      
+         @param system              the System to be integrated
+         @param atomCoordinates     atom coordinates
+         @param velocities          velocities
+         @param forces              forces
+         @param masses              atom masses
+         @param tolerance           the constraint tolerance
+         @param forcesAreValid      whether the forces are valid (duh!)
+         @param allAtoms            a list of all atoms not involved in a Drude-like pair
+         @param allPairs            a list of all Drude-like pairs, and their KT values, in the system
+         @param maxPairDistance     the maximum separation allowed for a Drude-like pair
+      
+         --------------------------------------------------------------------------------------- */
+      void ABstep(OpenMM::ContextImpl &context, const OpenMM::System& system, std::vector<OpenMM::Vec3>& atomCoordinates,
                   std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, double tolerance, bool &forcesAreValid,
                   const std::vector<int> & allAtoms, const std::vector<std::tuple<int, int, double>> & allPairs, double maxPairDistance);
       
@@ -85,4 +107,4 @@ class ReferenceVelocityVerletDynamics : public ReferenceDynamics {
 
 } // namespace OpenMM
 
-#endif // __ReferenceVelocityVerletDynamics_H__
+#endif // __ReferenceNoseHooverDynamics_H__

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -88,10 +88,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcGayBerneForceKernel(name, platform);
     if (name == IntegrateVerletStepKernel::Name())
         return new ReferenceIntegrateVerletStepKernel(name, platform, data);
-    if (name == IntegrateVelocityVerletStepKernel::Name())
-        return new ReferenceIntegrateVelocityVerletStepKernel(name, platform, data);
-    if (name == NoseHooverChainKernel::Name())
-        return new ReferenceNoseHooverChainKernel(name, platform);
+    if (name == IntegrateNoseHooverStepKernel::Name())
+        return new ReferenceIntegrateNoseHooverStepKernel(name, platform, data);
     if (name == IntegrateLangevinStepKernel::Name())
         return new ReferenceIntegrateLangevinStepKernel(name, platform, data);
     if (name == IntegrateLangevinMiddleStepKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -57,6 +57,7 @@
 #include "ReferenceLJCoulombIxn.h"
 #include "ReferenceMonteCarloBarostat.h"
 #include "ReferenceNoseHooverChain.h"
+#include "ReferenceNoseHooverDynamics.h"
 #include "ReferenceProperDihedralBond.h"
 #include "ReferenceRbDihedralBond.h"
 #include "ReferenceRMSDForce.h"
@@ -64,7 +65,6 @@
 #include "ReferenceTabulatedFunction.h"
 #include "ReferenceVariableStochasticDynamics.h"
 #include "ReferenceVariableVerletDynamics.h"
-#include "ReferenceVelocityVerletDynamics.h"
 #include "ReferenceVerletDynamics.h"
 #include "ReferenceVirtualSites.h"
 #include "openmm/CMMotionRemover.h"
@@ -2153,19 +2153,22 @@ double ReferenceIntegrateVerletStepKernel::computeKineticEnergy(ContextImpl& con
     return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
 }
 
-ReferenceIntegrateVelocityVerletStepKernel::~ReferenceIntegrateVelocityVerletStepKernel() {
+ReferenceIntegrateNoseHooverStepKernel::~ReferenceIntegrateNoseHooverStepKernel() {
+    if (chainPropagator)
+        delete chainPropagator;
     if (dynamics)
         delete dynamics;
 }
 
-void ReferenceIntegrateVelocityVerletStepKernel::initialize(const System& system, const NoseHooverIntegrator& integrator) {
+void ReferenceIntegrateNoseHooverStepKernel::initialize(const System& system, const NoseHooverIntegrator& integrator) {
     int numParticles = system.getNumParticles();
     masses.resize(numParticles);
     for (int i = 0; i < numParticles; ++i)
         masses[i] = system.getParticleMass(i);
+    this->chainPropagator = new ReferenceNoseHooverChain();
 }
 
-void ReferenceIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) {
+void ReferenceIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) {
     double stepSize = integrator.getStepSize();
     vector<Vec3>& posData = extractPositions(context);
     vector<Vec3>& velData = extractVelocities(context);
@@ -2174,20 +2177,214 @@ void ReferenceIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, c
         // Recreate the computation objects with the new parameters.
         if (dynamics)
             delete dynamics;
-        dynamics = new ReferenceVelocityVerletDynamics(context.getSystem().getNumParticles(), stepSize);
+        dynamics = new ReferenceNoseHooverDynamics(context.getSystem().getNumParticles(), stepSize);
         dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
         prevStepSize = stepSize;
     }
-    dynamics->update(context, context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance(), forcesAreValid,
+
+    // BA
+    dynamics->BAstep(context, context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance(), forcesAreValid,
+                     integrator.getAllThermostatedIndividualParticles(), integrator.getAllThermostatedPairs(), integrator.getMaximumPairDistance());
+    // O
+    int numChains = integrator.getNumThermostats();
+    for(int chain = 0; chain < numChains; ++chain) {
+        const auto &thermostatChain = integrator.getThermostat(chain);
+        std::pair<double, double> KEs = computeMaskedKineticEnergy(context, thermostatChain, true);
+        std::pair<double, double> scaleFactors = propagateChain(context, thermostatChain, KEs, stepSize);
+        scaleVelocities(context, thermostatChain, scaleFactors);
+    }
+    // AB
+    dynamics->ABstep(context, context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance(), forcesAreValid,
                      integrator.getAllThermostatedIndividualParticles(), integrator.getAllThermostatedPairs(), integrator.getMaximumPairDistance());
+
     data.time += stepSize;
     data.stepCount++;
 }
 
-double ReferenceIntegrateVelocityVerletStepKernel::computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator) {
+double ReferenceIntegrateNoseHooverStepKernel::computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator) {
     return computeShiftedKineticEnergy(context, masses, 0);
 }
 
+std::pair<double, double> ReferenceIntegrateNoseHooverStepKernel::propagateChain(ContextImpl& context, const NoseHooverChain &nhc,
+                                                                                     std::pair<double, double> kineticEnergy, double timeStep) {
+    double absKE = kineticEnergy.first;
+    double relKE = kineticEnergy.second;
+    if (absKE < 1e-8) return {1.0, 1.0};  // (catches the problem of zero velocities in the first dynamics step, where we have nothing to scale)
+    // Get the variables describing the NHC
+    int chainLength = nhc.getChainLength();
+    int chainID = nhc.getChainID();
+    int numDOFs = nhc.getNumDegreesOfFreedom();
+    int numMTS = nhc.getNumMultiTimeSteps();
+
+    // Get the state of the NHC from the context
+    auto& allChainPositions = extractNoseHooverPositions(context);
+    auto& allChainVelocities = extractNoseHooverVelocities(context);
+
+    int nAtoms = nhc.getThermostatedAtoms().size();
+    double absScale = 0;
+    if (nAtoms) {
+        if (allChainPositions.size() < 2*chainID+1){
+            allChainPositions.resize(2*chainID+1);
+        }
+        if (allChainVelocities.size() < 2*chainID+1){
+            allChainVelocities.resize(2*chainID+1);
+        }
+        auto& chainPositions = allChainPositions.at(2*chainID);
+        auto& chainVelocities = allChainVelocities.at(2*chainID);
+        if (chainPositions.size() < chainLength){
+            chainPositions.resize(chainLength, 0);
+        }
+        if (chainVelocities.size() < chainLength){
+            chainVelocities.resize(chainLength, 0);
+        }
+        double temperature = nhc.getTemperature();
+        double collisionFrequency = nhc.getCollisionFrequency();
+        absScale = chainPropagator->propagate(absKE, chainVelocities, chainPositions, numDOFs,
+                                              temperature, collisionFrequency, timeStep,
+                                              numMTS, nhc.getYoshidaSuzukiWeights());
+    }
+    double relScale = 0;
+    int nPairs = nhc.getThermostatedPairs().size();
+    if (nPairs) {
+        if (allChainPositions.size() < 2*chainID+2){
+            allChainPositions.resize(2*chainID+2);
+        }
+        if (allChainVelocities.size() < 2*chainID+2){
+            allChainVelocities.resize(2*chainID+2);
+        }
+        auto& chainPositions = allChainPositions.at(2*chainID+1);
+        auto& chainVelocities = allChainVelocities.at(2*chainID+1);
+        if (chainPositions.size() < chainLength){
+            chainPositions.resize(chainLength, 0);
+        }
+        if (chainVelocities.size() < chainLength){
+            chainVelocities.resize(chainLength, 0);
+        }
+        double temperature = nhc.getRelativeTemperature();
+        double collisionFrequency = nhc.getRelativeCollisionFrequency();
+        relScale = chainPropagator->propagate(relKE, chainVelocities, chainPositions, 3*nPairs,
+                                              temperature, collisionFrequency, timeStep,
+                                              numMTS, nhc.getYoshidaSuzukiWeights());
+    }
+    return {absScale, relScale};
+}
+
+double ReferenceIntegrateNoseHooverStepKernel::computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &nhc) {
+    double potentialEnergy = 0;
+    double kineticEnergy = 0;
+    int chainLength = nhc.getChainLength();
+    int chainID = nhc.getChainID();
+    int nAtoms = nhc.getThermostatedAtoms().size();
+    int nPairs = nhc.getThermostatedPairs().size();
+    auto& nhcPositions = extractNoseHooverPositions(context);
+    auto& nhcVelocities = extractNoseHooverVelocities(context);
+    if (nAtoms) {
+        double temperature = nhc.getTemperature();
+        double collisionFrequency = nhc.getCollisionFrequency();
+        double kT = temperature * BOLTZ;
+        int numDOFs = nhc.getNumDegreesOfFreedom();
+        for(int i = 0; i < chainLength; ++i) {
+            double prefac = i ? 1 : numDOFs;
+            double mass = prefac * kT / (collisionFrequency * collisionFrequency);
+            double velocity = nhcVelocities[2*chainID][i];
+            // The kinetic energy of this bead
+            kineticEnergy += 0.5 * mass * velocity * velocity;
+            // The potential energy of this bead
+            double position = nhcPositions[2*chainID][i];
+            potentialEnergy += prefac * kT * position;
+        }
+    }
+    if (nPairs) {
+        double temperature = nhc.getRelativeTemperature();
+        double collisionFrequency = nhc.getRelativeCollisionFrequency();
+        double kT = temperature * BOLTZ;
+        int numDOFs = 3 * nPairs;
+        for(int i = 0; i < chainLength; ++i) {
+            double prefac = i ? 1 : numDOFs;
+            double mass = prefac * kT / (collisionFrequency * collisionFrequency);
+            double velocity = nhcVelocities[2*chainID+1][i];
+            // The kinetic energy of this bead
+            kineticEnergy += 0.5 * mass * velocity * velocity;
+            // The potential energy of this bead
+            double position = nhcPositions[2*chainID+1][i];
+            potentialEnergy += prefac * kT * position;
+        }
+    }
+    return kineticEnergy + potentialEnergy;
+}
+
+std::pair<double, double> ReferenceIntegrateNoseHooverStepKernel::computeMaskedKineticEnergy(ContextImpl& context,
+                                                                const NoseHooverChain &noseHooverChain, bool downloadValue) {
+    const std::vector<int>& atomsList = noseHooverChain.getThermostatedAtoms();
+    const std::vector<std::pair<int,int>>& pairsList = noseHooverChain.getThermostatedPairs();
+    std::vector<Vec3>& velocities = extractVelocities(context);
+    const System& system = context.getSystem();
+    int numParticles = system.getNumParticles();
+    std::vector<double> masses(numParticles);
+    for (int i = 0; i < numParticles; ++i)
+        masses[i] = system.getParticleMass(i);
+
+    double comKE = 0;
+    double relKE = 0;
+    // kinetic energy of individual atoms
+    for (const auto &m: atomsList){
+        comKE += 0.5 * masses[m] * velocities[m].dot(velocities[m]);
+    }
+    // center of mass kinetic energy of pairs
+    for (const auto &p: pairsList){
+        double m1 = masses[p.first];
+        double m2 = masses[p.second];
+        Vec3 v1 = velocities[p.first];
+        Vec3 v2 = velocities[p.second];
+        double invMass = 1.0 / (m1 + m2);
+        double redMass = m1 * m2 * invMass;
+        double fracM1 = m1 * invMass;
+        double fracM2 = m2 * invMass;
+        Vec3 comVelocity = fracM1 * v1 + fracM2 * v2;
+        Vec3 relVelocity = v2 - v1;
+
+        comKE += 0.5 * (m1 + m2) * comVelocity.dot(comVelocity);
+        relKE += 0.5 * redMass * relVelocity.dot(relVelocity);
+    }
+    // We ignore the downloadValue argument here and always return the correct value
+    return {comKE, relKE};
+}
+
+
+void ReferenceIntegrateNoseHooverStepKernel::scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactors) {
+    const auto& atoms = noseHooverChain.getThermostatedAtoms();
+    const auto& pairs = noseHooverChain.getThermostatedPairs();
+    std::vector<Vec3>& velocities = extractVelocities(context);
+    double absScale = scaleFactors.first;
+    double relScale = scaleFactors.second;
+
+    const System& system = context.getSystem();
+    int numParticles = system.getNumParticles();
+    std::vector<double> masses(numParticles);
+    for (int i = 0; i < numParticles; ++i)
+        masses[i] = system.getParticleMass(i);
+    // scale absolute velocities
+    for (const auto &a: atoms){
+        velocities[a] *= absScale;
+    }
+    // scale relative velocities and absolute center of mass velocities for each pair
+    for (const auto &p: pairs){
+        int p1 = p.first;
+        int p2 = p.second;
+        double m1 = masses[p.first];
+        double m2 = masses[p.second];
+        Vec3 v1 = velocities[p.first];
+        Vec3 v2 = velocities[p.second];
+        double invMass = 1.0 / (m1 + m2);
+        double fracM1 = m1 * invMass;
+        double fracM2 = m2 * invMass;
+        Vec3 comVelocity = fracM1 * v1 + fracM2 * v2;
+        Vec3 relVelocity = v2 - v1;
+        velocities[p1] = absScale * comVelocity - relScale * relVelocity * fracM2;
+        velocities[p2] = absScale * comVelocity + relScale * relVelocity * fracM1;
+    }
+}
+
 ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
     if (dynamics)
         delete dynamics;
@@ -2512,195 +2709,6 @@ void ReferenceApplyAndersenThermostatKernel::execute(ContextImpl& context) {
         context.getIntegrator().getStepSize());
 }
 
-ReferenceNoseHooverChainKernel::~ReferenceNoseHooverChainKernel() {
-    if (chainPropagator)
-        delete chainPropagator;
-}
-
-void ReferenceNoseHooverChainKernel::initialize() {
-    this->chainPropagator = new ReferenceNoseHooverChain();
-    //SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
-}
-
-std::pair<double, double> ReferenceNoseHooverChainKernel::propagateChain(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> kineticEnergy, double timeStep) {
-    double absKE = kineticEnergy.first;
-    double relKE = kineticEnergy.second;
-    if (absKE < 1e-8) return {1.0, 1.0};  // (catches the problem of zero velocities in the first dynamics step, where we have nothing to scale)
-    // Get the variables describing the NHC
-    int chainLength = nhc.getChainLength();
-    int chainID = nhc.getChainID();
-    int numDOFs = nhc.getNumDegreesOfFreedom();
-    int numMTS = nhc.getNumMultiTimeSteps();
-
-    // Get the state of the NHC from the context
-    auto& allChainPositions = extractNoseHooverPositions(context);
-    auto& allChainVelocities = extractNoseHooverVelocities(context);
-
-    int nAtoms = nhc.getThermostatedAtoms().size();
-    double absScale = 0;
-    if (nAtoms) {
-        if (allChainPositions.size() < 2*chainID+1){
-            allChainPositions.resize(2*chainID+1);
-        }
-        if (allChainVelocities.size() < 2*chainID+1){
-            allChainVelocities.resize(2*chainID+1);
-        }
-        auto& chainPositions = allChainPositions.at(2*chainID);
-        auto& chainVelocities = allChainVelocities.at(2*chainID);
-        if (chainPositions.size() < chainLength){
-            chainPositions.resize(chainLength, 0);
-        }
-        if (chainVelocities.size() < chainLength){
-            chainVelocities.resize(chainLength, 0);
-        }
-        double temperature = nhc.getTemperature();
-        double collisionFrequency = nhc.getCollisionFrequency();
-        absScale = chainPropagator->propagate(absKE, chainVelocities, chainPositions, numDOFs,
-                                              temperature, collisionFrequency, timeStep,
-                                              numMTS, nhc.getYoshidaSuzukiWeights());
-    }
-    double relScale = 0;
-    int nPairs = nhc.getThermostatedPairs().size();
-    if (nPairs) {
-        if (allChainPositions.size() < 2*chainID+2){
-            allChainPositions.resize(2*chainID+2);
-        }
-        if (allChainVelocities.size() < 2*chainID+2){
-            allChainVelocities.resize(2*chainID+2);
-        }
-        auto& chainPositions = allChainPositions.at(2*chainID+1);
-        auto& chainVelocities = allChainVelocities.at(2*chainID+1);
-        if (chainPositions.size() < chainLength){
-            chainPositions.resize(chainLength, 0);
-        }
-        if (chainVelocities.size() < chainLength){
-            chainVelocities.resize(chainLength, 0);
-        }
-        double temperature = nhc.getRelativeTemperature();
-        double collisionFrequency = nhc.getRelativeCollisionFrequency();
-        relScale = chainPropagator->propagate(relKE, chainVelocities, chainPositions, 3*nPairs,
-                                              temperature, collisionFrequency, timeStep,
-                                              numMTS, nhc.getYoshidaSuzukiWeights());
-    }
-    return {absScale, relScale};
-}
-
-double ReferenceNoseHooverChainKernel::computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &nhc) {
-    double potentialEnergy = 0;
-    double kineticEnergy = 0;
-    int chainLength = nhc.getChainLength();
-    int chainID = nhc.getChainID();
-    int nAtoms = nhc.getThermostatedAtoms().size();
-    int nPairs = nhc.getThermostatedPairs().size();
-    auto& nhcPositions = extractNoseHooverPositions(context);
-    auto& nhcVelocities = extractNoseHooverVelocities(context);
-    if (nAtoms) {
-        double temperature = nhc.getTemperature();
-        double collisionFrequency = nhc.getCollisionFrequency();
-        double kT = temperature * BOLTZ;
-        int numDOFs = nhc.getNumDegreesOfFreedom();
-        for(int i = 0; i < chainLength; ++i) {
-            double prefac = i ? 1 : numDOFs;
-            double mass = prefac * kT / (collisionFrequency * collisionFrequency);
-            double velocity = nhcVelocities[2*chainID][i];
-            // The kinetic energy of this bead
-            kineticEnergy += 0.5 * mass * velocity * velocity;
-            // The potential energy of this bead
-            double position = nhcPositions[2*chainID][i];
-            potentialEnergy += prefac * kT * position;
-        }
-    }
-    if (nPairs) {
-        double temperature = nhc.getRelativeTemperature();
-        double collisionFrequency = nhc.getRelativeCollisionFrequency();
-        double kT = temperature * BOLTZ;
-        int numDOFs = 3 * nPairs;
-        for(int i = 0; i < chainLength; ++i) {
-            double prefac = i ? 1 : numDOFs;
-            double mass = prefac * kT / (collisionFrequency * collisionFrequency);
-            double velocity = nhcVelocities[2*chainID+1][i];
-            // The kinetic energy of this bead
-            kineticEnergy += 0.5 * mass * velocity * velocity;
-            // The potential energy of this bead
-            double position = nhcPositions[2*chainID+1][i];
-            potentialEnergy += prefac * kT * position;
-        }
-    }
-    return kineticEnergy + potentialEnergy;
-}
-
-std::pair<double, double> ReferenceNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue) {
-    const std::vector<int>& atomsList = noseHooverChain.getThermostatedAtoms();
-    const std::vector<std::pair<int,int>>& pairsList = noseHooverChain.getThermostatedPairs();
-    std::vector<Vec3>& velocities = extractVelocities(context);
-    const System& system = context.getSystem();
-    int numParticles = system.getNumParticles();
-    std::vector<double> masses(numParticles);
-    for (int i = 0; i < numParticles; ++i)
-        masses[i] = system.getParticleMass(i);
-
-    double comKE = 0;
-    double relKE = 0;
-    // kinetic energy of individual atoms
-    for (const auto &m: atomsList){
-        comKE += 0.5 * masses[m] * velocities[m].dot(velocities[m]);
-    }
-    // center of mass kinetic energy of pairs
-    for (const auto &p: pairsList){
-        double m1 = masses[p.first];
-        double m2 = masses[p.second];
-        Vec3 v1 = velocities[p.first];
-        Vec3 v2 = velocities[p.second];
-        double invMass = 1.0 / (m1 + m2);
-        double redMass = m1 * m2 * invMass;
-        double fracM1 = m1 * invMass;
-        double fracM2 = m2 * invMass;
-        Vec3 comVelocity = fracM1 * v1 + fracM2 * v2;
-        Vec3 relVelocity = v2 - v1;
-
-        comKE += 0.5 * (m1 + m2) * comVelocity.dot(comVelocity);
-        relKE += 0.5 * redMass * relVelocity.dot(relVelocity);
-    }
-    // We ignore the downloadValue argument here and always return the correct value
-    return {comKE, relKE};
-}
-
-
-void ReferenceNoseHooverChainKernel::scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactors) {
-    const auto& atoms = noseHooverChain.getThermostatedAtoms();
-    const auto& pairs = noseHooverChain.getThermostatedPairs();
-    std::vector<Vec3>& velocities = extractVelocities(context);
-    double absScale = scaleFactors.first;
-    double relScale = scaleFactors.second;
-
-    const System& system = context.getSystem();
-    int numParticles = system.getNumParticles();
-    std::vector<double> masses(numParticles);
-    for (int i = 0; i < numParticles; ++i)
-        masses[i] = system.getParticleMass(i);
-    // scale absolute velocities
-    for (const auto &a: atoms){
-        velocities[a] *= absScale;
-    }
-    // scale relative velocities and absolute center of mass velocities for each pair
-    for (const auto &p: pairs){
-        int p1 = p.first;
-        int p2 = p.second;
-        double m1 = masses[p.first];
-        double m2 = masses[p.second];
-        Vec3 v1 = velocities[p.first];
-        Vec3 v2 = velocities[p.second];
-        double invMass = 1.0 / (m1 + m2);
-        double fracM1 = m1 * invMass;
-        double fracM2 = m2 * invMass;
-        Vec3 comVelocity = fracM1 * v1 + fracM2 * v2;
-        Vec3 relVelocity = v2 - v1;
-        velocities[p1] = absScale * comVelocity - relScale * relVelocity * fracM2;
-        velocities[p2] = absScale * comVelocity + relScale * relVelocity * fracM1;
-    }
-}
-
-
 ReferenceApplyMonteCarloBarostatKernel::~ReferenceApplyMonteCarloBarostatKernel() {
     if (barostat)
         delete barostat;