Merge pull request #3 from openmm/master

Update from upstream

libraries/lepton/include/lepton/ParsedExpression.h

@@ -116,6 +116,7 @@ private:
     static ExpressionTreeNode precalculateConstantSubexpressions(const ExpressionTreeNode& node);
     static ExpressionTreeNode substituteSimplerExpression(const ExpressionTreeNode& node);
     static ExpressionTreeNode differentiate(const ExpressionTreeNode& node, const std::string& variable);
+    static bool isConstant(const ExpressionTreeNode& node);
     static double getConstantValue(const ExpressionTreeNode& node);
     static ExpressionTreeNode renameNodeVariables(const ExpressionTreeNode& node, const std::map<std::string, std::string>& replacements);
     ExpressionTreeNode rootNode;

libraries/lepton/src/ParsedExpression.cpp

@@ -121,19 +121,33 @@ ExpressionTreeNode ParsedExpression::substituteSimplerExpression(const Expressio
     vector<ExpressionTreeNode> children(node.getChildren().size());
     for (int i = 0; i < (int) children.size(); i++)
         children[i] = substituteSimplerExpression(node.getChildren()[i]);
+    
+    // Collect some info on constant expressions in children
+    bool first_const = children.size() > 0 && isConstant(children[0]); // is first child constant?
+    bool second_const = children.size() > 1 && isConstant(children[1]); ; // is second child constant?   
+    double first, second; // if yes, value of first and second child
+    if (first_const)
+        first = getConstantValue(children[0]);
+    if (second_const)
+        second = getConstantValue(children[1]);
+
     switch (node.getOperation().getId()) {
         case Operation::ADD:
         {
-            double first = getConstantValue(children[0]);
-            double second = getConstantValue(children[1]);
-            if (first == 0.0) // Add 0
-                return children[1];
-            if (second == 0.0) // Add 0
-                return children[0];
-            if (first == first) // Add a constant
-                return ExpressionTreeNode(new Operation::AddConstant(first), children[1]);
-            if (second == second) // Add a constant
-                return ExpressionTreeNode(new Operation::AddConstant(second), children[0]);
+            if (first_const) {
+                if (first == 0.0) { // Add 0
+                    return children[1];
+                } else { // Add a constant
+                    return ExpressionTreeNode(new Operation::AddConstant(first), children[1]);
+                }
+            }
+            if (second_const) {
+                if (second == 0.0) { // Add 0
+                    return children[0];
+                } else { // Add a constant
+                    return ExpressionTreeNode(new Operation::AddConstant(second), children[0]);
+                }
+            }
             if (children[1].getOperation().getId() == Operation::NEGATE) // a+(-b) = a-b
                 return ExpressionTreeNode(new Operation::Subtract(), children[0], children[1].getChildren()[0]);
             if (children[0].getOperation().getId() == Operation::NEGATE) // (-a)+b = b-a
@@ -144,34 +158,35 @@ ExpressionTreeNode ParsedExpression::substituteSimplerExpression(const Expressio
         {
             if (children[0] == children[1])
                 return ExpressionTreeNode(new Operation::Constant(0.0)); // Subtracting anything from itself is 0
-            double first = getConstantValue(children[0]);
-            if (first == 0.0) // Subtract from 0
-                return ExpressionTreeNode(new Operation::Negate(), children[1]);
-            double second = getConstantValue(children[1]);
-            if (second == 0.0) // Subtract 0
-                return children[0];
-            if (second == second) // Subtract a constant
-                return ExpressionTreeNode(new Operation::AddConstant(-second), children[0]);
+            if (first_const) {
+                if (first == 0.0) // Subtract from 0
+                    return ExpressionTreeNode(new Operation::Negate(), children[1]);
+            }
+            if (second_const) {
+                if (second == 0.0) { // Subtract 0
+                    return children[0];
+                } else { // Subtract a constant
+                    return ExpressionTreeNode(new Operation::AddConstant(-second), children[0]);
+                }
+            }
             if (children[1].getOperation().getId() == Operation::NEGATE) // a-(-b) = a+b
                 return ExpressionTreeNode(new Operation::Add(), children[0], children[1].getChildren()[0]);
             break;
         }
         case Operation::MULTIPLY:
-        {
-            double first = getConstantValue(children[0]);
-            double second = getConstantValue(children[1]);
-            if (first == 0.0 || second == 0.0) // Multiply by 0
+        {   
+            if ((first_const && first == 0.0) || (second_const && second == 0.0)) // Multiply by 0
                 return ExpressionTreeNode(new Operation::Constant(0.0));
-            if (first == 1.0) // Multiply by 1
+            if (first_const && first == 1.0) // Multiply by 1
                 return children[1];
-            if (second == 1.0) // Multiply by 1
+            if (second_const && second == 1.0) // Multiply by 1
                 return children[0];
-            if (children[0].getOperation().getId() == Operation::CONSTANT) { // Multiply by a constant
+            if (first_const) { // Multiply by a constant
                 if (children[1].getOperation().getId() == Operation::MULTIPLY_CONSTANT) // Combine two multiplies into a single one
                     return ExpressionTreeNode(new Operation::MultiplyConstant(first*dynamic_cast<const Operation::MultiplyConstant*>(&children[1].getOperation())->getValue()), children[1].getChildren()[0]);
                 return ExpressionTreeNode(new Operation::MultiplyConstant(first), children[1]);
             }
-            if (children[1].getOperation().getId() == Operation::CONSTANT) { // Multiply by a constant
+            if (second_const) { // Multiply by a constant
                 if (children[0].getOperation().getId() == Operation::MULTIPLY_CONSTANT) // Combine two multiplies into a single one
                     return ExpressionTreeNode(new Operation::MultiplyConstant(second*dynamic_cast<const Operation::MultiplyConstant*>(&children[0].getOperation())->getValue()), children[0].getChildren()[0]);
                 return ExpressionTreeNode(new Operation::MultiplyConstant(second), children[0]);
@@ -202,18 +217,16 @@ ExpressionTreeNode ParsedExpression::substituteSimplerExpression(const Expressio
         {
             if (children[0] == children[1])
                 return ExpressionTreeNode(new Operation::Constant(1.0)); // Dividing anything from itself is 0
-            double numerator = getConstantValue(children[0]);
-            if (numerator == 0.0) // 0 divided by something
+            if (first_const && first == 0.0) // 0 divided by something
                 return ExpressionTreeNode(new Operation::Constant(0.0));
-            if (numerator == 1.0) // 1 divided by something
+            if (first_const && first == 1.0) // 1 divided by something
                 return ExpressionTreeNode(new Operation::Reciprocal(), children[1]);
-            double denominator = getConstantValue(children[1]);
-            if (denominator == 1.0) // Divide by 1
+            if (second_const && second == 1.0) // Divide by 1
                 return children[0];
-            if (children[1].getOperation().getId() == Operation::CONSTANT) {
+            if (second_const) {
                 if (children[0].getOperation().getId() == Operation::MULTIPLY_CONSTANT) // Combine a multiply and a divide into one multiply
-                    return ExpressionTreeNode(new Operation::MultiplyConstant(dynamic_cast<const Operation::MultiplyConstant*>(&children[0].getOperation())->getValue()/denominator), children[0].getChildren()[0]);
-                return ExpressionTreeNode(new Operation::MultiplyConstant(1.0/denominator), children[0]); // Replace a divide with a multiply
+                    return ExpressionTreeNode(new Operation::MultiplyConstant(dynamic_cast<const Operation::MultiplyConstant*>(&children[0].getOperation())->getValue()/second), children[0].getChildren()[0]);
+                return ExpressionTreeNode(new Operation::MultiplyConstant(1.0/second), children[0]); // Replace a divide with a multiply
             }
             if (children[0].getOperation().getId() == Operation::NEGATE && children[1].getOperation().getId() == Operation::NEGATE) // The two negations cancel
                 return ExpressionTreeNode(new Operation::Divide(), children[0].getChildren()[0], children[1].getChildren()[0]);
@@ -229,34 +242,34 @@ ExpressionTreeNode ParsedExpression::substituteSimplerExpression(const Expressio
         }
         case Operation::POWER:
         {
-            double base = getConstantValue(children[0]);
-            if (base == 0.0) // 0 to any power is 0
+            if (first_const && first == 0.0) // 0 to any power is 0
                 return ExpressionTreeNode(new Operation::Constant(0.0));
-            if (base == 1.0) // 1 to any power is 1
-                return ExpressionTreeNode(new Operation::Constant(1.0));
-            double exponent = getConstantValue(children[1]);
-            if (exponent == 0.0) // x^0 = 1
+            if (first_const && first == 1.0) // 1 to any power is 1
                 return ExpressionTreeNode(new Operation::Constant(1.0));
-            if (exponent == 1.0) // x^1 = x
-                return children[0];
-            if (exponent == -1.0) // x^-1 = recip(x)
-                return ExpressionTreeNode(new Operation::Reciprocal(), children[0]);
-            if (exponent == 2.0) // x^2 = square(x)
-                return ExpressionTreeNode(new Operation::Square(), children[0]);
-            if (exponent == 3.0) // x^3 = cube(x)
-                return ExpressionTreeNode(new Operation::Cube(), children[0]);
-            if (exponent == 0.5) // x^0.5 = sqrt(x)
-                return ExpressionTreeNode(new Operation::Sqrt(), children[0]);
-            if (exponent == exponent) // Constant power
-                return ExpressionTreeNode(new Operation::PowerConstant(exponent), children[0]);
+            if (second_const) { // Constant exponent
+                if (second == 0.0) // x^0 = 1
+                    return ExpressionTreeNode(new Operation::Constant(1.0));
+                if (second == 1.0) // x^1 = x
+                    return children[0];
+                if (second == -1.0) // x^-1 = recip(x)
+                    return ExpressionTreeNode(new Operation::Reciprocal(), children[0]);
+                if (second == 2.0) // x^2 = square(x)
+                    return ExpressionTreeNode(new Operation::Square(), children[0]);
+                if (second == 3.0) // x^3 = cube(x)
+                    return ExpressionTreeNode(new Operation::Cube(), children[0]);
+                if (second == 0.5) // x^0.5 = sqrt(x)
+                    return ExpressionTreeNode(new Operation::Sqrt(), children[0]);
+                // Constant power
+                return ExpressionTreeNode(new Operation::PowerConstant(second), children[0]);
+            }
             break;
         }
         case Operation::NEGATE:
         {
             if (children[0].getOperation().getId() == Operation::MULTIPLY_CONSTANT) // Combine a multiply and a negate into a single multiply
                 return ExpressionTreeNode(new Operation::MultiplyConstant(-dynamic_cast<const Operation::MultiplyConstant*>(&children[0].getOperation())->getValue()), children[0].getChildren()[0]);
-            if (children[0].getOperation().getId() == Operation::CONSTANT) // Negate a constant
-                return ExpressionTreeNode(new Operation::Constant(-getConstantValue(children[0])));
+            if (first_const) // Negate a constant
+                return ExpressionTreeNode(new Operation::Constant(-first));
             if (children[0].getOperation().getId() == Operation::NEGATE) // The two negations cancel
                 return children[0].getChildren()[0];
             break;
@@ -265,7 +278,7 @@ ExpressionTreeNode ParsedExpression::substituteSimplerExpression(const Expressio
         {
             if (children[0].getOperation().getId() == Operation::MULTIPLY_CONSTANT) // Combine two multiplies into a single one
                 return ExpressionTreeNode(new Operation::MultiplyConstant(dynamic_cast<const Operation::MultiplyConstant*>(&node.getOperation())->getValue()*dynamic_cast<const Operation::MultiplyConstant*>(&children[0].getOperation())->getValue()), children[0].getChildren()[0]);
-            if (children[0].getOperation().getId() == Operation::CONSTANT) // Multiply two constants
+            if (first_const) // Multiply two constants
                 return ExpressionTreeNode(new Operation::Constant(dynamic_cast<const Operation::MultiplyConstant*>(&node.getOperation())->getValue()*getConstantValue(children[0])));
             if (children[0].getOperation().getId() == Operation::NEGATE) // Combine a multiply and a negate into a single multiply
                 return ExpressionTreeNode(new Operation::MultiplyConstant(-dynamic_cast<const Operation::MultiplyConstant*>(&node.getOperation())->getValue()), children[0].getChildren()[0]);
@@ -303,10 +316,15 @@ ExpressionTreeNode ParsedExpression::differentiate(const ExpressionTreeNode& nod
     return node.getOperation().differentiate(node.getChildren(),childDerivs, variable);
 }
 
+bool ParsedExpression::isConstant(const ExpressionTreeNode& node) {
+    return (node.getOperation().getId() == Operation::CONSTANT);
+}
+
 double ParsedExpression::getConstantValue(const ExpressionTreeNode& node) {
-    if (node.getOperation().getId() == Operation::CONSTANT)
-        return dynamic_cast<const Operation::Constant&>(node.getOperation()).getValue();
-    return numeric_limits<double>::quiet_NaN();
+    if (node.getOperation().getId() != Operation::CONSTANT) {
+        throw Exception("getConstantValue called on a non-constant ExpressionNode");
+    }
+    return dynamic_cast<const Operation::Constant&>(node.getOperation()).getValue();
 }
 
 ExpressionProgram ParsedExpression::createProgram() const {

olla/include/openmm/kernels.h

@@ -64,6 +64,8 @@
 #include "openmm/VariableLangevinIntegrator.h"
 #include "openmm/VariableVerletIntegrator.h"
 #include "openmm/VerletIntegrator.h"
+#include "openmm/NoseHooverIntegrator.h"
+#include "openmm/NoseHooverChain.h"
 #include <iosfwd>
 #include <set>
 #include <string>
@@ -1058,6 +1060,41 @@ public:
     virtual double computeKineticEnergy(ContextImpl& context, const VerletIntegrator& integrator) = 0;
 };
 
+/**
+ * This kernel is invoked by NoseHooverIntegrator to take one time step.
+ */
+class IntegrateVelocityVerletStepKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "IntegrateVelocityVerletStep";
+    }
+    IntegrateVelocityVerletStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param integrator the NoseHooverIntegrator this kernel will be used for
+     */
+    virtual void initialize(const System& system, const NoseHooverIntegrator& integrator) = 0;
+    /**
+     * Execute the kernel.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the NoseHooverIntegrator this kernel is being used for
+     * @param forcesAreValid a reference to the parent integrator's boolean for keeping
+     *                       track of the validity of the current forces.
+     */
+    virtual void execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) = 0;
+    /**
+     * Compute the kinetic energy.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the NoseHooverIntegrator this kernel is being used for
+     */
+    virtual double computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator) = 0;
+};
+
 /**
  * This kernel is invoked by LangevinIntegrator to take one time step.
  */
@@ -1327,6 +1364,58 @@ 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.h

@@ -75,6 +75,8 @@
 #include "openmm/VariableVerletIntegrator.h"
 #include "openmm/Vec3.h"
 #include "openmm/VerletIntegrator.h"
+#include "openmm/NoseHooverIntegrator.h"
+#include "openmm/NoseHooverChain.h"
 #include "openmm/VirtualSite.h"
 #include "openmm/Platform.h"
 #include "openmm/serialization/XmlSerializer.h"

openmmapi/include/openmm/NoseHooverChain.h

+#ifndef OPENMM_NOSEHOOVERCHAIN_H_
+#define OPENMM_NOSEHOOVERCHAIN_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   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 "openmm/Force.h"
+#include <string>
+#include <vector>
+#include "openmm/OpenMMException.h"
+#include "internal/windowsExport.h"
+
+namespace OpenMM {
+
+/**
+ * This class defines a chain of Nose-Hoover particles to be used as a heat bath to
+ * scale the velocities of a collection of particles subject to thermostating.  The
+ * heat bath is propagated using the multi time step approach detailed in
+ *
+ * G. J. Martyna, M. E. Tuckerman, D. J. Tobias and M. L. Klein, Mol. Phys. 87, 1117 (1996).
+ *
+ * where the total number of timesteps used to propagate the chain in each step is
+ * the number of MTS steps multiplied by the number of terms in the Yoshida-Suzuki decomposition.
+ *
+ * Two types of NHC may be created.  The first is a simple thermostat that couples with a given subset
+ * of the atoms within a system, controling their absolute motion.  The second is more elaborate and
+ * can thermostat tethered pairs of atoms and in this case two thermostats are created: one that controls
+ * the absolute center of mass velocity of each pair and another that controls their motion relative to
+ * one another.
+ */
+
+class OPENMM_EXPORT NoseHooverChain {
+public:
+    /**
+     * Create a NoseHooverChain.
+     * 
+     * @param temperature                the temperature of the heat bath for absolute motion (in Kelvin)
+     * @param collisionFrequency         the collision frequency for absolute motion (in 1/ps)
+     * @param relativeTemperature        the temperature of the heat bath for relative motion(in Kelvin).
+     *                                   This is only used if the list of thermostated pairs is not empty.
+     * @param relativeCollisionFrequency the collision frequency for relative motion(in 1/ps).
+     *                                   This is only used if the list of thermostated pairs is not empty.
+     * @param numDOFs                    the number of degrees of freedom in the particles that
+     *                                   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 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
+     * @param thermostatedPairs                 the list of connected pairs to be thermostated; their absolute center of mass motion will
+     *                                          be thermostated independently from their motion relative to one another.
+     */
+    NoseHooverChain(double temperature, double relativeTemperature, double collisionFrequency,
+                    double relativeCollisionFrequency, int numDOFs, int chainLength,
+                    int numMTS, int numYoshidaSuzuki, int chainID,
+                    const std::vector<int>& thermostatedAtoms, const std::vector< std::pair< int, int > > &thermostatedPairs);
+    /**
+     * Get the temperature of the heat bath for treating absolute particle motion (in Kelvin).
+     *
+     * @return the temperature of the heat bath, measured in Kelvin.
+     */
+    double getTemperature() const {
+        return temp;
+    }
+    /**
+     * Set the temperature of the heat bath for treating absolute particle motion.
+     * This will affect any new Contexts you create, but not ones that already exist.
+     *
+     * @param temperature the temperature of the heat bath (in Kelvin)
+     */
+    void setTemperature(double temperature) {
+        temp = temperature;
+    }
+    /**
+     * Get the temperature of the heat bath for treating relative particle motion (in Kelvin).
+     *
+     * @return the temperature of the heat bath, measured in Kelvin.
+     */
+    double getRelativeTemperature() const {
+        return relativeTemp;
+    }
+    /**
+     * Set the temperature of the heat bath for treating relative motion if this thermostat has
+     * been set up to treat connected pairs of atoms.  This will affect any new Contexts you create,
+     * but not ones that already exist.
+     *
+     * @param temperature the temperature of the heat bath for relative motion (in Kelvin)
+     */
+    void setRelativeTemperature(double temperature) {
+        relativeTemp = temperature;
+    }
+    /**
+     * Get the collision frequency for treating absolute particle motion (in 1/ps).
+     *
+     * @return the collision frequency, measured in 1/ps.
+     */
+    double getCollisionFrequency() const {
+        return freq;
+    }
+    /**
+     * Set the collision frequency for treating absolute particle motion.
+     * This will affect any new Contexts you create, but not those that already exist.
+     *
+     * @param frequency the collision frequency (in 1/ps)
+     */
+    void setCollisionFrequency(double frequency) {
+        freq = frequency;
+    }
+    /**
+     * Get the collision frequency for treating relative particle motion (in 1/ps).
+     *
+     * @return the collision frequency, measured in 1/ps.
+     */
+    double getRelativeCollisionFrequency() const {
+        return relativeFreq;
+    }
+    /**
+     * Set the collision frequency for treating relative particle motion if this thermostat has
+     * been set up to handle connected pairs of atoms.  This will affect any new Contexts you create,
+     * but not those that already exist.
+     *
+     * @param frequency the collision frequency (in 1/ps)
+     */
+    void setRelativeCollisionFrequency(double frequency) {
+        relativeFreq = frequency;
+    }
+    /**
+     * Get the number of degrees of freedom in the particles controled by this heat bath.
+     *
+     * @return the number of degrees of freedom.
+     */
+    int getNumDegreesOfFreedom() const {
+        return numDOFs;
+    }
+    /**
+     * Set the number of degrees of freedom in the particles controled by this heat bath.
+     * This will affect any new Contexts you create, but not those that already exist.
+     *
+     * @param numDOF  the number of degrees of freedom.
+     */
+    void setNumDegreesOfFreedom(int numDOF) {
+        numDOFs = numDOF;
+    }
+    /**
+     * Get the chain length of this heat bath.
+     *
+     * @return the chain length.
+     */
+    int getChainLength() const {
+        return chainLength;
+    }
+    /**
+     * Get the number of steps used in the multi time step propagation.
+     *
+     * @returns the number of multi time steps.
+     */
+    int getNumMultiTimeSteps() const {
+        return numMTS;
+    }
+    /**
+     * Get the number of steps used in the Yoshida-Suzuki decomposition for
+     * multi time step propagation.
+     *
+     * @returns the number of multi time steps in the Yoshida-Suzuki decomposition.
+     */
+    int getNumYoshidaSuzukiTimeSteps() const {
+        return numYS;
+    }
+    /**
+     * Get the chain id used to identify this chain
+     *
+     * @returns the chain id 
+     */
+    int getChainID() const {
+        return chainID;
+    }
+    /**
+     * Get the atom ids of all atoms that are thermostated 
+     *
+     * @returns ids of all atoms that are being handled by this thermostat
+     */
+    const std::vector<int>& getThermostatedAtoms() const {
+        return thermostatedAtoms;
+    }
+    /**
+     * Set list of atoms that are handled by this thermostat
+     *
+     * @param atomIDs 
+     */
+    void setThermostatedAtoms(const std::vector<int>& atomIDs){
+        thermostatedAtoms = atomIDs;
+    }
+    /**
+     * Get the list of any connected pairs to be handled by this thermostat.
+     * If this is a regular thermostat, returns an empty vector.
+     *
+     * @returns list of connected pairs.
+     */
+    const std::vector< std::pair< int, int > >& getThermostatedPairs() const {
+        return thermostatedPairs;
+    }
+    /**
+     * In case this thermostat handles the kinetic energy of Drude particles 
+     * set the atom IDs of all parent atoms. 
+     *
+     * @param pairIDs the list of connected pairs to thermostat.
+     */
+    void setThermostatedPairs(const std::vector< std::pair< int, int > >& pairIDs){
+        thermostatedPairs = pairIDs;
+    }
+    /**
+     * Get the weights used in the Yoshida Suzuki multi time step decomposition (dimensionless) 
+     *
+     * @returns the weights for the Yoshida-Suzuki integration
+     */
+    std::vector<double> getYoshidaSuzukiWeights() const;
+    /**
+     * Returns whether or not this force makes use of periodic boundary
+     * conditions.
+     *
+     * @returns true if force uses PBC and false otherwise
+     */
+    bool usesPeriodicBoundaryConditions() const {
+        return false;
+    }
+private:
+    double temp, freq, relativeTemp, relativeFreq;
+    int numDOFs, chainLength, numMTS, numYS;
+    // The suffix used to distinguish NH chains, e.g. for Drude particles vs. regular particles.
+    int chainID;
+    std::vector<int> thermostatedAtoms;
+    std::vector<std::pair<int, int>> thermostatedPairs;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_NOSEHOOVERCHAIN_H_*/

openmmapi/include/openmm/NoseHooverIntegrator.h

+#ifndef OPENMM_NOSEHOOVERINTEGRATOR_H_
+#define OPENMM_NOSEHOOVERINTEGRATOR_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   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 "Integrator.h"
+#include "openmm/State.h"
+#include "openmm/Kernel.h"
+#include "NoseHooverChain.h"
+#include "internal/windowsExport.h"
+
+#include <tuple>
+
+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.
+ */
+
+class OPENMM_EXPORT NoseHooverIntegrator : public Integrator {
+public:
+    /**
+     * Create a NoseHooverIntegrator.  This version creates a bare velocity Verlet 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)
+     */
+    explicit NoseHooverIntegrator(double stepSize);
+    /**
+     * Create a NoseHooverIntegrator.
+     *
+     * @param temperature the target temperature for the system (in Kelvin).
+     * @param collisionFrequency the frequency of the interaction with the heat bath (in inverse picoseconds).
+     * @param stepSize the step size with which to integrate the system (in picoseconds)
+     * @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).
+     */
+    explicit NoseHooverIntegrator(double temperature, double collisionFrequency, double stepSize,
+                                  int chainLength = 3, int numMTS = 3, int numYoshidaSuzuki = 3);
+
+    virtual ~NoseHooverIntegrator();
+   /**
+     * Advance a simulation through time by taking a series of time steps.
+     * 
+     * @param steps   the number of time steps to take
+     */
+    void step(int steps);
+    /**
+     * Add a simple Nose-Hoover Chain thermostat to control the temperature of the full system
+     *
+     * @param temperature the target temperature for the system.
+     * @param collisionFrequency the frequency of the interaction with the heat bath (in 1/ps).
+     * @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).
+     */
+     int addThermostat(double temperature, double collisionFrequency,
+                       int chainLength, int numMTS, int numYoshidaSuzuki);
+    /**
+     * Add a Nose-Hoover Chain thermostat to control the temperature of a collection of atoms and/or pairs of
+     * connected atoms within the full system.  A list of atoms defining the atoms to be thermostated is
+     * provided and the thermostat will only control members of that list.  Additionally a list of pairs of
+     * connected atoms may be provided; in this case both the center of mass absolute motion of each pair is
+     * controlled as well as their motion relative to each other, which is independently thermostated.
+     * If both the list of thermostated particles and thermostated pairs are empty all particles will be thermostated.
+     *
+     * @param thermostatedParticles list of particle ids to be thermostated.
+     * @param thermostatedPairs a list of pairs of connected atoms whose absolute center of mass motion
+     *        and motion relative to one another will be independently thermostated.
+     * @param temperature the target temperature for each pair's absolute of center of mass motion.
+     * @param collisionFrequency the frequency of the interaction with the heat bath for the
+     *        pairs' center of mass motion (in 1/ps).
+     * @param relativeTemperature the target temperature for each pair's relative motion.
+     * @param relativeCollisionFrequency the frequency of the interaction with the heat bath for the
+     *        pairs' relative motion (in 1/ps).
+     * @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).
+     */
+     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);
+    /**
+     * Get the temperature of the i-th chain for controling absolute particle motion (in Kelvin).
+     * 
+     * @param chainID the index of the Nose-Hoover chain thermostat (default=0).
+     * 
+     * @return the temperature.
+     */
+    double getTemperature(int chainID=0) const;
+    /**
+     * set the (absolute motion) temperature of the i-th chain.
+     *
+     * @param temperature the temperature for the Nose-Hoover chain thermostat (in Kelvin).
+     * @param chainID The id of the Nose-Hoover chain thermostat for which the temperature is set (default=0).
+     */
+    void setTemperature(double temperature, int chainID=0);
+    /**
+     * Get the temperature of the i-th chain for controling pairs' relative particle motion (in Kelvin).
+     *
+     * @param chainID the index of the Nose-Hoover chain thermostat (default=0).
+     *
+     * @return the temperature.
+     */
+    double getRelativeTemperature(int chainID=0) const;
+    /**
+     * set the (relative pair motion) temperature of the i-th chain.
+     *
+     * @param temperature the temperature for the Nose-Hoover chain thermostat (in Kelvin).
+     * @param chainID The id of the Nose-Hoover chain thermostat for which the temperature is set (default=0).
+     */
+    void setRelativeTemperature(double temperature, int chainID=0);
+    /**
+     * Get the collision frequency for absolute motion of the i-th chain (in 1/picosecond).
+     * 
+     * @param chainID the index of the Nose-Hoover chain thermostat (default=0).
+     *
+     * @return the collision frequency.
+     */
+    double getCollisionFrequency(int chainID=0) const;
+    /**
+     * Set the collision frequency for absolute motion of the i-th chain.
+     *
+     * @param frequency the collision frequency in picosecond.
+     * @param chainID the index of the Nose-Hoover chain thermostat (default=0).
+     */
+    void setCollisionFrequency(double frequency, int chainID=0);
+    /**
+     * Get the collision frequency for pairs' relative motion of the i-th chain (in 1/picosecond).
+     *
+     * @param chainID the index of the Nose-Hoover chain thermostat (default=0).
+     *
+     * @return the collision frequency.
+     */
+    double getRelativeCollisionFrequency(int chainID=0) const;
+    /**
+     * Set the collision frequency for pairs' relative motion of the i-th chain.
+     *
+     * @param frequency the collision frequency in picosecond.
+     * @param chainID the index of the Nose-Hoover chain thermostat (default=0).
+     */
+    void setRelativeCollisionFrequency(double frequency, int chainID=0);
+    /**
+     * Compute the total (potential + kinetic) heat bath energy for all heat baths
+     * associated with this integrator, at the current time.
+     */
+    double computeHeatBathEnergy();
+    /**
+     * Get the number of Nose-Hoover chains registered with this integrator.
+     */
+    int getNumThermostats() const {
+        return noseHooverChains.size();
+    }
+    /**
+     * Get the NoseHooverChain thermostat 
+     *
+     * @param chainID the index of the Nose-Hoover chain thermostat (default=0).
+     */
+    const NoseHooverChain& getThermostat(int chainID=0) const ;
+    /**
+     * This will be called by the Context when the user modifies aspects of the context state, such
+     * as positions, velocities, or parameters.  This gives the Integrator a chance to discard cached
+     * information.  This is <i>only</i> called when the user modifies information using methods of the Context
+     * object.  It is <i>not</i> called when a ForceImpl object modifies state information in its updateContextState()
+     * method (unless the ForceImpl calls a Context method to perform the modification).
+     * 
+     * @param changed     this specifies what aspect of the Context was changed
+     */
+    virtual void stateChanged(State::DataType changed) {
+       if (State::Positions == changed) forcesAreValid = false;
+    }
+    /**
+     * Return false, if this integrator was set up with the 'default constructor' that thermostats the whole system,
+     * true otherwise. Required for serialization.
+     */
+    bool hasSubsystemThermostats() const {
+        return hasSubsystemThermostats_;
+    }
+    /**
+     * Gets the maximum distance (in nm) that a connected pair may stray from each other. If zero, there are no
+     * constraints on the intra-pair separation.
+     */
+    double getMaximumPairDistance() const { return maxPairDistance_; }
+    /**
+     * Sets the maximum distance (in nm) that a connected pair may stray from each other, implemented using a hard
+     * wall. If set to zero, the hard wall constraint is omited and the pairs are free to be separated by any distance.
+     */
+    void setMaximumPairDistance(double distance) { maxPairDistance_ = distance; }
+    /**
+     * Get a list of all individual atoms (i.e. not involved in a connected Drude-like pair) in the system.
+     */
+    const std::vector<int> & getAllThermostatedIndividualParticles() const { return allAtoms; }
+    /**
+     * Get a list of all connected Drude-like pairs, and their target relative temperature, in the system.
+     */
+    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.
+     * It will also get called again if the application calls reinitialize() on the Context.
+     */
+    void initialize(ContextImpl& context);
+    /**
+     * Goes through the list of thermostats, sets the number of DOFs, and checks for errors in the thermostats.
+     */
+    void initializeThermostats(const System& system);
+    /**
+     * This will be called by the Context when it is destroyed to let the Integrator do any necessary
+     * cleanup.  It will also get called again if the application calls reinitialize() on the Context.
+     */
+    void cleanup();
+    /**
+     * Get the names of all Kernels used by this Integrator.
+     */
+    std::vector<std::string> getKernelNames();
+    /**
+     * Compute the kinetic energy of the system at the current time.
+     */
+    virtual double computeKineticEnergy();
+
+    std::vector<NoseHooverChain> noseHooverChains;
+    std::vector<int> allAtoms;
+    std::vector<std::tuple<int, int, double>> allPairs;
+    bool forcesAreValid;
+    Kernel vvKernel, nhcKernel;
+    bool hasSubsystemThermostats_;
+    double maxPairDistance_;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_NOSEHOOVERINTEGRATOR_H_*/

openmmapi/src/ContextImpl.cpp

@@ -481,6 +481,10 @@ void ContextImpl::loadCheckpoint(istream& stream) {
     }
     updateStateDataKernel.getAs<UpdateStateDataKernel>().loadCheckpoint(*this, stream);
     hasSetPositions = true;
+    integrator.stateChanged(State::Positions);
+    integrator.stateChanged(State::Velocities);
+    integrator.stateChanged(State::Parameters);
+    integrator.stateChanged(State::Energy);
 }
 
 void ContextImpl::systemChanged() {

openmmapi/src/NoseHooverChain.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 "openmm/NoseHooverChain.h"
+#include "openmm/OpenMMException.h"
+
+using namespace OpenMM;
+
+NoseHooverChain::NoseHooverChain(double temperature, double relativeTemperature, double collisionFrequency,
+                                 double relativeCollisionFrequency,
+                                 int numDOFs_, int chainLength_, int numMTS_,
+                                 int numYoshidaSuzuki, int chainID_,
+                                 const std::vector<int>& thermostatedAtoms, const std::vector<std::pair<int, int> > &thermostatedPairs):
+        temp(temperature), relativeTemp(relativeTemperature), freq(collisionFrequency),
+        relativeFreq(relativeCollisionFrequency), numDOFs(numDOFs_),
+        chainLength(chainLength_), numMTS(numMTS_), numYS(numYoshidaSuzuki),
+        chainID(chainID_), thermostatedAtoms(thermostatedAtoms), thermostatedPairs(thermostatedPairs)
+{}
+
+
+std::vector<double> NoseHooverChain::getYoshidaSuzukiWeights() const {
+    switch (numYS) {
+        case 1:
+            return {1};
+        case 3:
+            return {0.828981543588751, -0.657963087177502, 0.828981543588751};
+        case 5:
+            return {0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065,
+                    0.2967324292201065};
+        default:
+            throw OpenMMException("The number of Yoshida-Suzuki weights must be 1,3, or 5.");
+    }
+}

openmmapi/src/NoseHooverIntegrator.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 "openmm/NoseHooverIntegrator.h"
+#include "openmm/Context.h"
+#include "openmm/Force.h"
+#include "openmm/System.h"
+#include "openmm/NoseHooverChain.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/CMMotionRemover.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/kernels.h"
+#include <iostream>
+#include <string>
+#include <algorithm>
+
+using namespace OpenMM;
+using std::string;
+using std::vector;
+
+NoseHooverIntegrator::NoseHooverIntegrator(double stepSize):
+    forcesAreValid(false),
+    hasSubsystemThermostats_(true)
+{
+    setStepSize(stepSize);
+    setConstraintTolerance(1e-5);
+    setMaximumPairDistance(0.0);
+}
+NoseHooverIntegrator::NoseHooverIntegrator(double temperature, double collisionFrequency, double stepSize,
+                                           int chainLength, int numMTS, int numYoshidaSuzuki) : 
+    forcesAreValid(false),
+    hasSubsystemThermostats_(false) {
+
+    setStepSize(stepSize);
+    setConstraintTolerance(1e-5);
+    setMaximumPairDistance(0.0);
+    addThermostat(temperature, collisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
+}
+
+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;
+    return addSubsystemThermostat(std::vector<int>(), std::vector<std::pair<int, int>>(), temperature,
+                                  collisionFrequency, temperature, collisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
+}
+
+int NoseHooverIntegrator::addSubsystemThermostat(const std::vector<int>& thermostatedParticles,
+                                                 const std::vector< std::pair< int, int> > &thermostatedPairs,
+                                                 double temperature, double collisionFrequency,
+                                                 double relativeTemperature, double relativeCollisionFrequency,
+                                                 int chainLength, int numMTS, int numYoshidaSuzuki) {
+    int chainID = noseHooverChains.size();
+    // check if one thermostat already applies to all atoms or pairs
+    if ( (chainID > 0) && (noseHooverChains[0].getThermostatedAtoms().size()*noseHooverChains[0].getThermostatedPairs().size() == 0) ) {
+        throw OpenMMException(
+            "Cannot add thermostat, since one of the thermostats already in the integrator applies to all particles. "
+            "To manually add thermostats, use the constructor that takes only the "
+            "step size as an input argument."
+        );
+    }
+    int nDOF = 0; // is set in initializeThermostats()
+    noseHooverChains.emplace_back(temperature, relativeTemperature,
+                                 collisionFrequency, relativeCollisionFrequency,
+                                 nDOF, chainLength, numMTS,
+                                 numYoshidaSuzuki, chainID,
+                                 thermostatedParticles, thermostatedPairs);
+    return chainID;
+}
+
+const NoseHooverChain& NoseHooverIntegrator::getThermostat(int chainID) const {
+    ASSERT_VALID_INDEX(chainID, noseHooverChains);
+    return noseHooverChains[chainID];
+}
+
+void NoseHooverIntegrator::initializeThermostats(const System &system) {
+
+    allAtoms.clear();
+    allPairs.clear();
+    std::set<int> allAtomsSet;
+    for (int atom = 0; atom < system.getNumParticles(); ++atom) allAtomsSet.insert(atom);
+
+    for (auto &thermostat : noseHooverChains) {
+        const auto &thermostatedParticles = thermostat.getThermostatedAtoms();
+        const auto &thermostatedPairs = thermostat.getThermostatedPairs();
+
+        for( auto& pair : thermostatedPairs ) {
+            allAtomsSet.erase(pair.first);
+            allAtomsSet.erase(pair.second);
+            allPairs.emplace_back(pair.first, pair.second, thermostat.getRelativeTemperature());
+        }
+
+        // figure out the number of DOFs
+        int nDOF = 3*(thermostatedParticles.size() + thermostatedPairs.size());
+        for (int constraintNum = 0; constraintNum < system.getNumConstraints(); constraintNum++) {
+            int particle1, particle2;
+            double distance;
+            system.getConstraintParameters(constraintNum, particle1, particle2, distance);
+            bool particle1_in_thermostatedParticles = ((std::find(thermostatedParticles.begin(),
+                                                                  thermostatedParticles.end(), particle1)
+                                                                    != thermostatedParticles.end())) ||
+                                                      (std::find_if(thermostatedPairs.begin(),
+                                                                    thermostatedPairs.end(),
+                                                                    [&particle1](const std::pair<int, int>& pair){
+                                                                           return pair.first == particle1 || pair.second == particle1;})
+                                                                      != thermostatedPairs.end());
+            bool particle2_in_thermostatedParticles = ((std::find(thermostatedParticles.begin(),
+                                                                  thermostatedParticles.end(), particle2)
+                                                                    != thermostatedParticles.end())) ||
+                                                      (std::find_if(thermostatedPairs.begin(),
+                                                                    thermostatedPairs.end(),
+                                                                    [&particle2](const std::pair<int, int>& pair){
+                                                                           return pair.first == particle2 || pair.second == particle2;})
+                                                                      != thermostatedPairs.end());
+            if ((system.getParticleMass(particle1) > 0) && (system.getParticleMass(particle2) > 0)){
+                if ((particle1_in_thermostatedParticles && !particle2_in_thermostatedParticles) ||
+                     (!particle1_in_thermostatedParticles && particle2_in_thermostatedParticles)){
+                    throw OpenMMException("Cannot add only one of particles " + std::to_string(particle1) + " and " + std::to_string(particle2)
+                                            + " to NoseHooverChain, because they are connected by a constraint.");
+                }
+                if (particle1_in_thermostatedParticles && particle2_in_thermostatedParticles){
+                    nDOF -= 1;
+                }
+            }
+        }
+
+        // remove 3 degrees of freedom from thermostats that act on absolute motions
+        int numForces = system.getNumForces();
+        if (thermostatedPairs.size() == 0){
+            for (int forceNum = 0; forceNum < numForces; ++forceNum) {
+                if (dynamic_cast<const CMMotionRemover*>(&system.getForce(forceNum))) nDOF -= 3;
+            }
+        }
+
+        // set number of DoFs for chain 
+        thermostat.setNumDegreesOfFreedom(nDOF);
+    }
+
+
+    for (int chain1 = 0; chain1 < noseHooverChains.size(); ++chain1){
+        const auto& nhc = noseHooverChains[chain1];
+
+       // make sure that thermostats do not overlap
+        for (int chain2 = 0; chain2 < chain1; ++chain2){
+            const auto& other_nhc = noseHooverChains[chain2];
+            for (auto &particle: nhc.getThermostatedAtoms()){
+                bool isParticleInOtherChain = (std::find(other_nhc.getThermostatedAtoms().begin(),
+                                                         other_nhc.getThermostatedAtoms().end(),
+                                                         particle) != other_nhc.getThermostatedAtoms().end()) ||
+                                              (std::find_if(other_nhc.getThermostatedPairs().begin(),
+                                                            other_nhc.getThermostatedPairs().end(),
+                                               [&particle](const std::pair<int, int>& pair){ return pair.first == particle || pair.second == particle;})
+                                                  != other_nhc.getThermostatedPairs().end());
+                if (isParticleInOtherChain){
+                    throw OpenMMException("Found particle " + std::to_string(particle) + "in a different NoseHooverChain, "
+                                          "but particles can only be thermostated by one thermostat.");
+                }
+            }
+            for (auto &pair: nhc.getThermostatedPairs()){
+                bool isParticleInOtherChain = (std::find(other_nhc.getThermostatedAtoms().begin(),
+                                                         other_nhc.getThermostatedAtoms().end(),
+                                                         pair.first) != other_nhc.getThermostatedAtoms().end()) ||
+                                              (std::find(other_nhc.getThermostatedAtoms().begin(),
+                                                         other_nhc.getThermostatedAtoms().end(),
+                                                         pair.second) != other_nhc.getThermostatedAtoms().end()) ||
+                                              (std::find_if(other_nhc.getThermostatedPairs().begin(),
+                                                            other_nhc.getThermostatedPairs().end(),
+                                               [&pair](const std::pair<int, int>& other_pair){
+                                                        return pair.first == other_pair.first || pair.first == other_pair.second ||
+                                                               pair.second == other_pair.first || pair.second == other_pair.second;})
+                                                  != other_nhc.getThermostatedPairs().end());
+                if (isParticleInOtherChain){
+                    throw OpenMMException("Found pair " + std::to_string(pair.first) + "," +
+                                          std::to_string(pair.second) + " in a different NoseHooverChain, "
+                                          "but particles can only be thermostated by one thermostat.");
+                }
+            }
+        }
+
+        // make sure that massless particles are not thermostated
+        for(auto particle: nhc.getThermostatedAtoms()){
+            double mass = system.getParticleMass(particle);
+            if (mass == 0.0) {
+                throw OpenMMException("Found a particle with no mass (" + std::to_string(particle) + ") in a thermostat. Massless particles cannot be thermostated.");
+            }
+        }
+    }
+
+    for(const auto& atom : allAtomsSet) allAtoms.push_back(atom);
+}
+
+double NoseHooverIntegrator::getTemperature(int chainID) const {
+    ASSERT_VALID_INDEX(chainID, noseHooverChains);
+    return noseHooverChains[chainID].getTemperature();
+}
+
+void NoseHooverIntegrator::setTemperature(double temperature, int chainID){
+    ASSERT_VALID_INDEX(chainID, noseHooverChains);
+    noseHooverChains[chainID].setTemperature(temperature);
+
+}
+
+double NoseHooverIntegrator::getRelativeTemperature(int chainID) const {
+    ASSERT_VALID_INDEX(chainID, noseHooverChains);
+    return noseHooverChains[chainID].getRelativeTemperature();
+}
+
+void NoseHooverIntegrator::setRelativeTemperature(double temperature, int chainID){
+    ASSERT_VALID_INDEX(chainID, noseHooverChains);
+    noseHooverChains[chainID].setRelativeTemperature(temperature);
+
+}
+
+double NoseHooverIntegrator::getCollisionFrequency(int chainID) const {
+    ASSERT_VALID_INDEX(chainID, noseHooverChains);
+    return noseHooverChains[chainID].getCollisionFrequency();
+}
+
+void NoseHooverIntegrator::setCollisionFrequency(double frequency, int chainID){
+    ASSERT_VALID_INDEX(chainID, noseHooverChains);
+    noseHooverChains[chainID].setCollisionFrequency(frequency);
+}
+
+double NoseHooverIntegrator::getRelativeCollisionFrequency(int chainID) const {
+    ASSERT_VALID_INDEX(chainID, noseHooverChains);
+    return noseHooverChains[chainID].getRelativeCollisionFrequency();
+}
+
+void NoseHooverIntegrator::setRelativeCollisionFrequency(double frequency, int chainID){
+    ASSERT_VALID_INDEX(chainID, noseHooverChains);
+    noseHooverChains[chainID].setRelativeCollisionFrequency(frequency);
+}
+
+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;
+        }
+    } else {
+        kE = vvKernel.getAs<IntegrateVelocityVerletStepKernel>().computeKineticEnergy(*context, *this);
+    }
+    return kE;
+}
+
+double NoseHooverIntegrator::computeHeatBathEnergy() {
+    double energy = 0;
+    for(auto &nhc : noseHooverChains) {
+        if (context && (nhc.getNumDegreesOfFreedom() > 0)) {
+            energy += nhcKernel.getAs<NoseHooverChainKernel>().computeHeatBathEnergy(*context, nhc);
+        }
+    }
+    return energy;
+}
+
+void NoseHooverIntegrator::initialize(ContextImpl& contextRef) {
+    if (owner != NULL && &contextRef.getOwner() != owner)
+        throw OpenMMException("This Integrator is already bound to a context");
+
+    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();
+    forcesAreValid = false;
+
+    // check for drude particles and build the Nose-Hoover Chains
+    for (auto& thermostat: noseHooverChains){
+        // if there are no thermostated particles or pairs in the lists this is a regular thermostat for the whole (non-Drude) system
+        if ( (thermostat.getThermostatedAtoms().size() == 0) && (thermostat.getThermostatedPairs().size() == 0) ){
+            std::vector<int> thermostatedParticles;
+            for(int particle = 0; particle < system.getNumParticles(); ++particle) {
+                double mass = system.getParticleMass(particle);
+                if ( (mass > 0) && (mass < 0.8) ){
+                    std::cout << "Warning: Found particles with mass between 0.0 and 0.8 dalton. Did you mean to make a DrudeNoseHooverIntegrator instead? "
+                                 "The thermostat you are about to use will not treat these particles as Drude particles!" << std::endl;
+                }
+                if(system.getParticleMass(particle) > 0) {
+                   thermostatedParticles.push_back(particle);
+                }
+            }
+            thermostat.setThermostatedAtoms(thermostatedParticles);
+        }
+    }
+
+    initializeThermostats(system);
+}
+
+void NoseHooverIntegrator::cleanup() {
+    vvKernel = Kernel();
+    nhcKernel = Kernel();
+}
+
+vector<string> NoseHooverIntegrator::getKernelNames() {
+    std::vector<std::string> names;
+    names.push_back(NoseHooverChainKernel::Name());
+    names.push_back(IntegrateVelocityVerletStepKernel::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) {
+        context->updateContextState();
+        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);
+        }
+    }
+}

platforms/cuda/include/CudaIntegrationUtilities.h

@@ -30,6 +30,7 @@
 #include "openmm/System.h"
 #include "CudaContext.h"
 #include "windowsExportCuda.h"
+#include <map>
 #include <iosfwd>
 
 namespace OpenMM {
@@ -120,6 +121,12 @@ public:
      * @param timeShift   the amount by which to shift the velocities in time
      */
     double computeKineticEnergy(double timeShift);
+    /**
+     * Get the data structure that holds the state of all Nose-Hoover chains
+     * 
+     * @return  vector of chain states
+     */
+    std::map<int, CudaArray>& getNoseHooverChainState();
 private:
     void applyConstraints(bool constrainVelocities, double tol);
     CudaContext& context;
@@ -168,6 +175,7 @@ private:
     double2 lastStepSize;
     struct ShakeCluster;
     struct ConstraintOrderer;
+    std::map<int, CudaArray> noseHooverChainState;
 };
 
 } // namespace OpenMM

platforms/cuda/include/CudaKernels.h

@@ -1366,6 +1366,44 @@ private:
     CUfunction kernel1, kernel2;
 };
 
+/*
+ * This kernel is invoked by NoseHooverIntegrator to take one time step.
+ */
+class CudaIntegrateVelocityVerletStepKernel : public IntegrateVelocityVerletStepKernel {
+public:
+    CudaIntegrateVelocityVerletStepKernel(std::string name, const Platform& platform, CudaContext& cu) :
+                                  IntegrateVelocityVerletStepKernel(name, platform), cu(cu) { }
+    ~CudaIntegrateVelocityVerletStepKernel() {}
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param integrator the NoseHooverIntegrator this kernel will be used for
+     */
+    void initialize(const System& system, const NoseHooverIntegrator& integrator);
+    /**
+     * Execute the kernel.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the VerletIntegrator this kernel is being used for
+     * @param forcesAreValid a reference to the parent integrator's boolean for keeping
+     *                       track of the validity of the current forces.
+     */
+    void execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid);
+    /**
+     * Compute the kinetic energy.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the NoseHooverIntegrator this kernel is being used for
+     */
+    double computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator);
+private:
+    CudaContext& cu;
+    float prevMaxPairDistance;
+    CudaArray maxPairDistanceBuffer, pairListBuffer, atomListBuffer, pairTemperatureBuffer;
+    CUfunction kernel1, kernel2, kernel3, kernelHardWall;
+};
+
 /**
  * This kernel is invoked by LangevinIntegrator to take one time step.
  */
@@ -1716,6 +1754,71 @@ private:
     CUfunction kernel;
 };
 
+/**
+ * This kernel is invoked by NoseHooverChain at the start of each time step to adjust the thermostat
+ * and update the associated particle velocities.
+ */
+class CudaNoseHooverChainKernel : public NoseHooverChainKernel {
+public:
+    CudaNoseHooverChainKernel(std::string name, const Platform& platform, CudaContext& cu) : NoseHooverChainKernel(name, platform), cu(cu) {
+    }
+    ~CudaNoseHooverChainKernel() {}
+    /**
+     * 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 kineticEnergies the {absolute, relative} kineticEnergy 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.
+     */
+    virtual std::pair<double, double> propagateChain(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> kineticEnergies, 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 scaleFactors the {absolute, relative} multiplicative factor by which velocities are scaled.
+     */
+    virtual void scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactors);
+
+private:
+    int sumWorkGroupSize;
+    CudaContext& cu;
+    CudaArray energyBuffer, scaleFactorBuffer, kineticEnergyBuffer, chainMasses, chainForces, heatBathEnergy;
+    std::map<int, CudaArray> atomlists, pairlists;
+    std::map<int, CUfunction> propagateKernels;
+    CUfunction reduceEnergyKernel;
+    CUfunction computeHeatBathEnergyKernel;
+    CUfunction computeAtomsKineticEnergyKernel;
+    CUfunction computePairsKineticEnergyKernel;
+    CUfunction scaleAtomsVelocitiesKernel;
+    CUfunction scalePairsVelocitiesKernel;
+};
+
 /**
  * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
  */

platforms/cuda/src/CudaIntegrationUtilities.cpp

@@ -707,6 +707,24 @@ int CudaIntegrationUtilities::prepareRandomNumbers(int numValues) {
 }
 
 void CudaIntegrationUtilities::createCheckpoint(ostream& stream) {
+    size_t numChains = noseHooverChainState.size();
+    bool useDouble = context.getUseDoublePrecision() || context.getUseMixedPrecision();
+    stream.write((char*) &numChains, sizeof(size_t));
+    for (auto &chainState: noseHooverChainState){
+        int chainID = chainState.first;
+        size_t chainLength = chainState.second.getSize();
+        stream.write((char*) &chainID, sizeof(int));
+        stream.write((char*) &chainLength, sizeof(size_t));
+        if (useDouble) {
+            vector<double2> stateVec;
+            chainState.second.download(stateVec);
+            stream.write((char*) stateVec.data(), sizeof(double2)*chainLength);
+        } else {
+            vector<float2> stateVec;
+            chainState.second.download(stateVec);
+            stream.write((char*) stateVec.data(), sizeof(float2)*chainLength);
+        }
+    }
     if (!random.isInitialized()) 
         return;
     stream.write((char*) &randomPos, sizeof(int));
@@ -719,6 +737,28 @@ void CudaIntegrationUtilities::createCheckpoint(ostream& stream) {
 }
 
 void CudaIntegrationUtilities::loadCheckpoint(istream& stream) {
+    size_t numChains, chainLength;
+    bool useDouble = context.getUseDoublePrecision() || context.getUseMixedPrecision();
+    stream.read((char*) &numChains, sizeof(size_t));
+    noseHooverChainState.clear();
+    for (size_t i=0; i<numChains; i++){
+        int chainID;
+        stream.read((char*) &chainID, sizeof(int));
+        stream.read((char*) &chainLength, sizeof(size_t));
+        if (useDouble) {
+            noseHooverChainState[chainID] = CudaArray();
+            noseHooverChainState[chainID].initialize<double2>(context, chainLength, "chainState" + std::to_string(chainID));
+            std::vector<double2> stateVec(chainLength);
+            stream.read((char*) &stateVec[0], sizeof(double2)*chainLength);
+            noseHooverChainState[chainID].upload(stateVec);
+        } else {
+            noseHooverChainState[chainID] = CudaArray();
+            noseHooverChainState[chainID].initialize<float2>(context, chainLength, "chainState" + std::to_string(chainID));
+            std::vector<float2> stateVec(chainLength);
+            stream.read((char*) &stateVec[0], sizeof(float2)*chainLength);
+            noseHooverChainState[chainID].upload(stateVec);
+        }
+    }
     if (!random.isInitialized()) 
         return;
     stream.read((char*) &randomPos, sizeof(int));
@@ -768,10 +808,13 @@ double CudaIntegrationUtilities::computeKineticEnergy(double timeShift) {
                 energy += (v.x*v.x+v.y*v.y+v.z*v.z)/v.w;
         }
     }
-    
     // Restore the velocities.
     
     if (timeShift != 0)
         posDelta.copyTo(context.getVelm());
     return 0.5*energy;
 }
+
+std::map<int, CudaArray>& CudaIntegrationUtilities::getNoseHooverChainState(){
+        return noseHooverChainState;
+};

platforms/cuda/src/CudaKernelFactory.cpp

@@ -132,6 +132,10 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaIntegrateCustomStepKernel(name, platform, cu);
     if (name == ApplyAndersenThermostatKernel::Name())
         return new CudaApplyAndersenThermostatKernel(name, platform, cu);
+    if (name == NoseHooverChainKernel::Name())
+        return new CudaNoseHooverChainKernel(name, platform, cu);
+    if (name == IntegrateVelocityVerletStepKernel::Name())
+        return new CudaIntegrateVelocityVerletStepKernel(name, platform, cu);
     if (name == ApplyMonteCarloBarostatKernel::Name())
         return new CudaApplyMonteCarloBarostatKernel(name, platform, cu);
     if (name == RemoveCMMotionKernel::Name())

platforms/cuda/src/CudaKernels.cpp

@@ -56,6 +56,7 @@
 #include <cmath>
 #include <iterator>
 #include <set>
+#include <assert.h>
 
 using namespace OpenMM;
 using namespace std;
@@ -386,7 +387,7 @@ void CudaUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context, cons
 
 void CudaUpdateStateDataKernel::createCheckpoint(ContextImpl& context, ostream& stream) {
     cu.setAsCurrent();
-    int version = 2;
+    int version = 3;
     stream.write((char*) &version, sizeof(int));
     int precision = (cu.getUseDoublePrecision() ? 2 : cu.getUseMixedPrecision() ? 1 : 0);
     stream.write((char*) &precision, sizeof(int));
@@ -418,7 +419,7 @@ void CudaUpdateStateDataKernel::loadCheckpoint(ContextImpl& context, istream& st
     cu.setAsCurrent();
     int version;
     stream.read((char*) &version, sizeof(int));
-    if (version != 2)
+    if (version != 3)
         throw OpenMMException("Checkpoint was created with a different version of OpenMM");
     int precision;
     stream.read((char*) &precision, sizeof(int));
@@ -7075,6 +7076,119 @@ double CudaIntegrateVerletStepKernel::computeKineticEnergy(ContextImpl& context,
     return cu.getIntegrationUtilities().computeKineticEnergy(0.5*integrator.getStepSize());
 }
 
+void CudaIntegrateVelocityVerletStepKernel::initialize(const System& system, const NoseHooverIntegrator& integrator) {
+    cu.getPlatformData().initializeContexts(system);
+    cu.setAsCurrent();
+    map<string, string> defines;
+    defines["BOLTZ"] = cu.doubleToString(BOLTZ);
+    CUmodule module = cu.createModule(CudaKernelSources::velocityVerlet, defines, "");
+    kernel1 = cu.getKernel(module, "integrateVelocityVerletPart1");
+    kernel2 = cu.getKernel(module, "integrateVelocityVerletPart2");
+    kernel3 = cu.getKernel(module, "integrateVelocityVerletPart3");
+    kernelHardWall = cu.getKernel(module, "integrateVelocityVerletHardWall");
+    prevMaxPairDistance = -1.0;
+    maxPairDistanceBuffer.initialize<float>(cu, 1, "maxPairDistanceBuffer");
+}
+
+void CudaIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) {
+    cu.setAsCurrent();
+    CudaIntegrationUtilities& integration = cu.getIntegrationUtilities();
+    int paddedNumAtoms = cu.getPaddedNumAtoms();
+    double dt = integrator.getStepSize();
+    cu.getIntegrationUtilities().setNextStepSize(dt);
+
+    if( !forcesAreValid ) context.calcForcesAndEnergy(true, false);
+
+    const auto& atomList = integrator.getAllThermostatedIndividualParticles();
+    const auto& pairList = integrator.getAllThermostatedPairs();
+    int numAtoms = atomList.size();
+    int numPairs = pairList.size();
+    int numParticles = numAtoms + 2*numPairs;
+    float maxPairDistance = integrator.getMaximumPairDistance();
+    // Make sure atom and pair metadata is uploaded and has the correct dimensions
+    if (prevMaxPairDistance != maxPairDistance) {
+        std::vector<float> tmp(1, maxPairDistance);
+        maxPairDistanceBuffer.upload(tmp);
+        prevMaxPairDistance = maxPairDistance;
+    }
+    if (numAtoms !=0 && (!atomListBuffer.isInitialized() || atomListBuffer.getSize() != numAtoms)) {
+        if(atomListBuffer.isInitialized()) {
+            atomListBuffer.resize(atomList.size());
+        } else {
+            atomListBuffer.initialize<int>(cu, atomList.size(), "atomListBuffer");
+        }
+        atomListBuffer.upload(atomList);
+    }
+    if (numPairs !=0 && (!pairListBuffer.isInitialized() || pairListBuffer.getSize() != numPairs)) {
+        if (pairListBuffer.isInitialized()) {
+            pairListBuffer.resize(pairList.size());
+            pairTemperatureBuffer.resize(pairList.size());
+        } else {
+            pairListBuffer.initialize<int2>(cu, pairList.size(), "pairListBuffer");
+            pairTemperatureBuffer.initialize<float>(cu, pairList.size(), "pairTemperatureBuffer");
+        }
+        std::vector<int2> tmp;
+        std::vector<float> tmp2;
+        for(const auto &pair : pairList) {
+            tmp.push_back(make_int2(std::get<0>(pair), std::get<1>(pair)));
+            tmp2.push_back(std::get<2>(pair));
+        }
+        pairListBuffer.upload(tmp);
+        pairTemperatureBuffer.upload(tmp2);
+    }
+    //// Call the first integration kernel.
+
+    CUdeviceptr posCorrection = (cu.getUseMixedPrecision() ? cu.getPosqCorrection().getDevicePointer() : 0);
+    void* args1[] = {&numAtoms, &numPairs, &paddedNumAtoms, &cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(), &posCorrection,
+            &cu.getVelm().getDevicePointer(), &cu.getForce().getDevicePointer(), &integration.getPosDelta().getDevicePointer(),
+            &atomListBuffer.getDevicePointer(), &pairListBuffer.getDevicePointer()};
+    cu.executeKernel(kernel1, args1, std::max(numAtoms,numPairs), 128);
+
+    //// Apply constraints.
+
+    integration.applyConstraints(integrator.getConstraintTolerance());
+
+    //// Call the second integration kernel.
+    void* args2[] = {&numParticles, &cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(), &posCorrection,
+            &cu.getVelm().getDevicePointer(), &integration.getPosDelta().getDevicePointer()};
+    cu.executeKernel(kernel2, args2, numParticles, 128);
+
+    if (numPairs > 0) {
+        //// Enforce hard wall constraint
+        void* argsHardWall[] = {&numPairs, &maxPairDistanceBuffer.getDevicePointer(), 
+                 &cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(), &posCorrection,
+                 &cu.getVelm().getDevicePointer(), &pairListBuffer.getDevicePointer(), 
+                 &pairTemperatureBuffer.getDevicePointer()}; 
+        cu.executeKernel(kernelHardWall, argsHardWall, numPairs, 128);
+    }
+
+    integration.computeVirtualSites();
+
+    //// Update forces
+    context.calcForcesAndEnergy(true, false);
+    forcesAreValid = true;
+
+    //// Call the third integration kernel.
+    void* args3[] = {&numAtoms, &numPairs, &paddedNumAtoms, &cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(), &posCorrection,
+            &cu.getVelm().getDevicePointer(), &cu.getForce().getDevicePointer(), &integration.getPosDelta().getDevicePointer(),
+            &atomListBuffer.getDevicePointer(), &pairListBuffer.getDevicePointer()};
+    cu.executeKernel(kernel3, args3, std::max(numAtoms,numPairs), 128);
+
+    // TODO: Figure out if this is really needed.  The constraint velocities are accounted for
+    // in a finite difference sense in the step 3 kernel, when the velocities are updated.
+    integration.applyVelocityConstraints(integrator.getConstraintTolerance());
+
+    //// Update the time and step count.
+
+    cu.setTime(cu.getTime()+dt);
+    cu.setStepCount(cu.getStepCount()+1);
+    cu.reorderAtoms();
+}
+
+double CudaIntegrateVelocityVerletStepKernel::computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator) {
+    return cu.getIntegrationUtilities().computeKineticEnergy(0);
+}
+
 void CudaIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
     cu.getPlatformData().initializeContexts(system);
     cu.setAsCurrent();
@@ -8364,6 +8478,355 @@ void CudaApplyAndersenThermostatKernel::execute(ContextImpl& context) {
     cu.executeKernel(kernel, args, cu.getNumAtoms());
 }
 
+void CudaNoseHooverChainKernel::initialize() {
+    cu.setAsCurrent();
+
+    bool useDouble = cu.getUseDoublePrecision() || cu.getUseMixedPrecision();
+
+    map<string, string> defines;
+    sumWorkGroupSize = 512;
+    defines["WORK_GROUP_SIZE"] = cu.intToString(sumWorkGroupSize);
+    defines["MIXEDEXP"] = useDouble ? "exp" : "expf";
+    defines["BEGIN_YS_LOOP"] = "for(const real & ys : {1}) {";
+    defines["END_YS_LOOP"] = "}";
+    CUmodule module = cu.createModule(CudaKernelSources::noseHooverChain, defines, "-std=c++11");
+    propagateKernels[1] = cu.getKernel(module, "propagateNoseHooverChain");
+    defines["BEGIN_YS_LOOP"] = "for(const real & ys : {0.828981543588751, -0.657963087177502, 0.828981543588751}){";
+    module = cu.createModule(CudaKernelSources::noseHooverChain, defines, "-std=c++11");
+    propagateKernels[3] = cu.getKernel(module, "propagateNoseHooverChain");
+    defines["BEGIN_YS_LOOP"] = "for(const real & ys : {0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065, 0.2967324292201065}){";
+    module = cu.createModule(CudaKernelSources::noseHooverChain, defines, "-std=c++11");
+    propagateKernels[5] = cu.getKernel(module, "propagateNoseHooverChain");
+    module = cu.createModule(CudaKernelSources::noseHooverChain, defines, "-std=c++11");
+    reduceEnergyKernel = cu.getKernel(module, "reduceEnergyPair");
+
+    computeHeatBathEnergyKernel = cu.getKernel(module, "computeHeatBathEnergy");
+    computeAtomsKineticEnergyKernel = cu.getKernel(module, "computeAtomsKineticEnergy");
+    computePairsKineticEnergyKernel = cu.getKernel(module, "computePairsKineticEnergy");
+    scaleAtomsVelocitiesKernel = cu.getKernel(module, "scaleAtomsVelocities");
+    scalePairsVelocitiesKernel = cu.getKernel(module, "scalePairsVelocities");
+
+    int energyBufferSize = cu.getEnergyBuffer().getSize();
+    if (cu.getUseDoublePrecision() || cu.getUseMixedPrecision()) {
+        energyBuffer.initialize<double2>(cu, energyBufferSize, "energyBuffer");
+    } else {
+        energyBuffer.initialize<float2>(cu, energyBufferSize, "energyBuffer");
+    }
+}
+
+std::pair<double, double> CudaNoseHooverChainKernel::propagateChain(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> kineticEnergies, double timeStep) {
+    cu.setAsCurrent();
+
+    bool useDouble = cu.getUseDoublePrecision() || cu.getUseMixedPrecision();
+
+    int chainID = nhc.getChainID();
+    int nAtoms = nhc.getThermostatedAtoms().size();
+    int nPairs = nhc.getThermostatedPairs().size();
+    int chainLength = nhc.getChainLength();
+    int numYS = nhc.getNumYoshidaSuzukiTimeSteps();
+    int numMTS = nhc.getNumMultiTimeSteps();
+    int numDOFs = nhc.getNumDegreesOfFreedom();
+    double temperature = nhc.getTemperature();
+    double frequency = nhc.getCollisionFrequency();
+    double relativeTemperature = nhc.getRelativeTemperature();
+    double relativeFrequency = nhc.getRelativeCollisionFrequency();
+
+    if (numYS != 1 && numYS != 3 && numYS != 5) {
+        throw OpenMMException("Number of Yoshida Suzuki time steps has to be 1, 3, or 5.");
+    }
+
+    auto & chainState = cu.getIntegrationUtilities().getNoseHooverChainState();
+
+    if (!scaleFactorBuffer.isInitialized() ||scaleFactorBuffer.getSize() == 0) {
+        if(useDouble){
+            std::vector<double2> zeros{{0,0}};
+            if (scaleFactorBuffer.isInitialized()) {
+                scaleFactorBuffer.resize(1);
+            } else {
+                scaleFactorBuffer.initialize<double2>(cu, 1, "scaleFactorBuffer");
+            }
+            scaleFactorBuffer.upload(zeros);
+        } else {
+            std::vector<float2> zeros{{0,0}};
+            if (scaleFactorBuffer.isInitialized()) {
+                scaleFactorBuffer.resize(1);
+            } else {
+                scaleFactorBuffer.initialize<float2>(cu, 1, "scaleFactorBuffer");
+            }
+            scaleFactorBuffer.upload(zeros);
+        }
+    }
+    std::vector<double> zeros(chainLength,0);
+    if (!chainForces.isInitialized() || !chainMasses.isInitialized() ){
+        if(useDouble){
+            if (chainForces.isInitialized()) {
+                chainMasses.resize(chainLength);
+                chainForces.resize(chainLength);
+            } else {
+                chainMasses.initialize<double>(cu, chainLength, "chainMasses");
+                chainForces.initialize<double>(cu, chainLength, "chainForces");
+            }
+            chainMasses.upload(zeros);
+            chainForces.upload(zeros);
+        } else {
+            if (chainForces.isInitialized()) {
+                chainMasses.resize(chainLength);
+                chainForces.resize(chainLength);
+            } else {
+                chainMasses.initialize<float>(cu, chainLength, "chainMasses");
+                chainForces.initialize<float>(cu, chainLength, "chainForces");
+            }
+            chainMasses.upload(zeros);
+            chainForces.upload(zeros);
+        }
+    }
+    if (chainForces.getSize() < chainLength) chainMasses.resize(chainLength);
+    if (chainMasses.getSize() < chainLength) chainMasses.resize(chainLength);
+
+    float timeStepFloat = (float) timeStep;
+    // N.B. We ignore the incoming kineticEnergy and grab it from the device buffer instead
+    if (nAtoms) {
+        if (!chainState.count(2*chainID))  chainState[2*chainID] = CudaArray();
+        if (chainState.at(2*chainID).getSize() != chainLength) {
+            // We need to upload the CUDA array
+            if(useDouble){
+                if (chainState.at(2*chainID).isInitialized()) {
+                    chainState.at(2*chainID).resize(chainLength);
+                } else {
+                    chainState.at(2*chainID).initialize<double2>(cu, chainLength, "chainState" + std::to_string(2*chainID));
+                }
+                std::vector<double2> zeros(chainLength, make_double2(0, 0));
+                chainState.at(2*chainID).upload(zeros.data());
+            } else {
+                if (chainState.at(2*chainID).isInitialized()) {
+                    chainState.at(2*chainID).resize(chainLength);
+                } else {
+                    chainState.at(2*chainID).initialize<float2>(cu, chainLength, "chainState" + std::to_string(2*chainID));
+                }
+                std::vector<float2> zeros(chainLength, make_float2(0, 0));
+                chainState.at(2*chainID).upload(zeros.data());
+            }
+        }
+        int chainType = 0;
+        double kT = BOLTZ * temperature;
+        float kTfloat = (float) kT;
+        float frequencyFloat = (float) frequency;
+        void *args[] = {
+                    &chainState[2*chainID].getDevicePointer(),
+                    &kineticEnergyBuffer.getDevicePointer(),
+                    &scaleFactorBuffer.getDevicePointer(),
+                    &chainMasses.getDevicePointer(),
+                    &chainForces.getDevicePointer(),
+                    &chainType, &chainLength, &numMTS, &numDOFs,
+                    &timeStepFloat,
+                    useDouble ? (void*) &kT : (void*) &kTfloat, 
+                    &frequencyFloat
+                    };
+        cu.executeKernel(propagateKernels[numYS], args, 1, 1);
+    }
+    if (nPairs) {
+        if (!chainState.count(2*chainID+1)) chainState[2*chainID+1] = CudaArray();
+        if (chainState.at(2*chainID+1).getSize() != chainLength) {
+            // We need to upload the CUDA array
+            if(useDouble){
+                if (chainState.at(2*chainID+1).isInitialized()) {
+                    chainState.at(2*chainID+1).resize(chainLength);
+                } else {
+                    chainState.at(2*chainID+1).initialize<double2>(cu, chainLength, "chainState" + std::to_string(2*chainID+1));
+                }
+                std::vector<double2> zeros(chainLength, make_double2(0, 0));
+                chainState.at(2*chainID+1).upload(zeros.data());
+            } else {
+                if (chainState.at(2*chainID+1).isInitialized()) {
+                    chainState.at(2*chainID+1).resize(chainLength);
+                } else {
+                    chainState.at(2*chainID+1).initialize<float2>(cu, chainLength, "chainState" + std::to_string(2*chainID+1));
+                }
+                std::vector<float2> zeros(chainLength, make_float2(0, 0));
+                chainState.at(2*chainID+1).upload(zeros.data());
+            }
+        }
+        int chainType = 1;
+        double kT = BOLTZ * relativeTemperature;
+        int ndf = 3*nPairs;
+        float kTfloat = (float) kT;
+        float frequencyFloat = (float) relativeFrequency;
+        void *args[] = {
+                    &chainState[2*chainID+1].getDevicePointer(),
+                    &kineticEnergyBuffer.getDevicePointer(),
+                    &scaleFactorBuffer.getDevicePointer(),
+                    &chainMasses.getDevicePointer(),
+                    &chainForces.getDevicePointer(),
+                    &chainType, &chainLength, &numMTS, &ndf,
+                    &timeStepFloat,
+                    useDouble ? (void*) &kT : (void*) &kTfloat, 
+                    &frequencyFloat
+                    };
+        cu.executeKernel(propagateKernels[numYS], args, 1, 1);
+    }
+    return {0, 0};
+}
+
+double CudaNoseHooverChainKernel::computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &nhc) {
+    cu.setAsCurrent();
+
+    bool useDouble = cu.getUseDoublePrecision() || cu.getUseMixedPrecision();
+
+    int chainID = nhc.getChainID();
+    int chainLength = nhc.getChainLength();
+
+    auto & chainState = cu.getIntegrationUtilities().getNoseHooverChainState();
+
+    bool absChainIsValid = chainState.count(2*chainID) != 0 &&
+                           chainState[2*chainID].isInitialized() &&
+                           chainState[2*chainID].getSize() == chainLength;
+    bool relChainIsValid = chainState.count(2*chainID+1) != 0 &&
+                           chainState[2*chainID+1].isInitialized() &&
+                           chainState[2*chainID+1].getSize() == chainLength;
+
+    if (!absChainIsValid && !relChainIsValid) return 0.0;
+
+    if (!heatBathEnergy.isInitialized() || heatBathEnergy.getSize() == 0) {
+        if(useDouble){
+            std::vector<double> one(1);
+            heatBathEnergy.initialize<double>(cu, 1, "heatBathEnergy");
+            heatBathEnergy.upload(one);
+        } else {
+            std::vector<float> one(1);
+            heatBathEnergy.initialize<float>(cu, 1, "heatBathEnergy");
+            heatBathEnergy.upload(one);
+        }
+    }
+
+    cu.clearBuffer(heatBathEnergy);
+
+    if (absChainIsValid) {
+        int numDOFs = nhc.getNumDegreesOfFreedom();
+        double temperature = nhc.getTemperature();
+        double frequency = nhc.getCollisionFrequency();
+        double kT = BOLTZ * temperature;
+        float kTfloat = (float) kT;
+        float frequencyFloat = (float) frequency;
+
+        void * args[] = {&heatBathEnergy.getDevicePointer(), &chainLength, &numDOFs, 
+                         useDouble ? (void*) &kT : (void*) & kTfloat,
+                         &frequencyFloat, &chainState[2*chainID].getDevicePointer()};
+        cu.executeKernel(computeHeatBathEnergyKernel, args, 1, 1);
+    }
+    if (relChainIsValid) {
+        int numDOFs = 3 * nhc.getThermostatedPairs().size();
+        double temperature = nhc.getRelativeTemperature();
+        double frequency = nhc.getRelativeCollisionFrequency();
+        double kT = BOLTZ * temperature;
+        float kTfloat = (float) kT;
+        float frequencyFloat = (float) frequency;
+
+        void * args[] = {&heatBathEnergy.getDevicePointer(), &chainLength, &numDOFs, 
+                         useDouble ? (void*) &kT : (void*) & kTfloat,
+                         &frequencyFloat, &chainState[2*chainID+1].getDevicePointer()};
+        cu.executeKernel(computeHeatBathEnergyKernel, args, 1, 1);
+    }
+
+
+    void * pinnedBuffer = cu.getPinnedBuffer();
+    heatBathEnergy.download(pinnedBuffer);
+    if (useDouble){
+        return *((double*) pinnedBuffer);
+    } else {
+        return *((float*) pinnedBuffer);
+    }
+}
+
+std::pair<double, double> CudaNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &nhc, bool downloadValue) {
+    cu.setAsCurrent();
+
+    bool useDouble = cu.getUseDoublePrecision() || cu.getUseMixedPrecision();
+
+    int chainID = nhc.getChainID();
+    const auto & nhcAtoms = nhc.getThermostatedAtoms();
+    const auto & nhcPairs = nhc.getThermostatedPairs();
+    auto nAtoms = nhcAtoms.size();
+    auto nPairs = nhcPairs.size();
+    if (nAtoms) {
+        if (!atomlists.count(chainID)) { 
+            // We need to upload the CUDA array
+            atomlists[chainID] = CudaArray();
+            atomlists[chainID].initialize<int>(cu, nAtoms, "atomlist" + std::to_string(chainID));
+            atomlists[chainID].upload(nhcAtoms);
+        }
+        if (atomlists[chainID].getSize() != nAtoms) {
+            throw OpenMMException("Number of atoms changed. Cannot be handled by the same Nose-Hoover thermostat.");
+        }
+    }
+    if (nPairs) {
+        if (!pairlists.count(chainID)) { 
+            // We need to upload the CUDA array
+            pairlists[chainID] = CudaArray();
+            pairlists[chainID].initialize<int2>(cu, nPairs, "pairlist" + std::to_string(chainID));
+            std::vector<int2> int2vec;
+            for(const auto &p : nhcPairs) int2vec.push_back(make_int2(p.first, p.second));
+            pairlists[chainID].upload(int2vec);
+        }
+        if (pairlists[chainID].getSize() != nPairs) {
+            throw OpenMMException("Number of thermostated pairs changed. Cannot be handled by the same Nose-Hoover thermostat.");
+        }
+    }
+    if (!kineticEnergyBuffer.isInitialized() || kineticEnergyBuffer.getSize() == 0) {
+        if(useDouble){
+            std::vector<double2> zeros{{0,0}};
+            kineticEnergyBuffer.initialize<double2>(cu, 1, "kineticEnergyBuffer");
+            kineticEnergyBuffer.upload(zeros);
+        } else {
+            std::vector<float2> zeros{{0,0}};
+            kineticEnergyBuffer.initialize<float2>(cu, 1, "kineticEnergyBuffer");
+            kineticEnergyBuffer.upload(zeros);
+        }
+    }
+    cu.clearBuffer(energyBuffer);
+    if (nAtoms) {
+        void *args[] = {&energyBuffer.getDevicePointer(),&nAtoms, &cu.getVelm().getDevicePointer(), &atomlists[chainID].getDevicePointer()};
+        cu.executeKernel(computeAtomsKineticEnergyKernel, args, nAtoms);
+    }
+    if (nPairs) {
+        void *args[] = {&energyBuffer.getDevicePointer(),&nPairs, &cu.getVelm().getDevicePointer(), &pairlists[chainID].getDevicePointer()};
+        cu.executeKernel(computePairsKineticEnergyKernel, args, nPairs);
+    }
+
+    //taken from CudaContext::reduceEnergy(); the final kinetic energy will live in the kineticEnergy buffer
+    int bufferSize = energyBuffer.getSize();
+    void* args2[] = {&energyBuffer.getDevicePointer(), &kineticEnergyBuffer.getDevicePointer(), &bufferSize, &sumWorkGroupSize};
+    cu.executeKernel(reduceEnergyKernel, args2, sumWorkGroupSize, sumWorkGroupSize, sumWorkGroupSize*energyBuffer.getElementSize());
+
+    std::pair<double, double> KEs = {0, 0};
+    if (downloadValue) {
+        void * pinnedBuffer = cu.getPinnedBuffer();
+        kineticEnergyBuffer.download(pinnedBuffer);
+        KEs.first = useDouble ? *((double*) pinnedBuffer) : *((float*) pinnedBuffer);
+        KEs.second = useDouble ? *((double*) pinnedBuffer + 1) : *((float*) pinnedBuffer + 1);
+    }
+    return KEs;
+}
+
+void CudaNoseHooverChainKernel::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.
+    cu.setAsCurrent();
+
+    int chainID = nhc.getChainID();
+    auto nAtoms = nhc.getThermostatedAtoms().size();
+    auto nPairs = nhc.getThermostatedPairs().size();
+    if(nAtoms) {
+        void *args[] = {&scaleFactorBuffer.getDevicePointer(),
+                        &nAtoms, &cu.getVelm().getDevicePointer(), &atomlists[chainID].getDevicePointer()};
+        cu.executeKernel(scaleAtomsVelocitiesKernel, args, nAtoms);
+    }
+    if(nPairs) {
+        void *args[] = {&scaleFactorBuffer.getDevicePointer(),
+                        &nPairs, &cu.getVelm().getDevicePointer(), &pairlists[chainID].getDevicePointer()};
+        cu.executeKernel(scalePairsVelocitiesKernel, args, nPairs);
+    }
+}
+
 void CudaApplyMonteCarloBarostatKernel::initialize(const System& system, const Force& thermostat) {
     cu.setAsCurrent();
     savedPositions.initialize(cu, cu.getPaddedNumAtoms(), cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4), "savedPositions");

platforms/cuda/src/CudaPlatform.cpp

@@ -96,6 +96,7 @@ CudaPlatform::CudaPlatform() {
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
+    registerKernelFactory(IntegrateVelocityVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBAOABStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
@@ -103,6 +104,7 @@ 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/src/kernels/noseHooverChain.cu

+
+#include <initializer_list>
+
+extern "C" __global__ void propagateNoseHooverChain(mixed2* __restrict__ chainData, const mixed2 * __restrict__ energySum, mixed2* __restrict__ scaleFactor,
+                                                    mixed* __restrict__ chainMasses, mixed* __restrict__ chainForces, 
+                                                    int chainType, int chainLength, int numMTS, int numDOFs, float timeStep,
+                                                    mixed kT, float frequency){
+    const mixed & kineticEnergy = chainType ? energySum[0].y : energySum[0].x;
+    mixed &scale = chainType ? scaleFactor[0].y : scaleFactor[0].x;
+    scale = (mixed) 1;
+    if(kineticEnergy < 1e-8) return;
+    for (int bead = 0; bead < chainLength; ++bead) chainMasses[bead] = kT / (frequency * frequency);
+    chainMasses[0] *= numDOFs;
+    mixed KE2 = 2.0f * kineticEnergy;
+    mixed timeOverMTS = timeStep / numMTS;
+    chainForces[0] = (KE2 - numDOFs * kT) / chainMasses[0];
+    for (int bead = 0; bead < chainLength - 1; ++bead) {
+        chainForces[bead + 1] = (chainMasses[bead] * chainData[bead].y * chainData[bead].y - kT) / chainMasses[bead + 1];
+    }
+    for (int mts = 0; mts < numMTS; ++mts) {
+        BEGIN_YS_LOOP
+            mixed wdt = ys * timeOverMTS;
+            chainData[chainLength-1].y += 0.25f * wdt * chainForces[chainLength-1];
+            for (int bead = chainLength - 2; bead >= 0; --bead) {
+                mixed aa = MIXEDEXP(-0.125f * wdt * chainData[bead + 1].y);
+                chainData[bead].y = aa * (chainData[bead].y * aa + 0.25f * wdt * chainForces[bead]);
+            }
+            // update particle velocities
+            mixed aa = MIXEDEXP(-0.5f * 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;
+            }
+            // 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 = MIXEDEXP(-0.125f * wdt * chainData[bead + 1].y);
+                chainData[bead].y = aa * (aa * chainData[bead].y + 0.25f * 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];
+        END_YS_LOOP
+    } // MTS loop
+}
+
+
+/**
+ * Compute total (potential + kinetic) energy of the Nose-Hoover beads
+ */
+extern "C" __global__ void computeHeatBathEnergy(mixed* __restrict__ heatBathEnergy, int chainLength, int numDOFs,
+                                                 mixed kT, float frequency, const mixed2* __restrict__ chainData){
+    // Note that this is always incremented; make sure it's zeroed properly before the first call
+    mixed &energy = heatBathEnergy[0];
+
+    for(int i = 0; i < chainLength; ++i) {
+        mixed prefac = i ? 1 : numDOFs;
+        mixed mass = prefac * kT / (frequency * frequency);
+        mixed velocity = chainData[i].y; 
+        // The kinetic energy of this bead
+        energy += 0.5f * mass * velocity * velocity;
+        // The potential energy of this bead
+        mixed position = chainData[i].x;
+        energy += prefac * kT * position;
+    }
+}
+
+extern "C" __global__ void computeAtomsKineticEnergy(mixed2 * __restrict__ energyBuffer, int numAtoms,
+                                                     const mixed4* __restrict__ velm, const int *__restrict__ atoms){
+    mixed2 energy = make_mixed2(0,0);
+    //energy = 1; return;
+    for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < numAtoms; index += blockDim.x*gridDim.x) {
+        int atom = atoms[index];
+        mixed4 v = velm[atom];
+        mixed mass = v.w == 0 ? 0 : 1 / v.w;
+        energy.x += 0.5f * mass * (v.x*v.x + v.y*v.y + v.z*v.z);
+    }
+    energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] = energy;
+}
+
+extern "C" __global__ void computePairsKineticEnergy(mixed2 * __restrict__ energyBuffer, int numPairs,
+                                                     const mixed4* __restrict__ velm, const int2 *__restrict__ pairs){
+    mixed2 energy = make_mixed2(0,0);
+    for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < numPairs; index += blockDim.x*gridDim.x) {
+        int2 pair = pairs[index];
+        int atom1 = pair.x;
+        int atom2 = pair.y;
+        mixed4 v1 = velm[atom1];
+        mixed4 v2 = velm[atom2];
+        mixed m1 = v1.w == 0 ? 0 : 1 / v1.w;
+        mixed m2 = v2.w == 0 ? 0 : 1 / v2.w;
+        mixed4 cv;
+        cv.x = (m1*v1.x + m2*v2.x) / (m1 + m2);
+        cv.y = (m1*v1.y + m2*v2.y) / (m1 + m2);
+        cv.z = (m1*v1.z + m2*v2.z) / (m1 + m2);
+        mixed4 rv;
+        rv.x = v2.x - v1.x;
+        rv.y = v2.y - v1.y;
+        rv.z = v2.z - v1.z;
+        energy.x += 0.5f * (m1 + m2) * (cv.x*cv.x + cv.y*cv.y + cv.z*cv.z);
+        energy.y += 0.5f * (m1 * m2 / (m1 + m2)) * (rv.x*rv.x + rv.y*rv.y + rv.z*rv.z);
+    }
+    // The atoms version of this has been called already, so accumulate instead of assigning here
+    energyBuffer[blockIdx.x*blockDim.x+threadIdx.x].x += energy.x;
+    energyBuffer[blockIdx.x*blockDim.x+threadIdx.x].y += energy.y;
+}
+
+extern "C" __global__ void scaleAtomsVelocities(mixed2* __restrict__ scaleFactor, int numAtoms,
+                                                mixed4* __restrict__ velm, const int *__restrict__ atoms){
+    const mixed &scale = scaleFactor[0].x;
+    for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < numAtoms; index += blockDim.x*gridDim.x) {
+        int atom = atoms[index];
+        mixed4 &v = velm[atom];
+        v.x *= scale;
+        v.y *= scale;
+        v.z *= scale;
+    }
+}
+
+extern "C" __global__ void scalePairsVelocities(mixed2 * __restrict__ scaleFactor, int numPairs,
+                                                mixed4* __restrict__ velm, const int2 *__restrict__ pairs){
+    const mixed &absScale = scaleFactor[0].x;
+    const mixed &relScale = scaleFactor[0].y;
+    for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < numPairs; index += blockDim.x*gridDim.x) {
+        int atom1 = pairs[index].x;
+        int atom2 = pairs[index].y;
+        mixed4 v1 = velm[atom1];
+        mixed4 v2 = velm[atom2];
+        mixed m1 = v1.w == 0 ? 0 : 1 / v1.w;
+        mixed m2 = v2.w == 0 ? 0 : 1 / v2.w;
+        mixed4 cv;
+        cv.x = (m1*v1.x + m2*v2.x) / (m1 + m2);
+        cv.y = (m1*v1.y + m2*v2.y) / (m1 + m2);
+        cv.z = (m1*v1.z + m2*v2.z) / (m1 + m2);
+        mixed4 rv;
+        rv.x = v2.x - v1.x;
+        rv.y = v2.y - v1.y;
+        rv.z = v2.z - v1.z;
+        v1.x = absScale * cv.x - relScale * rv.x * m2 / (m1 + m2);
+        v1.y = absScale * cv.y - relScale * rv.y * m2 / (m1 + m2);
+        v1.z = absScale * cv.z - relScale * rv.z * m2 / (m1 + m2);
+        v2.x = absScale * cv.x + relScale * rv.x * m1 / (m1 + m2);
+        v2.y = absScale * cv.y + relScale * rv.y * m1 / (m1 + m2);
+        v2.z = absScale * cv.z + relScale * rv.z * m1 / (m1 + m2);
+        velm[atom1] = v1;
+        velm[atom2] = v2;
+    }
+}
+
+/**
+ * Sum the energy buffer containing a pair of energies stored as mixed2.  This is copied from utilities.cu with small modifications
+ */
+extern "C" __global__ void reduceEnergyPair(const mixed2* __restrict__ energyBuffer, mixed2* __restrict__ result, int bufferSize, int workGroupSize) {
+    __shared__ mixed2 tempBuffer[WORK_GROUP_SIZE];
+    const unsigned int thread = threadIdx.x;
+    mixed2 sum = make_mixed2(0,0);
+    for (unsigned int idx = thread; idx < bufferSize; idx += blockDim.x) {
+        sum.x += energyBuffer[idx].x;
+        sum.y += energyBuffer[idx].y;
+    }
+    tempBuffer[thread] = sum;
+    for (int i = 1; i < workGroupSize; i *= 2) {
+        __syncthreads();
+        if (thread%(i*2) == 0 && thread+i < workGroupSize) {
+            tempBuffer[thread].x += tempBuffer[thread+i].x;
+            tempBuffer[thread].y += tempBuffer[thread+i].y;
+        }
+    }
+    if (thread == 0)
+        *result = tempBuffer[0];
+}

