Merge https://github.com/openmm/openmm

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

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -33,7 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "openmm/AndersenThermostat.h"
-#include "openmm/BAOABLangevinIntegrator.h"
+#include "openmm/LangevinMiddleIntegrator.h"
 #include "openmm/BrownianIntegrator.h"
 #include "openmm/CMAPTorsionForce.h"
 #include "openmm/CMMotionRemover.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.
  */
@@ -1092,40 +1129,36 @@ public:
 };
 
 /**
- * This kernel is invoked by BAOABLangevinIntegrator to take one time step.
+ * This kernel is invoked by LangevinMiddleIntegrator to take one time step.
  */
-class IntegrateBAOABStepKernel : public KernelImpl {
+class IntegrateLangevinMiddleStepKernel : public KernelImpl {
 public:
     static std::string Name() {
-        return "IntegrateBAOABStep";
+        return "IntegrateLangevinMiddleStep";
     }
-    IntegrateBAOABStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    IntegrateLangevinMiddleStepKernel(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 BAOABLangevinIntegrator this kernel will be used for
+     * @param integrator the LangevinMiddleIntegrator this kernel will be used for
      */
-    virtual void initialize(const System& system, const BAOABLangevinIntegrator& integrator) = 0;
+    virtual void initialize(const System& system, const LangevinMiddleIntegrator& integrator) = 0;
     /**
      * Execute the kernel.
      * 
      * @param context    the context in which to execute this kernel
-     * @param integrator the BAOABLangevinIntegrator this kernel is being used for
-     * @param forcesAreValid if the context has been modified since the last time step, this will be
-     *                       false to show that cached forces are invalid and must be recalculated.
-     *                       On exit, this should specify whether the cached forces are valid at the
-     *                       end of the step.
+     * @param integrator the LangevinMiddleIntegrator this kernel is being used for
      */
-    virtual void execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid) = 0;
+    virtual void execute(ContextImpl& context, const LangevinMiddleIntegrator& integrator) = 0;
     /**
      * Compute the kinetic energy.
      * 
      * @param context    the context in which to execute this kernel
-     * @param integrator the BAOABLangevinIntegrator this kernel is being used for
+     * @param integrator the LangevinMiddleIntegrator this kernel is being used for
      */
-    virtual double computeKineticEnergy(ContextImpl& context, const BAOABLangevinIntegrator& integrator) = 0;
+    virtual double computeKineticEnergy(ContextImpl& context, const LangevinMiddleIntegrator& integrator) = 0;
 };
 
 /**
@@ -1327,6 +1360,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

@@ -33,7 +33,6 @@
  * -------------------------------------------------------------------------- */
 
 #include "openmm/AndersenThermostat.h"
-#include "openmm/BAOABLangevinIntegrator.h"
 #include "openmm/BrownianIntegrator.h"
 #include "openmm/CMAPTorsionForce.h"
 #include "openmm/CMMotionRemover.h"
@@ -57,6 +56,7 @@
 #include "openmm/HarmonicBondForce.h"
 #include "openmm/Integrator.h"
 #include "openmm/LangevinIntegrator.h"
+#include "openmm/LangevinMiddleIntegrator.h"
 #include "openmm/LocalEnergyMinimizer.h"
 #include "openmm/MonteCarloAnisotropicBarostat.h"
 #include "openmm/MonteCarloBarostat.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/CustomIntegrator.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2011-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -50,8 +50,8 @@ namespace OpenMM {
  * with the particle positions and momenta.
  *
  * To create an integration algorithm, you first define a set of variables the
- * integrator will compute.  Variables come in two types: <i>global</i> variables
- * have a single value, while <i>per-DOF</i> variables have a value for every
+ * integrator will compute.  Variables come in two types: global variables
+ * have a single value, while per-DOF variables have a value for every
  * degree of freedom (x, y, or z coordinate of a particle).  You can define as
  * many variables as you want of each type.  The value of any variable can be
  * computed by the integration algorithm, or set directly by calling a method on
@@ -112,14 +112,14 @@ namespace OpenMM {
  * evaluated, a different value will be used.  When used in a per-DOF
  * expression, a different value will be used for every degree of freedom.
  * Note, however, that if this variable appears multiple times in a single
- * expression, the <i>same</i> value is used everywhere it appears in that
+ * expression, the same value is used everywhere it appears in that
  * expression.</li>
  * <li>gaussian: (either global or per-DOF, read-only) This is a Gaussian
  * distributed random number with mean 0 and variance 1.  Every time an expression
  * is evaluated, a different value will be used.  When used in a per-DOF
  * expression, a different value will be used for every degree of freedom.
  * Note, however, that if this variable appears multiple times in a single
- * expression, the <i>same</i> value is used everywhere it appears in that
+ * expression, the same value is used everywhere it appears in that
  * expression.</li>
  * <li>A global variable is created for every adjustable parameter defined
  * in the integrator's Context.</li>
@@ -218,6 +218,43 @@ namespace OpenMM {
  * integrator.addComputePerDof("angularMomentum", "m*cross(x, v)");
  * </pre></tt>
  * 
+ * Here are two more examples that may be useful as starting points for writing
+ * your own integrators.  The first one implements the algorithm used by the
+ * standard VerletIntegrator class.  This is a leapfrog algorithm, in contrast
+ * to the velocity Verlet algorithm shown above, so it only requires applying
+ * constraints once in each time step.
+ * 
+ * <tt><pre>
+ * CustomIntegrator integrator(dt);
+ * integrator.addPerDofVariable("x0", 0);
+ * integrator.addUpdateContextState();
+ * integrator.addComputePerDof("x0", "x");
+ * integrator.addComputePerDof("v", "v+dt*f/m");
+ * integrator.addComputePerDof("x", "x+dt*v");
+ * integrator.addConstrainPositions();
+ * integrator.addComputePerDof("v", "(x-x0)/dt");
+ * </pre></tt>
+ * 
+ * The second one implements the algorithm used by the standard
+ * LangevinMiddleIntegrator class.  kB is Boltzmann's constant.
+ * 
+ * <tt><pre>
+ * CustomIntegrator integrator(dt);
+ * integrator.addGlobalVariable("a", exp(-friction*dt));
+ * integrator.addGlobalVariable("b", sqrt(1-exp(-2*friction*dt)));
+ * integrator.addGlobalVariable("kT", kB*temperature);
+ * integrator.addPerDofVariable("x1", 0);
+ * integrator.addUpdateContextState();
+ * integrator.addComputePerDof("v", "v + dt*f/m");
+ * integrator.addConstrainVelocities();
+ * integrator.addComputePerDof("x", "x + 0.5*dt*v");
+ * integrator.addComputePerDof("v", "a*v + b*sqrt(kT/m)*gaussian");
+ * integrator.addComputePerDof("x", "x + 0.5*dt*v");
+ * integrator.addComputePerDof("x1", "x");
+ * integrator.addConstrainPositions();
+ * integrator.addComputePerDof("v", "v + (x-x1)/dt");
+ * </pre></tt>
+ * 
  * Another feature of CustomIntegrator is that it can use derivatives of the
  * potential energy with respect to context parameters.  These derivatives are
  * typically computed by custom forces, and are only computed if a Force object
@@ -230,7 +267,7 @@ namespace OpenMM {
  *
  * An Integrator has one other job in addition to evolving the equations of motion:
  * it defines how to compute the kinetic energy of the system.  Depending on the
- * integration method used, simply summing mv<sup>2</sup>/2 over all degrees of
+ * integration method used, simply summing (mv^2)/2 over all degrees of
  * freedom may not give the correct answer.  For example, in a leapfrog integrator
  * the velocities are "delayed" by half a time step, so the above formula would
  * give the kinetic energy half a time step ago, not at the current time.
@@ -250,7 +287,7 @@ namespace OpenMM {
  * The kinetic energy expression may depend on the following pre-defined variables:
  * x, v, f, m, dt.  It also may depend on user-defined global and per-DOF variables,
  * and on the values of adjustable parameters defined  in the integrator's Context.
- * It may <i>not</i> depend on any other variable, such as the potential energy,
+ * It may not depend on any other variable, such as the potential energy,
  * the force from a single force group, or a random number.
  *
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following

openmmapi/include/openmm/BAOABLangevinIntegrator.h → openmmapi/include/openmm/LangevinMiddleIntegrator.h

-#ifndef OPENMM_BAOABLANGEVININTEGRATOR_H_
-#define OPENMM_BAOABLANGEVININTEGRATOR_H_
+#ifndef OPENMM_LANGEVINMIDDLEINTEGRATOR_H_
+#define OPENMM_LANGEVINMIDDLEINTEGRATOR_H_
 
 /* -------------------------------------------------------------------------- *
  *                                   OpenMM                                   *
@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -40,21 +40,21 @@ namespace OpenMM {
 
 /**
  * This is an Integrator which simulates a System using Langevin dynamics, with
- * the BAOAB discretization of Leimkuhler and Matthews (http://dx.doi.org/10.1093/amrx/abs010).
+ * the LFMiddle discretization (http://dx.doi.org/10.1021/acs.jpca.9b02771).
  * This method tend to produce more accurate configurational sampling than other
  * discretizations, such as the one used in LangevinIntegrator.
  */
 
-class OPENMM_EXPORT BAOABLangevinIntegrator : public Integrator {
+class OPENMM_EXPORT LangevinMiddleIntegrator : public Integrator {
 public:
     /**
-     * Create a BAOABLangevinIntegrator.
+     * Create a LangevinMiddleIntegrator.
      * 
      * @param temperature    the temperature of the heat bath (in Kelvin)
      * @param frictionCoeff  the friction coefficient which couples the system to the heat bath (in inverse picoseconds)
      * @param stepSize       the step size with which to integrate the system (in picoseconds)
      */
-    BAOABLangevinIntegrator(double temperature, double frictionCoeff, double stepSize);
+    LangevinMiddleIntegrator(double temperature, double frictionCoeff, double stepSize);
     /**
      * Get the temperature of the heat bath (in Kelvin).
      *
@@ -128,10 +128,6 @@ protected:
      * cleanup.  It will also get called again if the application calls reinitialize() on the Context.
      */
     void cleanup();
-    /**
-     * When the user modifies the state, we need to mark that the forces need to be recalculated.
-     */
-    void stateChanged(State::DataType changed);
     /**
      * Get the names of all Kernels used by this Integrator.
      */
@@ -147,10 +143,9 @@ protected:
 private:
     double temperature, friction;
     int randomNumberSeed;
-    bool forcesAreValid;
     Kernel kernel;
 };
 
 } // namespace OpenMM
 
-#endif /*OPENMM_BAOABLANGEVININTEGRATOR_H_*/
+#endif /*OPENMM_LANGEVINMIDDLEINTEGRATOR_H_*/

openmmapi/include/openmm/NonbondedForce.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2018 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -568,6 +568,32 @@ public:
                nonbondedMethod == NonbondedForce::PME ||
                nonbondedMethod == NonbondedForce::LJPME;
     }
+    /**
+     * Get whether periodic boundary conditions should be applied to exceptions.  Usually this is not
+     * appropriate, because exceptions are normally used to represent bonded interactions (1-2, 1-3, and
+     * 1-4 pairs), but there are situations when it does make sense.  For example, you may want to simulate
+     * an infinite chain where one end of a molecule is bonded to the opposite end of the next periodic
+     * copy.
+     * 
+     * Regardless of this value, periodic boundary conditions are only applied to exceptions if they also
+     * are applied to other interactions.  If the nonbonded method is NoCutoff or CutoffNonPeriodic, this
+     * value is ignored.  Also note that cutoffs are never applied to exceptions, again because they are
+     * normally used to represent bonded interactions.
+     */
+    bool getExceptionsUsePeriodicBoundaryConditions() const;
+    /**
+     * Set whether periodic boundary conditions should be applied to exceptions.  Usually this is not
+     * appropriate, because exceptions are normally used to represent bonded interactions (1-2, 1-3, and
+     * 1-4 pairs), but there are situations when it does make sense.  For example, you may want to simulate
+     * an infinite chain where one end of a molecule is bonded to the opposite end of the next periodic
+     * copy.
+     * 
+     * Regardless of this value, periodic boundary conditions are only applied to exceptions if they also
+     * get applied to other interactions.  If the nonbonded method is NoCutoff or CutoffNonPeriodic, this
+     * value is ignored.  Also note that cutoffs are never applied to exceptions, again because they are
+     * normally used to represent bonded interactions.
+     */
+    void setExceptionsUsePeriodicBoundaryConditions(bool periodic);
 protected:
     ForceImpl* createImpl() const;
 private:
@@ -578,7 +604,7 @@ private:
     class ExceptionOffsetInfo;
     NonbondedMethod nonbondedMethod;
     double cutoffDistance, switchingDistance, rfDielectric, ewaldErrorTol, alpha, dalpha;
-    bool useSwitchingFunction, useDispersionCorrection;
+    bool useSwitchingFunction, useDispersionCorrection, exceptionsUsePeriodic;
     int recipForceGroup, nx, ny, nz, dnx, dny, dnz;
     void addExclusionsToSet(const std::vector<std::set<int> >& bonded12, std::set<int>& exclusions, int baseParticle, int fromParticle, int currentLevel) const;
     int getGlobalParameterIndex(const std::string& parameter) const;

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/BAOABLangevinIntegrator.cpp → openmmapi/src/LangevinMiddleIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -29,7 +29,7 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-#include "openmm/BAOABLangevinIntegrator.h"
+#include "openmm/LangevinMiddleIntegrator.h"
 #include "openmm/Context.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
@@ -40,51 +40,47 @@ using namespace OpenMM;
 using std::string;
 using std::vector;
 
-BAOABLangevinIntegrator::BAOABLangevinIntegrator(double temperature, double frictionCoeff, double stepSize) {
+LangevinMiddleIntegrator::LangevinMiddleIntegrator(double temperature, double frictionCoeff, double stepSize) {
     setTemperature(temperature);
     setFriction(frictionCoeff);
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
     setRandomNumberSeed(0);
-    forcesAreValid = false;
 }
 
-void BAOABLangevinIntegrator::initialize(ContextImpl& contextRef) {
+void LangevinMiddleIntegrator::initialize(ContextImpl& contextRef) {
     if (owner != NULL && &contextRef.getOwner() != owner)
         throw OpenMMException("This Integrator is already bound to a context");
     context = &contextRef;
     owner = &contextRef.getOwner();
-    kernel = context->getPlatform().createKernel(IntegrateBAOABStepKernel::Name(), contextRef);
-    kernel.getAs<IntegrateBAOABStepKernel>().initialize(contextRef.getSystem(), *this);
+    kernel = context->getPlatform().createKernel(IntegrateLangevinMiddleStepKernel::Name(), contextRef);
+    kernel.getAs<IntegrateLangevinMiddleStepKernel>().initialize(contextRef.getSystem(), *this);
 }
 
-void BAOABLangevinIntegrator::cleanup() {
+void LangevinMiddleIntegrator::cleanup() {
     kernel = Kernel();
 }
 
-void BAOABLangevinIntegrator::stateChanged(State::DataType changed) {
-    forcesAreValid = false;
-}
-
-vector<string> BAOABLangevinIntegrator::getKernelNames() {
+vector<string> LangevinMiddleIntegrator::getKernelNames() {
     std::vector<std::string> names;
-    names.push_back(IntegrateBAOABStepKernel::Name());
+    names.push_back(IntegrateLangevinMiddleStepKernel::Name());
     return names;
 }
 
-double BAOABLangevinIntegrator::computeKineticEnergy() {
-    return kernel.getAs<IntegrateBAOABStepKernel>().computeKineticEnergy(*context, *this);
+double LangevinMiddleIntegrator::computeKineticEnergy() {
+    return kernel.getAs<IntegrateLangevinMiddleStepKernel>().computeKineticEnergy(*context, *this);
 }
 
-bool BAOABLangevinIntegrator::kineticEnergyRequiresForce() const {
+bool LangevinMiddleIntegrator::kineticEnergyRequiresForce() const {
     return false;
 }
 
-void BAOABLangevinIntegrator::step(int steps) {
+void LangevinMiddleIntegrator::step(int steps) {
     if (context == NULL)
         throw OpenMMException("This Integrator is not bound to a context!");  
     for (int i = 0; i < steps; ++i) {
         context->updateContextState();
-        kernel.getAs<IntegrateBAOABStepKernel>().execute(*context, *this, forcesAreValid);
+        context->calcForcesAndEnergy(true, false);
+        kernel.getAs<IntegrateLangevinMiddleStepKernel>().execute(*context, *this);
     }
 }

openmmapi/src/NonbondedForce.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2018 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -48,7 +48,7 @@ using std::stringstream;
 using std::vector;
 
 NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), switchingDistance(-1.0), rfDielectric(78.3),
-        ewaldErrorTol(5e-4), alpha(0.0), dalpha(0.0), useSwitchingFunction(false), useDispersionCorrection(true), recipForceGroup(-1),
+        ewaldErrorTol(5e-4), alpha(0.0), dalpha(0.0), useSwitchingFunction(false), useDispersionCorrection(true), exceptionsUsePeriodic(false), recipForceGroup(-1),
         nx(0), ny(0), nz(0), dnx(0), dny(0), dnz(0) {
 }
 
@@ -347,3 +347,11 @@ void NonbondedForce::setReciprocalSpaceForceGroup(int group) {
 void NonbondedForce::updateParametersInContext(Context& context) {
     dynamic_cast<NonbondedForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context));
 }
+
+bool NonbondedForce::getExceptionsUsePeriodicBoundaryConditions() const {
+    return exceptionsUsePeriodic;
+}
+
+void NonbondedForce::setExceptionsUsePeriodicBoundaryConditions(bool periodic) {
+    exceptionsUsePeriodic = periodic;
+}

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/common/CMakeLists.txt

+# Encode the kernel sources into a C++ class.
+
+SET(KERNEL_SOURCE_DIR "${CMAKE_CURRENT_SOURCE_DIR}/src")
+SET(KERNEL_SOURCE_CLASS CommonKernelSources)
+SET(KERNELS_CPP ${CMAKE_CURRENT_BINARY_DIR}/src/${KERNEL_SOURCE_CLASS}.cpp)
+SET(KERNELS_H ${CMAKE_CURRENT_BINARY_DIR}/src/${KERNEL_SOURCE_CLASS}.h)
+INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}/src)
+FILE(GLOB COMMON_KERNELS ${KERNEL_SOURCE_DIR}/kernels/*.cc)
+ADD_CUSTOM_COMMAND(OUTPUT ${KERNELS_CPP} ${KERNELS_H}
+    COMMAND ${CMAKE_COMMAND}
+    ARGS -D KERNEL_SOURCE_DIR=${KERNEL_SOURCE_DIR} -D KERNELS_CPP=${KERNELS_CPP} -D KERNELS_H=${KERNELS_H} -D KERNEL_SOURCE_CLASS=${KERNEL_SOURCE_CLASS} -D KERNEL_FILE_EXTENSION=cc -P ${CMAKE_SOURCE_DIR}/cmake_modules/EncodeKernelFiles.cmake
+    DEPENDS ${COMMON_KERNELS}
+)
+SET_SOURCE_FILES_PROPERTIES(${KERNELS_CPP} ${KERNELS_H} PROPERTIES GENERATED TRUE)
+ADD_CUSTOM_TARGET(CommonKernels DEPENDS ${KERNELS_CPP} ${KERNELS_H})
+
+# Install headers
+
+FILE(GLOB CORE_HEADERS include/openmm/common/*.h)
+INSTALL_FILES(/include/openmm/common FILES ${CORE_HEADERS})

platforms/common/include/openmm/common/ArrayInterface.h

+#ifndef OPENMM_ARRAYINTERFACE_H_
+#define OPENMM_ARRAYINTERFACE_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: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/OpenMMException.h"
+#include "openmm/common/windowsExportCommon.h"
+#include <vector>
+
+namespace OpenMM {
+
+class ComputeContext;
+
+/**
+ * This abstract class defines the interface for arrays stored on a computing device.
+ */
+
+class OPENMM_EXPORT_COMMON ArrayInterface {
+public:
+    virtual ~ArrayInterface() {
+    }
+    /**
+     * Initialize this array.
+     *
+     * @param context           the context for which to create the array
+     * @param size              the number of elements in the array
+     * @param elementSize       the size of each element in bytes
+     * @param name              the name of the array
+     */
+    virtual void initialize(ComputeContext& context, int size, int elementSize, const std::string& name) = 0;
+    /**
+     * Initialize this object.  The template argument is the data type of each array element.
+     *
+     * @param context           the context for which to create the array
+     * @param size              the number of elements in the array
+     * @param name              the name of the array
+     */
+    template <class T>
+    void initialize(ComputeContext& context, int size, const std::string& name) {
+        initialize(context, size, sizeof(T), name);
+    }
+    /**
+     * Recreate the internal storage to have a different size.
+     */
+    virtual void resize(int size) = 0;
+    /**
+     * Get whether this array has been initialized.
+     */
+    virtual bool isInitialized() const = 0;
+    /**
+     * Get the number of elements in the array.
+     */
+    virtual int getSize() const = 0;
+    /**
+     * Get the size of each element in bytes.
+     */
+    virtual int getElementSize() const = 0;
+    /**
+     * Get the name of the array.
+     */
+    virtual const std::string& getName() const = 0;
+    /**
+     * Get the context this array belongs to.
+     */
+    virtual ComputeContext& getContext() = 0;
+    /**
+     * Copy the values in a vector to the device memory.
+     * 
+     * @param data      the data in host memory to copy
+     * @param convert   if true, automatic conversions between single and double
+     *                  precision will be performed as necessary
+     */
+    template <class T>
+    void upload(const std::vector<T>& data, bool convert=false) {
+        if (convert && data.size() == getSize() && sizeof(T) != getElementSize()) {
+            if (sizeof(T) == 2*getElementSize()) {
+                // Convert values from double to single precision.
+                const double* d = reinterpret_cast<const double*>(&data[0]);
+                std::vector<float> v(getElementSize()*getSize()/sizeof(float));
+                for (int i = 0; i < v.size(); i++)
+                    v[i] = (float) d[i];
+                upload(&v[0], true);
+                return;
+            }
+            if (2*sizeof(T) == getElementSize()) {
+                // Convert values from single to double precision.
+                const float* d = reinterpret_cast<const float*>(&data[0]);
+                std::vector<double> v(getElementSize()*getSize()/sizeof(double));
+                for (int i = 0; i < v.size(); i++)
+                    v[i] = (double) d[i];
+                upload(&v[0], true);
+                return;
+            }
+        }
+        if (sizeof(T) != getElementSize() || data.size() != getSize())
+            throw OpenMMException("Error uploading array "+getName()+": The specified vector does not match the size of the array");
+        upload(&data[0], true);
+    }
+    /**
+     * Copy the values in the array to a vector.
+     */
+    template <class T>
+    void download(std::vector<T>& data) const {
+        if (sizeof(T) != getElementSize())
+            throw OpenMMException("Error downloading array "+getName()+": The specified vector has the wrong element size");
+        if (data.size() != getSize())
+            data.resize(getSize());
+        download(&data[0], true);
+    }
+    /**
+     * Copy the values from host memory to the array.
+     * 
+     * @param data     the data to copy
+     * @param blocking if true, this call will block until the transfer is complete.  Subclasses often
+     *                 have restrictions on non-blocking copies, such as that the source data must be
+     *                 in page-locked memory.
+     */
+    virtual void upload(const void* data, bool blocking=true) = 0;
+    /**
+     * Copy the values in the array to host memory.
+     * 
+     * @param data     the destination to copy the value to
+     * @param blocking if true, this call will block until the transfer is complete.  Subclasses often
+     *                 have restrictions on non-blocking copies, such as that the destination must be
+     *                 in page-locked memory.
+     */
+    virtual void download(void* data, bool blocking=true) const = 0;
+    /**
+     * Copy the values in this array to a second array.
+     * 
+     * @param dest     the destination array to copy to
+     */
+    virtual void copyTo(ArrayInterface& dest) const = 0;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_ARRAYINTERFACE_H_*/

platforms/common/include/openmm/common/BondedUtilities.h

+#ifndef OPENMM_BONDEDUTILITIES_H_
+#define OPENMM_BONDEDUTILITIES_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) 2011-2019 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/common/ArrayInterface.h"
+#include <string>
+#include <vector>
+
+namespace OpenMM {
+
+/**
+ * This abstract class defines an interface for computing bonded interactions.  Call
+ * getBondedUtilities() on a ComputeContext to get the BondedUtilities object for that
+ * context.
+ * 
+ * This class provides a generic mechanism for evaluating bonded interactions.  You write only
+ * the source code needed to compute one interaction, and this object takes care of creating
+ * and executing a complete kernel that loops over bonds, evaluates each one, and accumulates
+ * the resulting forces and energies.  This offers two advantages.  First, it simplifies the
+ * task of writing a new Force.  Second, it allows multiple forces to be evaluated by a single
+ * kernel, which reduces overhead and improves performance.
+ * 
+ * A "bonded interaction" means an interaction that affects a small, fixed set of particles.
+ * The interaction energy may depend on the positions of only those particles, and the list of
+ * particles forming a "bond" may not change with time.  Examples of bonded interactions
+ * include HarmonicBondForce, HarmonicAngleForce, and PeriodicTorsionForce.
+ * 
+ * To create a bonded interaction, call addInteraction().  You pass to it a block of source
+ * code for evaluating the interaction.  The inputs and outputs for that source code are as
+ * follows:
+ * 
+ * <ol>
+ * <li>The index of the bond being evaluated will have been stored in the unsigned int variable "index".</li>
+ * <li>The indices of the atoms forming that bond will have been stored in the unsigned int variables "atom1",
+ * "atom2", ....</li>
+ * <li>The positions of those atoms will have been stored in the real4 variables "pos1", "pos2", ....</li>
+ * <li>A real variable called "energy" will exist.  Your code should add the potential energy of the
+ * bond to that variable.</li>
+ * <li>Your code should define real3 variables called "force1", "force2", ... that contain the force to
+ * apply to each atom.</li>
+ * </ol>
+ * 
+ * As a simple example, the following source code would be used to implement a pairwise interaction of
+ * the form E=r^2:
+ * 
+ * <tt><pre>
+ * real4 delta = pos2-pos1;
+ * energy += delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
+ * real3 force1 = 2.0f*delta;
+ * real3 force2 = -2.0f*delta;
+ * </pre></tt>
+ * 
+ * Interactions will often depend on parameters or other data.  Call addArgument() to provide the data
+ * to this class.  It will be passed to the interaction kernel as an argument, and you can refer to it
+ * from your interaction code.
+ */
+
+class OPENMM_EXPORT_COMMON BondedUtilities {
+public:
+    virtual ~BondedUtilities() {
+    }
+    /**
+     * Add a bonded interaction.
+     *
+     * @param atoms    this should have one entry for each bond, and that entry should contain the list
+     *                 of atoms involved in the bond.  Every entry must have the same number of atoms.
+     * @param source   the code to evaluate the interaction
+     * @param group    the force group in which the interaction should be calculated
+     */
+    virtual void addInteraction(const std::vector<std::vector<int> >& atoms, const std::string& source, int group) = 0;
+    /**
+     * Add an argument that should be passed to the interaction kernel.
+     * 
+     * @param data    the array containing the data to pass
+     * @param type    the data type contained in the memory (e.g. "float4")
+     * @return the name that will be used for the argument.  Any code you pass to addInteraction() should
+     * refer to it by this name.
+     */
+    virtual std::string addArgument(ArrayInterface& data, const std::string& type) = 0;
+    /**
+     * Register that the interaction kernel will be computing the derivative of the potential energy
+     * with respect to a parameter.
+     * 
+     * @param param   the name of the parameter
+     * @return the variable that will be used to accumulate the derivative.  Any code you pass to addInteraction() should
+     * add its contributions to this variable.
+     */
+    virtual std::string addEnergyParameterDerivative(const std::string& param) = 0;
+    /**
+     * Add some code that should be included in the program, before the start of the kernel.
+     * This can be used, for example, to define functions that will be called by the kernel.
+     * 
+     * @param source   the code to include
+     */
+    virtual void addPrefixCode(const std::string& source) = 0;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_BONDEDUTILITIES_H_*/

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

+#ifndef OPENMM_COMMONKERNELS_H_
+#define OPENMM_COMMONKERNELS_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) 2008-2020 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/common/ComputeArray.h"
+#include "openmm/common/ComputeContext.h"
+#include "openmm/common/ComputeParameterSet.h"
+#include "openmm/Platform.h"
+#include "openmm/kernels.h"
+#include "openmm/internal/CompiledExpressionSet.h"
+#include "openmm/internal/CustomIntegratorUtilities.h"
+#include "lepton/CompiledExpression.h"
+
+namespace OpenMM {
+
+
+/**
+ * This kernel is invoked by HarmonicBondForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CommonCalcHarmonicBondForceKernel : public CalcHarmonicBondForceKernel {
+public:
+    CommonCalcHarmonicBondForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcHarmonicBondForceKernel(name, platform),
+            hasInitializedKernel(false), cc(cc), system(system) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the HarmonicBondForce this kernel will be used for
+     */
+    void initialize(const System& system, const HarmonicBondForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the HarmonicBondForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const HarmonicBondForce& force);
+private:
+    class ForceInfo;
+    int numBonds;
+    bool hasInitializedKernel;
+    ComputeContext& cc;
+    ForceInfo* info;
+    const System& system;
+    ComputeArray params;
+};
+
+/**
+ * This kernel is invoked by CustomBondForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CommonCalcCustomBondForceKernel : public CalcCustomBondForceKernel {
+public:
+    CommonCalcCustomBondForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcCustomBondForceKernel(name, platform),
+            hasInitializedKernel(false), cc(cc), system(system), params(NULL) {
+    }
+    ~CommonCalcCustomBondForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomBondForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomBondForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomBondForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomBondForce& force);
+private:
+    class ForceInfo;
+    int numBonds;
+    bool hasInitializedKernel;
+    ComputeContext& cc;
+    ForceInfo* info;
+    const System& system;
+    ComputeParameterSet* params;
+    ComputeArray globals;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+};
+
+/**
+ * This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CommonCalcHarmonicAngleForceKernel : public CalcHarmonicAngleForceKernel {
+public:
+    CommonCalcHarmonicAngleForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcHarmonicAngleForceKernel(name, platform),
+            hasInitializedKernel(false), cc(cc), system(system) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the HarmonicAngleForce this kernel will be used for
+     */
+    void initialize(const System& system, const HarmonicAngleForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the HarmonicAngleForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& force);
+private:
+    class ForceInfo;
+    int numAngles;
+    bool hasInitializedKernel;
+    ComputeContext& cc;
+    ForceInfo* info;
+    const System& system;
+    ComputeArray params;
+};
+
+/**
+ * This kernel is invoked by CustomAngleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CommonCalcCustomAngleForceKernel : public CalcCustomAngleForceKernel {
+public:
+    CommonCalcCustomAngleForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcCustomAngleForceKernel(name, platform),
+            hasInitializedKernel(false), cc(cc), system(system), params(NULL) {
+    }
+    ~CommonCalcCustomAngleForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomAngleForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomAngleForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomAngleForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomAngleForce& force);
+private:
+    class ForceInfo;
+    int numAngles;
+    bool hasInitializedKernel;
+    ComputeContext& cc;
+    ForceInfo* info;
+    const System& system;
+    ComputeParameterSet* params;
+    ComputeArray globals;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+};
+
+/**
+ * This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CommonCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
+public:
+    CommonCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcPeriodicTorsionForceKernel(name, platform),
+            hasInitializedKernel(false), cc(cc), system(system) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the PeriodicTorsionForce this kernel will be used for
+     */
+    void initialize(const System& system, const PeriodicTorsionForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the PeriodicTorsionForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force);
+private:
+    class ForceInfo;
+    int numTorsions;
+    bool hasInitializedKernel;
+    ComputeContext& cc;
+    ForceInfo* info;
+    const System& system;
+    ComputeArray params;
+};
+
+/**
+ * This kernel is invoked by RBTorsionForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CommonCalcRBTorsionForceKernel : public CalcRBTorsionForceKernel {
+public:
+    CommonCalcRBTorsionForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcRBTorsionForceKernel(name, platform),
+            hasInitializedKernel(false), cc(cc), system(system) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the RBTorsionForce this kernel will be used for
+     */
+    void initialize(const System& system, const RBTorsionForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the RBTorsionForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const RBTorsionForce& force);
+private:
+    class ForceInfo;
+    int numTorsions;
+    bool hasInitializedKernel;
+    ComputeContext& cc;
+    ForceInfo* info;
+    const System& system;
+    ComputeArray params1;
+    ComputeArray params2;
+};
+
+/**
+ * This kernel is invoked by CustomTorsionForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CommonCalcCustomTorsionForceKernel : public CalcCustomTorsionForceKernel {
+public:
+    CommonCalcCustomTorsionForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcCustomTorsionForceKernel(name, platform),
+            hasInitializedKernel(false), cc(cc), system(system), params(NULL) {
+    }
+    ~CommonCalcCustomTorsionForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomTorsionForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomTorsionForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomTorsionForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomTorsionForce& force);
+private:
+    class ForceInfo;
+    int numTorsions;
+    bool hasInitializedKernel;
+    ComputeContext& cc;
+    ForceInfo* info;
+    const System& system;
+    ComputeParameterSet* params;
+    ComputeArray globals;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+};
+
+/**
+ * This kernel is invoked by CMAPTorsionForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CommonCalcCMAPTorsionForceKernel : public CalcCMAPTorsionForceKernel {
+public:
+    CommonCalcCMAPTorsionForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcCMAPTorsionForceKernel(name, platform),
+            hasInitializedKernel(false), cc(cc), system(system) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CMAPTorsionForce this kernel will be used for
+     */
+    void initialize(const System& system, const CMAPTorsionForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CMAPTorsionForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force);
+private:
+    class ForceInfo;
+    int numTorsions;
+    bool hasInitializedKernel;
+    ComputeContext& cc;
+    ForceInfo* info;
+    const System& system;
+    std::vector<mm_int2> mapPositionsVec;
+    ComputeArray coefficients;
+    ComputeArray mapPositions;
+    ComputeArray torsionMaps;
+};
+
+/**
+ * This kernel is invoked by CustomExternalForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CommonCalcCustomExternalForceKernel : public CalcCustomExternalForceKernel {
+public:
+    CommonCalcCustomExternalForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcCustomExternalForceKernel(name, platform),
+            hasInitializedKernel(false), cc(cc), system(system), params(NULL) {
+    }
+    ~CommonCalcCustomExternalForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomExternalForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomExternalForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomExternalForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomExternalForce& force);
+private:
+    class ForceInfo;
+    int numParticles;
+    bool hasInitializedKernel;
+    ComputeContext& cc;
+    ForceInfo* info;
+    const System& system;
+    ComputeParameterSet* params;
+    ComputeArray globals;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+};
+
+/**
+ * This kernel is invoked by CustomCompoundBondForce to calculate the forces acting on the system.
+ */
+class CommonCalcCustomCompoundBondForceKernel : public CalcCustomCompoundBondForceKernel {
+public:
+    CommonCalcCustomCompoundBondForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcCustomCompoundBondForceKernel(name, platform),
+            cc(cc), params(NULL), system(system) {
+    }
+    ~CommonCalcCustomCompoundBondForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomCompoundBondForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomCompoundBondForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomCompoundBondForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force);
+
+private:
+    class ForceInfo;
+    int numBonds;
+    ComputeContext& cc;
+    ForceInfo* info;
+    ComputeParameterSet* params;
+    ComputeArray globals;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+    std::vector<ComputeArray> tabulatedFunctions;
+    const System& system;
+};
+
+/**
+ * This kernel is invoked by CustomCentroidBondForce to calculate the forces acting on the system.
+ */
+class CommonCalcCustomCentroidBondForceKernel : public CalcCustomCentroidBondForceKernel {
+public:
+    CommonCalcCustomCentroidBondForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcCustomCentroidBondForceKernel(name, platform),
+            cc(cc), params(NULL), system(system) {
+    }
+    ~CommonCalcCustomCentroidBondForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomCentroidBondForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomCentroidBondForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomCentroidBondForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force);
+
+private:
+    class ForceInfo;
+    int numGroups, numBonds;
+    bool needEnergyParamDerivs;
+    ComputeContext& cc;
+    ForceInfo* info;
+    ComputeParameterSet* params;
+    ComputeArray globals, groupParticles, groupWeights, groupOffsets;
+    ComputeArray groupForces, bondGroups, centerPositions;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+    std::vector<ComputeArray> tabulatedFunctions;
+    std::vector<void*> groupForcesArgs;
+    ComputeKernel computeCentersKernel, groupForcesKernel, applyForcesKernel;
+    const System& system;
+};
+
+/**
+ * This kernel is invoked by CustomNonbondedForce to calculate the forces acting on the system.
+ */
+class CommonCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
+public:
+    CommonCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcCustomNonbondedForceKernel(name, platform),
+            cc(cc), params(NULL), forceCopy(NULL), system(system), hasInitializedKernel(false) {
+    }
+    ~CommonCalcCustomNonbondedForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomNonbondedForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomNonbondedForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomNonbondedForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force);
+private:
+    class ForceInfo;
+    void initInteractionGroups(const CustomNonbondedForce& force, const std::string& interactionSource, const std::vector<std::string>& tableTypes);
+    ComputeContext& cc;
+    ForceInfo* info;
+    ComputeParameterSet* params;
+    ComputeArray globals, interactionGroupData, filteredGroupData, numGroupTiles;
+    ComputeKernel interactionGroupKernel, prepareNeighborListKernel, buildNeighborListKernel;
+    std::vector<void*> interactionGroupArgs;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+    std::vector<ComputeArray> tabulatedFunctions;
+    double longRangeCoefficient;
+    std::vector<double> longRangeCoefficientDerivs;
+    bool hasInitializedLongRangeCorrection, hasInitializedKernel, hasParamDerivs, useNeighborList;
+    int numGroupThreadBlocks;
+    CustomNonbondedForce* forceCopy;
+    const System& system;
+};
+
+/**
+ * This kernel is invoked by GBSAOBCForce to calculate the forces acting on the system.
+ */
+class CommonCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
+public:
+    CommonCalcGBSAOBCForceKernel(std::string name, const Platform& platform, ComputeContext& cc) : CalcGBSAOBCForceKernel(name, platform), cc(cc),
+            hasCreatedKernels(false) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the GBSAOBCForce this kernel will be used for
+     */
+    void initialize(const System& system, const GBSAOBCForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the GBSAOBCForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force);
+private:
+    class ForceInfo;
+    double prefactor, surfaceAreaFactor, cutoff;
+    bool hasCreatedKernels;
+    int maxTiles;
+    ComputeContext& cc;
+    ForceInfo* info;
+    ComputeArray params, charges, bornSum, bornRadii, bornForce, obcChain;
+    ComputeKernel computeBornSumKernel, reduceBornSumKernel, force1Kernel, reduceBornForceKernel;
+};
+
+/**
+ * This kernel is invoked by CustomGBForce to calculate the forces acting on the system.
+ */
+class CommonCalcCustomGBForceKernel : public CalcCustomGBForceKernel {
+public:
+    CommonCalcCustomGBForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcCustomGBForceKernel(name, platform),
+            hasInitializedKernels(false), cc(cc), params(NULL), computedValues(NULL), energyDerivs(NULL), energyDerivChain(NULL), system(system) {
+    }
+    ~CommonCalcCustomGBForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomGBForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomGBForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomGBForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomGBForce& force);
+private:
+    class ForceInfo;
+    double cutoff;
+    bool hasInitializedKernels, needParameterGradient, needEnergyParamDerivs;
+    int maxTiles, numComputedValues;
+    ComputeContext& cc;
+    ForceInfo* info;
+    ComputeParameterSet* params;
+    ComputeParameterSet* computedValues;
+    ComputeParameterSet* energyDerivs;
+    ComputeParameterSet* energyDerivChain;
+    std::vector<ComputeParameterSet*> dValuedParam;
+    std::vector<ComputeArray> dValue0dParam;
+    ComputeArray longEnergyDerivs, globals, valueBuffers;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+    std::vector<ComputeArray> tabulatedFunctions;
+    std::vector<bool> pairValueUsesParam, pairEnergyUsesParam, pairEnergyUsesValue;
+    const System& system;
+    ComputeKernel pairValueKernel, perParticleValueKernel, pairEnergyKernel, perParticleEnergyKernel, gradientChainRuleKernel;
+    std::string pairValueSrc, pairEnergySrc;
+    std::map<std::string, std::string> pairValueDefines, pairEnergyDefines;
+};
+
+/**
+ * This kernel is invoked by CustomHbondForce to calculate the forces acting on the system.
+ */
+class CommonCalcCustomHbondForceKernel : public CalcCustomHbondForceKernel {
+public:
+    CommonCalcCustomHbondForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcCustomHbondForceKernel(name, platform),
+            hasInitializedKernel(false), cc(cc), donorParams(NULL), acceptorParams(NULL), system(system) {
+    }
+    ~CommonCalcCustomHbondForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomHbondForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomHbondForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomHbondForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomHbondForce& force);
+private:
+    class ForceInfo;
+    int numDonors, numAcceptors;
+    bool hasInitializedKernel;
+    ComputeContext& cc;
+    ForceInfo* info;
+    ComputeParameterSet* donorParams;
+    ComputeParameterSet* acceptorParams;
+    ComputeArray globals;
+    ComputeArray donors;
+    ComputeArray acceptors;
+    ComputeArray donorBufferIndices;
+    ComputeArray acceptorBufferIndices;
+    ComputeArray donorExclusions;
+    ComputeArray acceptorExclusions;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+    std::vector<ComputeArray> tabulatedFunctions;
+    const System& system;
+    ComputeKernel donorKernel, acceptorKernel;
+};
+
+/**
+ * This kernel is invoked by CustomManyParticleForce to calculate the forces acting on the system.
+ */
+class CommonCalcCustomManyParticleForceKernel : public CalcCustomManyParticleForceKernel {
+public:
+    CommonCalcCustomManyParticleForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcCustomManyParticleForceKernel(name, platform),
+            hasInitializedKernel(false), cc(cc), params(NULL), system(system) {
+    }
+    ~CommonCalcCustomManyParticleForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomManyParticleForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomManyParticleForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomManyParticleForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force);
+
+private:
+    class ForceInfo;
+    ComputeContext& cc;
+    ForceInfo* info;
+    bool hasInitializedKernel;
+    NonbondedMethod nonbondedMethod;
+    int maxNeighborPairs, forceWorkgroupSize, findNeighborsWorkgroupSize;
+    ComputeParameterSet* params;
+    ComputeArray globals, particleTypes,  orderIndex, particleOrder;
+    ComputeArray exclusions, exclusionStartIndex, blockCenter, blockBoundingBox;
+    ComputeArray neighborPairs, numNeighborPairs, neighborStartIndex, numNeighborsForAtom, neighbors;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+    std::vector<ComputeArray> tabulatedFunctions;
+    const System& system;
+    ComputeKernel forceKernel, blockBoundsKernel, neighborsKernel, startIndicesKernel, copyPairsKernel;
+    ComputeEvent event;
+};
+
+/**
+ * This kernel is invoked by GayBerneForce to calculate the forces acting on the system.
+ */
+class CommonCalcGayBerneForceKernel : public CalcGayBerneForceKernel {
+public:
+    CommonCalcGayBerneForceKernel(std::string name, const Platform& platform, ComputeContext& cc) : CalcGayBerneForceKernel(name, platform), cc(cc),
+            hasInitializedKernels(false) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the GayBerneForce this kernel will be used for
+     */
+    void initialize(const System& system, const GayBerneForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the GayBerneForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const GayBerneForce& force);
+private:
+    class ForceInfo;
+    class ReorderListener;
+    void sortAtoms();
+    ComputeContext& cc;
+    ForceInfo* info;
+    bool hasInitializedKernels;
+    int numRealParticles, maxNeighborBlocks;
+    GayBerneForce::NonbondedMethod nonbondedMethod;
+    ComputeArray sortedParticles, axisParticleIndices, sigParams, epsParams;
+    ComputeArray scale, exceptionParticles, exceptionParams;
+    ComputeArray aMatrix, bMatrix, gMatrix;
+    ComputeArray exclusions, exclusionStartIndex, blockCenter, blockBoundingBox;
+    ComputeArray neighbors, neighborIndex, neighborBlockCount;
+    ComputeArray sortedPos, torque;
+    std::vector<bool> isRealParticle;
+    std::vector<std::pair<int, int> > exceptionAtoms;
+    std::vector<std::pair<int, int> > excludedPairs;
+    ComputeKernel framesKernel, blockBoundsKernel, neighborsKernel, forceKernel, torqueKernel;
+    ComputeEvent event;
+};
+
+/**
+ * This kernel is invoked by VerletIntegrator to take one time step.
+ */
+class CommonIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
+public:
+    CommonIntegrateVerletStepKernel(std::string name, const Platform& platform, ComputeContext& cc) : IntegrateVerletStepKernel(name, platform), cc(cc),
+            hasInitializedKernels(false) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param integrator the VerletIntegrator this kernel will be used for
+     */
+    void initialize(const System& system, const VerletIntegrator& integrator);
+    /**
+     * Execute the kernel.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param integrator the VerletIntegrator this kernel is being used for
+     */
+    void execute(ContextImpl& context, const VerletIntegrator& integrator);
+    /**
+     * Compute the kinetic energy.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the VerletIntegrator this kernel is being used for
+     */
+    double computeKineticEnergy(ContextImpl& context, const VerletIntegrator& integrator);
+private:
+    ComputeContext& cc;
+    bool hasInitializedKernels;
+    ComputeKernel kernel1, kernel2;
+};
+
+/**
+ * This kernel is invoked by LangevinIntegrator to take one time step.
+ */
+class CommonIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
+public:
+    CommonIntegrateLangevinStepKernel(std::string name, const Platform& platform, ComputeContext& cc) : IntegrateLangevinStepKernel(name, platform), cc(cc),
+            hasInitializedKernels(false) {
+    }
+    /**
+     * Initialize the kernel, setting up the particle masses.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param integrator the LangevinIntegrator this kernel will be used for
+     */
+    void initialize(const System& system, const LangevinIntegrator& integrator);
+    /**
+     * Execute the kernel.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param integrator the LangevinIntegrator this kernel is being used for
+     */
+    void execute(ContextImpl& context, const LangevinIntegrator& integrator);
+    /**
+     * Compute the kinetic energy.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the LangevinIntegrator this kernel is being used for
+     */
+    double computeKineticEnergy(ContextImpl& context, const LangevinIntegrator& integrator);
+private:
+    ComputeContext& cc;
+    double prevTemp, prevFriction, prevStepSize;
+    bool hasInitializedKernels;
+    ComputeArray params;
+    ComputeKernel kernel1, kernel2;
+};
+
+/**
+ * This kernel is invoked by LangevinMiddleIntegrator to take one time step.
+ */
+class CommonIntegrateLangevinMiddleStepKernel : public IntegrateLangevinMiddleStepKernel {
+public:
+    CommonIntegrateLangevinMiddleStepKernel(std::string name, const Platform& platform, ComputeContext& cc) : IntegrateLangevinMiddleStepKernel(name, platform), cc(cc),
+            hasInitializedKernels(false) {
+    }
+    /**
+     * Initialize the kernel, setting up the particle masses.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param integrator the LangevinMiddleIntegrator this kernel will be used for
+     */
+    void initialize(const System& system, const LangevinMiddleIntegrator& integrator);
+    /**
+     * Execute the kernel.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the LangevinMiddleIntegrator this kernel is being used for
+     */
+    void execute(ContextImpl& context, const LangevinMiddleIntegrator& integrator);
+    /**
+     * Compute the kinetic energy.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the LangevinMiddleIntegrator this kernel is being used for
+     */
+    double computeKineticEnergy(ContextImpl& context, const LangevinMiddleIntegrator& integrator);
+private:
+    ComputeContext& cc;
+    double prevTemp, prevFriction, prevStepSize;
+    bool hasInitializedKernels;
+    ComputeArray params, oldDelta;
+    ComputeKernel kernel1, kernel2, kernel3;
+};
+
+/**
+ * This kernel is invoked by BrownianIntegrator to take one time step.
+ */
+class CommonIntegrateBrownianStepKernel : public IntegrateBrownianStepKernel {
+public:
+    CommonIntegrateBrownianStepKernel(std::string name, const Platform& platform, ComputeContext& cc) : IntegrateBrownianStepKernel(name, platform), cc(cc),
+            hasInitializedKernels(false), prevTemp(-1), prevFriction(-1), prevStepSize(-1) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param integrator the BrownianIntegrator this kernel will be used for
+     */
+    void initialize(const System& system, const BrownianIntegrator& integrator);
+    /**
+     * Execute the kernel.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param integrator the BrownianIntegrator this kernel is being used for
+     */
+    void execute(ContextImpl& context, const BrownianIntegrator& integrator);
+    /**
+     * Compute the kinetic energy.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the BrownianIntegrator this kernel is being used for
+     */
+    double computeKineticEnergy(ContextImpl& context, const BrownianIntegrator& integrator);
+private:
+    ComputeContext& cc;
+    double prevTemp, prevFriction, prevStepSize;
+    bool hasInitializedKernels;
+    ComputeKernel kernel1, kernel2;
+};
+
+/**
+ * This kernel is invoked by VerletIntegrator to take one time step.
+ */
+class CommonIntegrateVariableVerletStepKernel : public IntegrateVariableVerletStepKernel {
+public:
+    CommonIntegrateVariableVerletStepKernel(std::string name, const Platform& platform, ComputeContext& cc) : IntegrateVariableVerletStepKernel(name, platform), cc(cc),
+            hasInitializedKernels(false) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param integrator the VariableVerletIntegrator this kernel will be used for
+     */
+    void initialize(const System& system, const VariableVerletIntegrator& integrator);
+    /**
+     * Execute the kernel.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param integrator the VariableVerletIntegrator this kernel is being used for
+     * @param maxTime    the maximum time beyond which the simulation should not be advanced
+     * @return the size of the step that was taken
+     */
+    double execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
+    /**
+     * Compute the kinetic energy.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the VariableVerletIntegrator this kernel is being used for
+     */
+    double computeKineticEnergy(ContextImpl& context, const VariableVerletIntegrator& integrator);
+private:
+    ComputeContext& cc;
+    bool hasInitializedKernels;
+    int blockSize;
+    ComputeKernel kernel1, kernel2, selectSizeKernel;
+};
+
+/**
+ * This kernel is invoked by VariableLangevinIntegrator to take one time step.
+ */
+class CommonIntegrateVariableLangevinStepKernel : public IntegrateVariableLangevinStepKernel {
+public:
+    CommonIntegrateVariableLangevinStepKernel(std::string name, const Platform& platform, ComputeContext& cc) : IntegrateVariableLangevinStepKernel(name, platform), cc(cc),
+            hasInitializedKernels(false) {
+    }
+    /**
+     * Initialize the kernel, setting up the particle masses.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param integrator the VariableLangevinIntegrator this kernel will be used for
+     */
+    void initialize(const System& system, const VariableLangevinIntegrator& integrator);
+    /**
+     * Execute the kernel.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param integrator the VariableLangevinIntegrator this kernel is being used for
+     * @param maxTime    the maximum time beyond which the simulation should not be advanced
+     * @return the size of the step that was taken
+     */
+    double execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
+    /**
+     * Compute the kinetic energy.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the VariableLangevinIntegrator this kernel is being used for
+     */
+    double computeKineticEnergy(ContextImpl& context, const VariableLangevinIntegrator& integrator);
+private:
+    ComputeContext& cc;
+    bool hasInitializedKernels;
+    int blockSize;
+    ComputeArray params;
+    ComputeKernel kernel1, kernel2, selectSizeKernel;
+    double prevTemp, prevFriction, prevErrorTol;
+};
+
+/**
+ * This kernel is invoked by CustomIntegrator to take one time step.
+ */
+class CommonIntegrateCustomStepKernel : public IntegrateCustomStepKernel {
+public:
+    enum GlobalTargetType {DT, VARIABLE, PARAMETER};
+    CommonIntegrateCustomStepKernel(std::string name, const Platform& platform, ComputeContext& cc) : IntegrateCustomStepKernel(name, platform), cc(cc),
+            hasInitializedKernels(false), needsEnergyParamDerivs(false) {
+    }
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param integrator the CustomIntegrator this kernel will be used for
+     */
+    void initialize(const System& system, const CustomIntegrator& integrator);
+    /**
+     * Execute the kernel.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the CustomIntegrator this kernel is being used for
+     * @param forcesAreValid if the context has been modified since the last time step, this will be
+     *                       false to show that cached forces are invalid and must be recalculated.
+     *                       On exit, this should specify whether the cached forces are valid at the
+     *                       end of the step.
+     */
+    void execute(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid);
+    /**
+     * Compute the kinetic energy.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the CustomIntegrator this kernel is being used for
+     * @param forcesAreValid if the context has been modified since the last time step, this will be
+     *                       false to show that cached forces are invalid and must be recalculated.
+     *                       On exit, this should specify whether the cached forces are valid at the
+     *                       end of the step.
+     */
+    double computeKineticEnergy(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid);
+    /**
+     * Get the values of all global variables.
+     *
+     * @param context   the context in which to execute this kernel
+     * @param values    on exit, this contains the values
+     */
+    void getGlobalVariables(ContextImpl& context, std::vector<double>& values) const;
+    /**
+     * Set the values of all global variables.
+     *
+     * @param context   the context in which to execute this kernel
+     * @param values    a vector containing the values
+     */
+    void setGlobalVariables(ContextImpl& context, const std::vector<double>& values);
+    /**
+     * Get the values of a per-DOF variable.
+     *
+     * @param context   the context in which to execute this kernel
+     * @param variable  the index of the variable to get
+     * @param values    on exit, this contains the values
+     */
+    void getPerDofVariable(ContextImpl& context, int variable, std::vector<Vec3>& values) const;
+    /**
+     * Set the values of a per-DOF variable.
+     *
+     * @param context   the context in which to execute this kernel
+     * @param variable  the index of the variable to get
+     * @param values    a vector containing the values
+     */
+    void setPerDofVariable(ContextImpl& context, int variable, const std::vector<Vec3>& values);
+private:
+    class ReorderListener;
+    class GlobalTarget;
+    class DerivFunction;
+    std::string createPerDofComputation(const std::string& variable, const Lepton::ParsedExpression& expr, CustomIntegrator& integrator,
+        const std::string& forceName, const std::string& energyName, std::vector<const TabulatedFunction*>& functions,
+        std::vector<std::pair<std::string, std::string> >& functionNames);
+    void prepareForComputation(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid);
+    Lepton::ExpressionTreeNode replaceDerivFunctions(const Lepton::ExpressionTreeNode& node, OpenMM::ContextImpl& context);
+    void findExpressionsForDerivs(const Lepton::ExpressionTreeNode& node, std::vector<std::pair<Lepton::ExpressionTreeNode, std::string> >& variableNodes);
+    void recordGlobalValue(double value, GlobalTarget target, CustomIntegrator& integrator);
+    void recordChangedParameters(ContextImpl& context);
+    bool evaluateCondition(int step);
+    ComputeContext& cc;
+    double energy;
+    float energyFloat;
+    int numGlobalVariables, sumWorkGroupSize;
+    bool hasInitializedKernels, deviceGlobalsAreCurrent, modifiesParameters, hasAnyConstraints, needsEnergyParamDerivs;
+    std::vector<bool> deviceValuesAreCurrent;
+    mutable std::vector<bool> localValuesAreCurrent;
+    ComputeArray globalValues, sumBuffer, summedValue;
+    ComputeArray uniformRandoms, randomSeed, perDofEnergyParamDerivs;
+    std::vector<ComputeArray> tabulatedFunctions, perDofValues;
+    std::map<int, double> savedEnergy;
+    std::map<int, ComputeArray> savedForces;
+    std::set<int> validSavedForces;
+    mutable std::vector<std::vector<mm_float4> > localPerDofValuesFloat;
+    mutable std::vector<std::vector<mm_double4> > localPerDofValuesDouble;
+    std::map<std::string, double> energyParamDerivs;
+    std::vector<std::string> perDofEnergyParamDerivNames;
+    std::vector<double> localPerDofEnergyParamDerivs;
+    std::vector<double> localGlobalValues;
+    std::vector<double> initialGlobalVariables;
+    std::vector<std::vector<ComputeKernel> > kernels;
+    ComputeKernel randomKernel, kineticEnergyKernel, sumKineticEnergyKernel;
+    std::vector<CustomIntegrator::ComputationType> stepType;
+    std::vector<CustomIntegratorUtilities::Comparison> comparisons;
+    std::vector<std::vector<Lepton::CompiledExpression> > globalExpressions;
+    CompiledExpressionSet expressionSet;
+    std::vector<bool> needsGlobals, needsForces, needsEnergy;
+    std::vector<bool> computeBothForceAndEnergy, invalidatesForces, merged;
+    std::vector<int> forceGroupFlags, blockEnd, requiredGaussian, requiredUniform;
+    std::vector<int> stepEnergyVariableIndex, globalVariableIndex, parameterVariableIndex;
+    int gaussianVariableIndex, uniformVariableIndex, dtVariableIndex;
+    std::vector<std::string> parameterNames;
+    std::vector<GlobalTarget> stepTarget;
+};
+
+class CommonIntegrateCustomStepKernel::GlobalTarget {
+public:
+    CommonIntegrateCustomStepKernel::GlobalTargetType type;
+    int variableIndex;
+    GlobalTarget() {
+    }
+    GlobalTarget(CommonIntegrateCustomStepKernel::GlobalTargetType type, int variableIndex) : type(type), variableIndex(variableIndex) {
+    }
+};
+
+/**
+ * This kernel is invoked to remove center of mass motion from the system.
+ */
+class CommonRemoveCMMotionKernel : public RemoveCMMotionKernel {
+public:
+    CommonRemoveCMMotionKernel(std::string name, const Platform& platform, ComputeContext& cc) : RemoveCMMotionKernel(name, platform), cc(cc) {
+    }
+    /**
+     * Initialize the kernel, setting up the particle masses.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CMMotionRemover this kernel will be used for
+     */
+    void initialize(const System& system, const CMMotionRemover& force);
+    /**
+     * Execute the kernel.
+     *
+     * @param context    the context in which to execute this kernel
+     */
+    void execute(ContextImpl& context);
+private:
+    ComputeContext& cc;
+    int frequency;
+    ComputeArray cmMomentum;
+    ComputeKernel kernel1, kernel2;
+};
+
+/**
+ * This kernel is invoked by RMSDForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CommonCalcRMSDForceKernel : public CalcRMSDForceKernel {
+public:
+    CommonCalcRMSDForceKernel(std::string name, const Platform& platform, ComputeContext& cc) : CalcRMSDForceKernel(name, platform), cc(cc) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the RMSDForce this kernel will be used for
+     */
+    void initialize(const System& system, const RMSDForce& force);
+    /**
+     * Record the reference positions and particle indices.
+     */
+    void recordParameters(const RMSDForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * This is the internal implementation of execute(), templatized on whether we're
+     * using single or double precision.
+     */
+    template <class REAL>
+    double executeImpl(ContextImpl& context);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the RMSDForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const RMSDForce& force);
+private:
+    class ForceInfo;
+    ComputeContext& cc;
+    ForceInfo* info;
+    int blockSize;
+    double sumNormRef;
+    ComputeArray referencePos, particles, buffer;
+    ComputeKernel kernel1, kernel2;
+};
+
+/**
+ * This kernel is invoked by AndersenThermostat at the start of each time step to adjust the particle velocities.
+ */
+class CommonApplyAndersenThermostatKernel : public ApplyAndersenThermostatKernel {
+public:
+    CommonApplyAndersenThermostatKernel(std::string name, const Platform& platform, ComputeContext& cc) : ApplyAndersenThermostatKernel(name, platform), cc(cc) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param thermostat the AndersenThermostat this kernel will be used for
+     */
+    void initialize(const System& system, const AndersenThermostat& thermostat);
+    /**
+     * Execute the kernel.
+     *
+     * @param context    the context in which to execute this kernel
+     */
+    void execute(ContextImpl& context);
+private:
+    ComputeContext& cc;
+    int randomSeed;
+    ComputeArray atomGroups;
+    ComputeKernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_COMMONKERNELS_H_*/

platforms/common/include/openmm/common/ComputeArray.h

+#ifndef OPENMM_COMPUTEARRAY_H_
+#define OPENMM_COMPUTEARRAY_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: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/common/ArrayInterface.h"
+
+namespace OpenMM {
+
+/**
+ * This is an implementation of ArrayInterface that acts as a wrapper around a platform-specific
+ * array implementation (typically CudaArray or OpenCLArray).  This class can be used in code that
+ * is not platform-specific, and an appropriate implementation array is created automatically
+ * based on the ComputeContext.
+ */
+
+class OPENMM_EXPORT_COMMON ComputeArray : public ArrayInterface {
+public:
+    /**
+     * Create an uninitialized ComputeArray object.  It cannot be used until initialize() is called on it.
+     */
+    ComputeArray();
+    /**
+     * Release all resources allocated by this object.
+     */
+    ~ComputeArray();
+    /**
+     * Get the internal array this object is wrapping.
+     */
+    ArrayInterface& getArray();
+    /**
+     * Initialize this array.
+     *
+     * @param context           the context for which to create the array
+     * @param size              the number of elements in the array
+     * @param elementSize       the size of each element in bytes
+     * @param name              the name of the array
+     */
+    void initialize(ComputeContext& context, int size, int elementSize, const std::string& name);
+    /**
+     * Initialize this object.  The template argument is the data type of each array element.
+     *
+     * @param context           the context for which to create the array
+     * @param size              the number of elements in the array
+     * @param name              the name of the array
+     */
+    template <class T>
+    void initialize(ComputeContext& context, int size, const std::string& name) {
+        initialize(context, size, sizeof(T), name);
+    }
+    /**
+     * Recreate the internal storage to have a different size.
+     */
+    void resize(int size);
+    /**
+     * Get whether this array has been initialized.
+     */
+    bool isInitialized() const;
+    /**
+     * Get the number of elements in the array.
+     */
+    int getSize() const;
+    /**
+     * Get the size of each element in bytes.
+     */
+    int getElementSize() const;
+    /**
+     * Get the name of the array.
+     */
+    const std::string& getName() const;
+    /**
+     * Get the context this array belongs to.
+     */
+    ComputeContext& getContext();
+    /**
+     * Copy the values in a vector to the Buffer.
+     */
+    template <class T>
+    void upload(const std::vector<T>& data, bool convert=false) {
+        ArrayInterface::upload(data, convert);
+    }
+    /**
+     * Copy the values in the Buffer to a vector.
+     */
+    template <class T>
+    void download(std::vector<T>& data) const {
+        ArrayInterface::download(data);
+    }
+    /**
+     * Copy the values from host memory to the array.
+     * 
+     * @param data     the data to copy
+     * @param blocking if true, this call will block until the transfer is complete.  Subclasses often
+     *                 have restrictions on non-blocking copies, such as that the source data must be
+     *                 in page-locked memory.
+     */
+    void upload(const void* data, bool blocking=true);
+    /**
+     * Copy the values in the array to host memory.
+     * 
+     * @param data     the destination to copy the value to
+     * @param blocking if true, this call will block until the transfer is complete.  Subclasses often
+     *                 have restrictions on non-blocking copies, such as that the destination must be
+     *                 in page-locked memory.
+     */
+    void download(void* data, bool blocking=true) const;
+    /**
+     * Copy the values in this array to a second array.
+     * 
+     * @param dest     the destination array to copy to
+     */
+    void copyTo(ArrayInterface& dest) const;
+private:
+    ArrayInterface* impl;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_COMPUTEARRAY_H_*/

platforms/common/include/openmm/common/ComputeContext.h

+#ifndef OPENMM_COMPUTECONTEXT_H_
+#define OPENMM_COMPUTECONTEXT_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: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#ifdef _MSC_VER
+    // Prevent Windows from defining macros that interfere with other code.
+    #define NOMINMAX
+#endif
+#include "openmm/common/ArrayInterface.h"
+#include "openmm/common/BondedUtilities.h"
+#include "openmm/common/ComputeEvent.h"
+#include "openmm/common/ComputeForceInfo.h"
+#include "openmm/common/ComputeProgram.h"
+#include "openmm/common/ComputeVectorTypes.h"
+#include "openmm/common/IntegrationUtilities.h"
+#include "openmm/common/NonbondedUtilities.h"
+#include "openmm/Vec3.h"
+#include <pthread.h>
+#include <map>
+#include <queue>
+#include <string>
+#include <vector>
+
+namespace OpenMM {
+
+class ExpressionUtilities;
+class System;
+class ThreadPool;
+
+/**
+ * This abstract class defines the interface by which platforms compile and execute
+ * kernels.  It also manages the arrays use for storing standard information, like
+ * positions and forces.
+ */
+
+class OPENMM_EXPORT_COMMON ComputeContext {
+public:
+    class WorkTask;
+    class WorkThread;
+    class ReorderListener;
+    class ForcePreComputation;
+    class ForcePostComputation;
+    ComputeContext(const System& system);
+    virtual ~ComputeContext();
+    /**
+     * Add a ComputeForceInfo to this context.  Force kernels call this during initialization
+     * to provide information about particular forces.
+     */
+    virtual void addForce(ComputeForceInfo* force);
+    /**
+     * Get all ComputeForceInfos that have been added to this context.
+     */
+    std::vector<ComputeForceInfo*>& getForceInfos() {
+        return forces;
+    }
+    /**
+     * Request that the context provide at least a particular number of force buffers.
+     * This is only meaningful for devices that do not support 64 bit atomic operations.
+     * On other devices, this will typically have no effect.  Force kernels should call
+     * this during initialization.
+     */
+    virtual void requestForceBuffers(int minBuffers) {
+    }
+    /**
+     * Set this as the current context for the calling thread.  This should be called before
+     * doing any computation when you do not know what other code has just been executing on
+     * the thread.  Platforms that rely on binding contexts to threads (such as CUDA) need to
+     * implement this.
+     */
+    virtual void setAsCurrent() {
+    }
+    /**
+     * Get the number of contexts being used for the current simulation.
+     * This is relevant when a simulation is parallelized across multiple devices.  In that case,
+     * one ComputeContext is created for each device.
+     */
+    virtual int getNumContexts() const = 0;
+    /**
+     * Get the index of this context in the list of ones being used for the current simulation.
+     * This is relevant when a simulation is parallelized across multiple devices.  In that case,
+     * one ComputeContext is created for each device.
+     */
+    virtual int getContextIndex() const = 0;
+    /**
+     * Construct an uninitialized array of the appropriate class for this platform.  The returned
+     * value should be created on the heap with the "new" operator.
+     */
+    virtual ArrayInterface* createArray() = 0;
+    /**
+     * Construct a ComputeEvent object of the appropriate class for this platform.
+     */
+    virtual ComputeEvent createEvent() = 0;
+    /**
+     * Compile source code to create a ComputeProgram.
+     *
+     * @param source             the source code of the program
+     * @param defines            a set of preprocessor definitions (name, value) to define when compiling the program
+     */
+    virtual ComputeProgram compileProgram(const std::string source, const std::map<std::string, std::string>& defines=std::map<std::string, std::string>()) = 0;
+    /**
+     * Set all elements of an array to 0.
+     */
+    virtual void clearBuffer(ArrayInterface& array) = 0;
+    /**
+     * Register an array that should be automatically cleared (all elements set to 0) at the start of each force or energy computation.
+     */
+    virtual void addAutoclearBuffer(ArrayInterface& array) = 0;
+    /**
+     * Get whether the device being used is a CPU.  In some cases, different algorithms
+     * may be more efficient on CPUs and GPUs.
+     */
+    virtual bool getIsCPU() const = 0;
+    /**
+     * Get the SIMD width of the device being used.
+     */
+    virtual int getSIMDWidth() const = 0;
+    /**
+     * Get whether the device being used supports 64 bit atomic operations on global memory.
+     */
+    virtual bool getSupports64BitGlobalAtomics() const = 0;
+    /**
+     * Get whether the device being used supports double precision math.
+     */
+    virtual bool getSupportsDoublePrecision() const = 0;
+    /**
+     * Get whether double precision is being used.
+     */
+    virtual bool getUseDoublePrecision() const = 0;
+    /**
+     * Get whether mixed precision is being used.
+     */
+    virtual bool getUseMixedPrecision() const = 0;
+    /**
+     * Get the current simulation time.
+     */
+    double getTime() {
+        return time;
+    }
+    /**
+     * Set the current simulation time.
+     */
+    void setTime(double t) {
+        time = t;
+    }
+    /**
+     * Get the number of integration steps that have been taken.
+     */
+    int getStepCount() {
+        return stepCount;
+    }
+    /**
+     * Set the number of integration steps that have been taken.
+     */
+    void setStepCount(int steps) {
+        stepCount = steps;
+    }
+    /**
+     * Get the number of times forces or energy has been computed.
+     */
+    int getComputeForceCount() {
+        return computeForceCount;
+    }
+    /**
+     * Set the number of times forces or energy has been computed.
+     */
+    void setComputeForceCount(int count) {
+        computeForceCount = count;
+    }
+    /**
+     * Get the number of time steps since the atoms were reordered.
+     */
+    int getStepsSinceReorder() const {
+        return stepsSinceReorder;
+    }
+    /**
+     * Set the number of time steps since the atoms were reordered.
+     */
+    void setStepsSinceReorder(int steps) {
+        stepsSinceReorder = steps;
+    }
+    /**
+     * Get whether atoms were reordered during the most recent force/energy computation.
+     */
+    bool getAtomsWereReordered() const {
+        return atomsWereReordered;
+    }
+    /**
+     * Set whether atoms were reordered during the most recent force/energy computation.
+     */
+    void setAtomsWereReordered(bool wereReordered) {
+        atomsWereReordered = wereReordered;
+    }
+    /**
+     * Reorder the internal arrays of atoms to try to keep spatially contiguous atoms close
+     * together in the arrays.
+     */
+    void reorderAtoms();
+    /**
+     * Add a listener that should be called whenever atoms get reordered.  The OpenCLContext
+     * assumes ownership of the object, and deletes it when the context itself is deleted.
+     */
+    void addReorderListener(ReorderListener* listener);
+    /**
+     * Get the list of ReorderListeners.
+     */
+    std::vector<ReorderListener*>& getReorderListeners() {
+        return reorderListeners;
+    }
+    /**
+     * Add a pre-computation that should be called at the very start of force and energy evaluations.
+     * The OpenCLContext assumes ownership of the object, and deletes it when the context itself is deleted.
+     */
+    void addPreComputation(ForcePreComputation* computation);
+    /**
+     * Get the list of ForcePreComputations.
+     */
+    std::vector<ForcePreComputation*>& getPreComputations() {
+        return preComputations;
+    }
+    /**
+     * Add a post-computation that should be called at the very end of force and energy evaluations.
+     * The OpenCLContext assumes ownership of the object, and deletes it when the context itself is deleted.
+     */
+    void addPostComputation(ForcePostComputation* computation);
+    /**
+     * Get the list of ForcePostComputations.
+     */
+    std::vector<ForcePostComputation*>& getPostComputations() {
+        return postComputations;
+    }
+    /**
+     * Get the flag that marks whether the current force evaluation is valid.
+     */
+    bool getForcesValid() const {
+        return forcesValid;
+    }
+    /**
+     * Get the flag that marks whether the current force evaluation is valid.
+     */
+    void setForcesValid(bool valid) {
+        forcesValid = valid;
+    }
+    /**
+     * Get the number of atoms.
+     */
+    int getNumAtoms() const {
+        return numAtoms;
+    }
+    /**
+     * Get the number of atoms, rounded up to a multiple of TileSize.  This is the actual size of
+     * most arrays with one element per atom.
+     */
+    int getPaddedNumAtoms() const {
+        return paddedNumAtoms;
+    }
+    /**
+     * Get the standard number of thread blocks to use when executing kernels.
+     */
+    virtual int getNumThreadBlocks() const = 0;
+    /**
+     * Get the maximum number of threads in a thread block supported by this device.
+     */
+    virtual int getMaxThreadBlockSize() const = 0;
+    /**
+     * Get the array which contains the position (the xyz components) and charge (the w component) of each atom.
+     */
+    virtual ArrayInterface& getPosq() = 0;
+    /**
+     * Get the array which contains a correction to the position of each atom.  This only exists if getUseMixedPrecision() returns true.
+     */
+    virtual ArrayInterface& getPosqCorrection() = 0;
+    /**
+     * Get the array which contains the velocity (the xyz components) and inverse mass (the w component) of each atom.
+     */
+    virtual ArrayInterface& getVelm() = 0;
+    /**
+     * On devices that do not support 64 bit atomics, this returns an array containing buffers of type real4 in which
+     * forces can be accumulated.  Do not call this if getSupports64BitGlobalAtomics() returns true.  The returned value
+     * in that case is undefined, and it may throw an exception.
+     */
+    virtual ArrayInterface& getForceBuffers() = 0;
+    /**
+     * Get the array which contains a contribution to each force represented as 64 bit fixed point.
+     */
+    virtual ArrayInterface& getLongForceBuffer() = 0;
+    /**
+     * Get the array which contains the buffer in which energy is computed.
+     */
+    virtual ArrayInterface& getEnergyBuffer() = 0;
+    /**
+     * Get the array which contains the buffer in which derivatives of the energy with respect to parameters are computed.
+     */
+    virtual ArrayInterface& getEnergyParamDerivBuffer() = 0;
+    /**
+     * Get a pointer to a block of pinned memory that can be used for asynchronous transfers between host and device.
+     * This is guaranteed to be at least as large as any of the arrays returned by methods of this class.
+     * 
+     * Because this buffer is freely available to all code, care is needed to avoid conflicts.  Only access this
+     * buffer from the main thread, and make sure all transfers are complete before you invoke any other code that
+     * might make use of it
+     */
+    virtual void* getPinnedBuffer() = 0;
+    /**
+     * Get a shared ThreadPool that code can use to parallelize operations.
+     * 
+     * Because this object is freely available to all code, care is needed to avoid conflicts.  Only use it
+     * from the main thread, and make sure all operations are complete before you invoke any other code that
+     * might make use of it
+     */
+    virtual ThreadPool& getThreadPool() = 0;
+    /**
+     * Get the host-side vector which contains the index of each atom.
+     */
+    const std::vector<int>& getAtomIndex() const {
+        return atomIndex;
+    }
+    /**
+     * Get the array which contains the index of each atom.
+     */
+    virtual ArrayInterface& getAtomIndexArray() = 0;
+    /**
+     * Get the number of cells by which the positions are offset.
+     */
+    std::vector<mm_int4>& getPosCellOffsets() {
+        return posCellOffsets;
+    }
+    /**
+     * Replace all occurrences of a list of substrings.
+     *
+     * @param input   a string to process
+     * @param replacements a set of strings that should be replaced with new strings wherever they appear in the input string
+     * @return a new string produced by performing the replacements
+     */
+    std::string replaceStrings(const std::string& input, const std::map<std::string, std::string>& replacements) const;
+    /**
+     * Convert a number to a string in a format suitable for including in a kernel.
+     * This takes into account whether the context uses single or double precision.
+     */
+    std::string doubleToString(double value) const;
+    /**
+     * Convert a number to a string in a format suitable for including in a kernel.
+     */
+    std::string intToString(int value) const;
+    /**
+     * Get whether the periodic box is triclinic.
+     */
+    virtual bool getBoxIsTriclinic() const = 0;
+    /**
+     * Get the vectors defining the periodic box.
+     */
+    virtual void getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) const = 0;
+    /**
+     * Set the vectors defining the periodic box.
+     */
+    virtual void setPeriodicBoxVectors(const Vec3& a, const Vec3& b, const Vec3& c) = 0; 
+    /**
+     * Get the IntegrationUtilities for this context.
+     */
+    virtual IntegrationUtilities& getIntegrationUtilities() = 0;
+    /**
+     * Get the ExpressionUtilities for this context.
+     */
+    virtual ExpressionUtilities& getExpressionUtilities() = 0;
+    /**
+     * Get the BondedUtilities for this context.
+     */
+    virtual BondedUtilities& getBondedUtilities() = 0;
+    /**
+     * Get the NonbondedUtilities for this context.
+     */
+    virtual NonbondedUtilities& getNonbondedUtilities() = 0;
+    /**
+     * This should be called by the Integrator from its own initialize() method.
+     * It ensures all contexts are fully initialized.
+     */
+    virtual void initializeContexts() = 0;
+    /**
+     * Get the thread used by this context for executing parallel computations.
+     */
+    WorkThread& getWorkThread() {
+        return *thread;
+    }
+    /**
+     * Get the names of all parameters with respect to which energy derivatives are computed.
+     */
+    virtual const std::vector<std::string>& getEnergyParamDerivNames() const = 0;
+    /**
+     * Get a workspace data structure used for accumulating the values of derivatives of the energy
+     * with respect to parameters.
+     */
+    virtual std::map<std::string, double>& getEnergyParamDerivWorkspace() = 0;
+    /**
+     * Register that the derivative of potential energy with respect to a context parameter
+     * will need to be calculated.  If this is called multiple times for a single parameter,
+     * it is only added to the list once.
+     * 
+     * @param param    the name of the parameter to add
+     */
+    virtual void addEnergyParameterDerivative(const std::string& param) = 0;
+    /**
+     * Mark that the current molecule definitions (and hence the atom order) may be invalid.
+     * This should be called whenever force field parameters change.  It will cause the definitions
+     * and order to be revalidated.
+     * 
+     * If you know which force has changed, calling the alternate form that takes a ComputeForceInfo
+     * is more efficient.
+     */
+    void invalidateMolecules();
+    /**
+     * Mark that the current molecule definitions from one particular force (and hence the atom order)
+     * may be invalid.  This should be called whenever force field parameters change.  It will cause the
+     * definitions and order to be revalidated.
+     */
+    bool invalidateMolecules(ComputeForceInfo* force);
+    /**
+     * Wait until all work that has been queued (kernel executions, asynchronous data transfers, etc.)
+     * has been submitted to the device.  This does not mean it has necessarily been completed.
+     * Calling this periodically may improve the responsiveness of the computer's GUI, but at the
+     * expense of reduced simulation performance.
+     */
+    virtual void flushQueue() = 0;
+protected:
+    struct Molecule;
+    struct MoleculeGroup;
+    class VirtualSiteInfo;
+    void findMoleculeGroups();
+    /**
+     * This is the internal implementation of reorderAtoms(), templatized by the numerical precision in use.
+     */
+    template <class Real, class Real4, class Mixed, class Mixed4>
+    void reorderAtomsImpl();
+    const System& system;
+    double time;
+    int numAtoms, paddedNumAtoms, stepCount, computeForceCount, stepsSinceReorder;
+    bool atomsWereReordered, forcesValid;
+    std::vector<ComputeForceInfo*> forces;
+    std::vector<Molecule> molecules;
+    std::vector<MoleculeGroup> moleculeGroups;
+    std::vector<int> atomIndex;
+    std::vector<mm_int4> posCellOffsets;
+    std::vector<ReorderListener*> reorderListeners;
+    std::vector<ForcePreComputation*> preComputations;
+    std::vector<ForcePostComputation*> postComputations;
+    WorkThread* thread;
+};
+
+struct ComputeContext::Molecule {
+    std::vector<int> atoms;
+    std::vector<int> constraints;
+    std::vector<std::vector<int> > groups;
+};
+
+struct ComputeContext::MoleculeGroup {
+    std::vector<int> atoms;
+    std::vector<int> instances;
+    std::vector<int> offsets;
+};
+
+/**
+ * This abstract class defines a task to be executed on the worker thread.
+ */
+class OPENMM_EXPORT_COMMON ComputeContext::WorkTask {
+public:
+    virtual void execute() = 0;
+    virtual ~WorkTask() {
+    }
+};
+
+class OPENMM_EXPORT_COMMON ComputeContext::WorkThread {
+public:
+    struct ThreadData;
+    WorkThread();
+    ~WorkThread();
+    /**
+     * Request that a task be executed on the worker thread.  The argument should have been allocated on the
+     * heap with the "new" operator.  After its execute() method finishes, the object will be deleted automatically.
+     */
+    void addTask(ComputeContext::WorkTask* task);
+    /**
+     * Get whether the worker thread is idle, waiting for a task to be added.
+     */
+    bool isWaiting();
+    /**
+     * Get whether the worker thread has exited.
+     */
+    bool isFinished();
+    /**
+     * Block until all tasks have finished executing and the worker thread is idle.
+     */
+    void flush();
+private:
+    std::queue<ComputeContext::WorkTask*> tasks;
+    bool waiting, finished;
+    pthread_mutex_t queueLock;
+    pthread_cond_t waitForTaskCondition, queueEmptyCondition;
+    pthread_t thread;
+};
+
+/**
+ * This abstract class defines a function to be executed whenever atoms get reordered.
+ * Objects that need to know when reordering happens should create a ReorderListener
+ * and register it by calling addReorderListener().
+ */
+class OPENMM_EXPORT_COMMON ComputeContext::ReorderListener {
+public:
+    virtual void execute() = 0;
+    virtual ~ReorderListener() {
+    }
+};
+
+/**
+ * This abstract class defines a function to be executed at the very beginning of force and
+ * energy evaluation, before any other calculation has been done.  It is useful for operations
+ * that need to be performed at a nonstandard point in the process.  After creating a
+ * ForcePreComputation, register it by calling addForcePreComputation().
+ */
+class OPENMM_EXPORT_COMMON ComputeContext::ForcePreComputation {
+public:
+    virtual ~ForcePreComputation() {
+    }
+    /**
+     * @param includeForce  true if forces should be computed
+     * @param includeEnergy true if potential energy should be computed
+     * @param groups        a set of bit flags for which force groups to include
+     */
+    virtual void computeForceAndEnergy(bool includeForces, bool includeEnergy, int groups) = 0;
+};
+
+/**
+ * This abstract class defines a function to be executed at the very end of force and
+ * energy evaluation, after all other calculations have been done.  It is useful for operations
+ * that need to be performed at a nonstandard point in the process.  After creating a
+ * ForcePostComputation, register it by calling addForcePostComputation().
+ */
+class OPENMM_EXPORT_COMMON ComputeContext::ForcePostComputation {
+public:
+    virtual ~ForcePostComputation() {
+    }
+    /**
+     * @param includeForce  true if forces should be computed
+     * @param includeEnergy true if potential energy should be computed
+     * @param groups        a set of bit flags for which force groups to include
+     * @return an optional contribution to add to the potential energy.
+     */
+    virtual double computeForceAndEnergy(bool includeForces, bool includeEnergy, int groups) = 0;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_COMPUTECONTEXT_H_*/