platforms/cuda/src/kernels/utilities.cu

@@ -99,4 +99,4 @@ __global__ void setCharges(real* __restrict__ charges, real4* __restrict__ posq,
     for (int i = blockDim.x*blockIdx.x+threadIdx.x; i < numAtoms; i += blockDim.x*gridDim.x)
         posq[i].w = charges[atomOrder[i]];
 }
-}
\ No newline at end of file
+}

platforms/cuda/src/kernels/velocityVerlet.cu

+/**
+ * Perform the first step of Velocity Verlet integration.
+ * 
+ * update displacements (posDelta) and velocities (velm)
+ */
+
+extern "C" __global__ void integrateVelocityVerletPart1(int numAtoms, int numPairs, int paddedNumAtoms, const mixed2* __restrict__ dt, const real4* __restrict__ posq,
+        const real4* __restrict__ posqCorrection, mixed4* __restrict__ velm, const long long* __restrict__ force, mixed4* __restrict__ posDelta,
+        const int* __restrict__ atomList, const int2* __restrict__ pairList) {
+    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;
+    for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < numAtoms; index += blockDim.x*gridDim.x) {
+        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;
+            velm[atom] = velocity;
+        }
+    }
+    for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < numPairs; index += blockDim.x*gridDim.x) {
+        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;
+        comFrc.x = force[atom1] + force[atom2];
+        comFrc.y = force[atom1 + paddedNumAtoms] + force[atom2 + paddedNumAtoms];
+        comFrc.z = force[atom1 + paddedNumAtoms*2] + force[atom2 + paddedNumAtoms*2];
+        mixed3 relFrc;
+        relFrc.x = mass1fract*force[atom2] - mass2fract*force[atom1];
+        relFrc.y = mass1fract*force[atom2+paddedNumAtoms] - mass2fract*force[atom1+paddedNumAtoms];
+        relFrc.z = mass1fract*force[atom2+paddedNumAtoms*2] - mass2fract*force[atom1+paddedNumAtoms*2];
+        comVel.x += comFrc.x * scale * invTotMass;
+        comVel.y += comFrc.y * scale * invTotMass;
+        comVel.z += comFrc.z * scale * invTotMass;
+        relVel.x += relFrc.x * scale * invRedMass;
+        relVel.y += relFrc.y * scale * invRedMass;
+        relVel.z += relFrc.z * scale * 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;
+        }
+    }
+}
+
+/**
+ * Perform the second step of Velocity Verlet integration.
+ * 
+ * apply displacements to positions (posq) after constraints have been enforced
+ */
+
+extern "C" __global__ void integrateVelocityVerletPart2(int numAtoms, mixed2* __restrict__ dt, real4* __restrict__ posq,
+        real4* __restrict__ posqCorrection, mixed4* __restrict__ velm, const mixed4* __restrict__ posDelta) {
+    mixed2 stepSize = dt[0];
+    int index = blockIdx.x*blockDim.x+threadIdx.x;
+    if (index == 0)
+        dt[0].x = stepSize.y;
+    for (; index < numAtoms; index += blockDim.x*gridDim.x) {
+        mixed4 velocity = velm[index];
+        if (velocity.w != 0.0) {
+#ifdef USE_MIXED_PRECISION
+            real4 pos1 = posq[index];
+            real4 pos2 = posqCorrection[index];
+            mixed4 pos = make_mixed4(pos1.x+(mixed)pos2.x, pos1.y+(mixed)pos2.y, pos1.z+(mixed)pos2.z, pos1.w);
+#else
+            real4 pos = posq[index];
+#endif
+            mixed4 delta = posDelta[index];
+            pos.x += delta.x;
+            pos.y += delta.y;
+            pos.z += delta.z;
+#ifdef USE_MIXED_PRECISION
+            posq[index] = make_real4((real) pos.x, (real) pos.y, (real) pos.z, (real) pos.w);
+            posqCorrection[index] = make_real4(pos.x-(real) pos.x, pos.y-(real) pos.y, pos.z-(real) pos.z, 0);
+#else
+            posq[index] = pos;
+#endif
+        }
+    }
+}
+
+
+/**
+ * Perform the third step of Velocity Verlet integration.
+ *
+ * modify the velocities (velm) after the force update
+ */
+
+extern "C" __global__ void integrateVelocityVerletPart3(int numAtoms, int numPairs, int paddedNumAtoms, mixed2* __restrict__ dt, real4* __restrict__ posq,
+        real4* __restrict__ posqCorrection, mixed4* __restrict__ velm,  const long long* __restrict__ force, const mixed4* __restrict__ posDelta,
+        const int* __restrict__ atomList, const int2* __restrict__ pairList) {
+    mixed2 stepSize = dt[0];
+#if __CUDA_ARCH__ >= 130
+    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 = blockIdx.x*blockDim.x+threadIdx.x;
+    if (index == 0)
+        dt[0].x = stepSize.y;
+
+    for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < numAtoms; index += blockDim.x*gridDim.x) {
+        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;
+#if __CUDA_ARCH__ < 130
+            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;
+        }
+    }
+    for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < numPairs; index += blockDim.x*gridDim.x) {
+        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;
+        comFrc.x = force[atom1] + force[atom2];
+        comFrc.y = force[atom1 + paddedNumAtoms] + force[atom2 + paddedNumAtoms];
+        comFrc.z = force[atom1 + paddedNumAtoms*2] + force[atom2 + paddedNumAtoms*2];
+        mixed3 relFrc;
+        relFrc.x = mass1fract*force[atom2] - mass2fract*force[atom1];
+        relFrc.y = mass1fract*force[atom2+paddedNumAtoms] - mass2fract*force[atom1+paddedNumAtoms];
+        relFrc.z = mass1fract*force[atom2+paddedNumAtoms*2] - mass2fract*force[atom1+paddedNumAtoms*2];
+        comVel.x += comFrc.x * scale * invTotMass;
+        comVel.y += comFrc.y * scale * invTotMass;
+        comVel.z += comFrc.z * scale * invTotMass;
+        relVel.x += relFrc.x * scale * invRedMass;
+        relVel.y += relFrc.y * scale * invRedMass;
+        relVel.z += relFrc.z * scale * 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;
+#if __CUDA_ARCH__ < 130
+            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;
+#if __CUDA_ARCH__ < 130
+            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;
+        }
+    }
+}
+
+
+
+/**
+ * Apply the hard wall constraint 
+ */
+
+extern "C" __global__ void integrateVelocityVerletHardWall(int numPairs, const float* __restrict__ maxPairDistance, mixed2* __restrict__ dt, real4* __restrict__ posq,
+        real4* __restrict__ posqCorrection, mixed4* __restrict__ velm, 
+        const int2* __restrict__ pairList, const float* __restrict__ pairTemperature) {
+    mixed dtPos = dt[0].y;
+    mixed maxDelta = (mixed) maxPairDistance[0];
+    // Apply hard wall constraints.
+    if (maxDelta > 0) {
+        for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < numPairs; index += blockDim.x*gridDim.x) {
+            const mixed hardWallScale = sqrt( ((mixed) pairTemperature[index]) * ((mixed) BOLTZ));
+            int2 atom = make_int2(pairList[index].x, pairList[index].y);
+#ifdef USE_MIXED_PRECISION
+            real4 posv1 = posq[atom.x];
+            real4 posc1 = posqCorrection[atom.x];
+            mixed4 pos1 = make_mixed4(posv1.x+(mixed)posc1.x, posv1.y+(mixed)posc1.y, posv1.z+(mixed)posc1.z, posv1.w);
+            real4 posv2 = posq[atom.y];
+            real4 posc2 = posqCorrection[atom.y];
+            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[atom.x];
+            real4 pos2 = posq[atom.y];
+#endif
+            mixed3 delta = make_mixed3(
+                mixed (pos1.x - pos2.x),
+                mixed (pos1.y - pos2.y),
+                mixed (pos1.z - pos2.z)
+            );
+            mixed r = sqrt(delta.x*delta.x + delta.y*delta.y + delta.z*delta.z);
+            mixed rInv = 1/r;
+            if (rInv*maxDelta < 1.0) {
+                // The constraint has been violated, so make the inter-particle distance "bounce"
+                // off the hard wall.
+                mixed3 bondDir = make_mixed3(delta.x * rInv, delta.y * rInv, delta.z * rInv);
+                mixed3 vel1 = make_mixed3(velm[atom.x].x, velm[atom.x].y, velm[atom.x].z);
+                mixed3 vel2 = make_mixed3(velm[atom.y].x, velm[atom.y].y, velm[atom.y].z);
+                mixed m1 = velm[atom.x].w != 0.0 ? 1.0/velm[atom.x].w : 0.0;
+                mixed m2 = velm[atom.y].w != 0.0 ? 1.0/velm[atom.y].w : 0.0;
+                mixed invTotMass = (m1 + m2 != 0.0) ? 1.0 /(m1 + m2) : 0.0;
+                mixed deltaR = r-maxDelta;
+                mixed deltaT = dtPos;
+                mixed dt = dtPos;
+
+                mixed dotvr1 = vel1.x*bondDir.x + vel1.y*bondDir.y + vel1.z*bondDir.z;
+                mixed3 vb1 = make_mixed3(bondDir.x*dotvr1, bondDir.y*dotvr1, bondDir.z*dotvr1);
+                mixed3 vp1 = make_mixed3(vel1.x-vb1.x, vel1.y-vb1.y, vel1.z-vb1.z);
+                if (m2 == 0) {
+                    // The parent particle is massless, so move only the Drude particle.
+
+                    if (dotvr1 != 0.0)
+                        deltaT = deltaR/fabs(dotvr1);
+                    if (deltaT > dtPos)
+                        deltaT = dtPos;
+                    dotvr1 = -dotvr1*hardWallScale/(fabs(dotvr1)*sqrt(m1));
+                    mixed dr = -deltaR + deltaT*dotvr1;
+                    pos1.x += bondDir.x*dr;
+                    pos1.y += bondDir.y*dr;
+                    pos1.z += bondDir.z*dr;
+                    velm[atom.x] = make_mixed4(vp1.x + bondDir.x*dotvr1, vp1.y + bondDir.y*dotvr1, vp1.z + bondDir.z*dotvr1, velm[atom.x].w);
+#ifdef USE_MIXED_PRECISION
+                    posq[atom.x] = make_real4((real) pos1.x, (real) pos1.y, (real) pos1.z, (real) pos1.w);
+                    posqCorrection[atom.x] = make_real4(pos1.x-(real) pos1.x, pos1.y-(real) pos1.y, pos1.z-(real) pos1.z, 0);
+#else
+                    posq[atom.x] = pos1;
+#endif
+                }
+                else {
+                    // Move both particles.
+                    mixed dotvr2 = vel2.x*bondDir.x + vel2.y*bondDir.y + vel2.z*bondDir.z;
+                    mixed3 vb2 = make_mixed3(bondDir.x*dotvr2, bondDir.y*dotvr2, bondDir.z*dotvr2);
+                    mixed3 vp2 = make_mixed3(vel2.x-vb2.x, vel2.y-vb2.y, vel2.z-vb2.z);
+                    mixed vbCMass = (m1*dotvr1 + m2*dotvr2)*invTotMass;
+                    dotvr1 -= vbCMass;
+                    dotvr2 -= vbCMass;
+                    if (dotvr1 != dotvr2)
+                        deltaT = deltaR/fabs(dotvr1-dotvr2);
+                    if (deltaT > dt)
+                        deltaT = dt;
+                    mixed vBond = hardWallScale/sqrt(m1);
+                    dotvr1 = -dotvr1*vBond*m2*invTotMass/fabs(dotvr1);
+                    dotvr2 = -dotvr2*vBond*m1*invTotMass/fabs(dotvr2);
+                    mixed dr1 = -deltaR*m2*invTotMass + deltaT*dotvr1;
+                    mixed dr2 = deltaR*m1*invTotMass + deltaT*dotvr2;
+                    dotvr1 += vbCMass;
+                    dotvr2 += vbCMass;
+                    pos1.x += bondDir.x*dr1;
+                    pos1.y += bondDir.y*dr1;
+                    pos1.z += bondDir.z*dr1;
+                    pos2.x += bondDir.x*dr2;
+                    pos2.y += bondDir.y*dr2;
+                    pos2.z += bondDir.z*dr2;
+                    velm[atom.x] = make_mixed4(vp1.x + bondDir.x*dotvr1, vp1.y + bondDir.y*dotvr1, vp1.z + bondDir.z*dotvr1, velm[atom.x].w);
+                    velm[atom.y] = make_mixed4(vp2.x + bondDir.x*dotvr2, vp2.y + bondDir.y*dotvr2, vp2.z + bondDir.z*dotvr2, velm[atom.y].w);
+#ifdef USE_MIXED_PRECISION
+                    posq[atom.x] = make_real4((real) pos1.x, (real) pos1.y, (real) pos1.z, (real) pos1.w);
+                    posq[atom.y] = make_real4((real) pos2.x, (real) pos2.y, (real) pos2.z, (real) pos2.w);
+                    posqCorrection[atom.x] = make_real4(pos1.x-(real) pos1.x, pos1.y-(real) pos1.y, pos1.z-(real) pos1.z, 0);
+                    posqCorrection[atom.y] = make_real4(pos2.x-(real) pos2.x, pos2.y-(real) pos2.y, pos2.z-(real) pos2.z, 0);
+#else
+                    posq[atom.x] = pos1;
+                    posq[atom.y] = pos2;
+#endif
+                }
+            }
+        }
+    } /* end of hard wall constraint part */
+}

platforms/cuda/tests/TestCudaNoseHooverIntegrator.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. Simmmonett                           *
+ * 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 "TestNoseHooverIntegrator.h"
+
+void runPlatformTests() {
+}

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() {
+}