Implemented parameter derivatives for CustomCompoundBondForce and CustomCentroidBondForce

openmmapi/include/openmm/CustomCentroidBondForce.h

@@ -98,6 +98,10 @@ namespace OpenMM {
  * bondParameters.push_back(k);
  * force->addBond(bondGroups, bondParameters);
  * </pre></tt>
+ * 
+ * This class also has the ability to compute derivatives of the potential energy with respect to global parameters.
+ * Call addEnergyParameterDerivative() to request that the derivative with respect to a particular parameter be
+ * computed.  You can then query its value in a Context by calling getState() on it.
  *
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
  * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
@@ -150,6 +154,13 @@ public:
     int getNumGlobalParameters() const {
         return globalParameters.size();
     }
+    /**
+     * Get the number of global parameters with respect to which the derivative of the energy
+     * should be computed.
+     */
+    int getNumEnergyParameterDerivatives() const {
+        return energyParameterDerivatives.size();
+    }
     /**
      * Get the number of tabulated functions that have been defined.
      */
@@ -229,6 +240,21 @@ public:
      * @param defaultValue   the default value of the parameter
      */
     void setGlobalParameterDefaultValue(int index, double defaultValue);
+    /**
+     * Request that this Force compute the derivative of its energy with respect to a global parameter.
+     * The parameter must have already been added with addGlobalParameter().
+     *
+     * @param name             the name of the parameter
+     */
+    void addEnergyParameterDerivative(const std::string& name);
+    /**
+     * Get the name of a global parameter with respect to which this Force should compute the
+     * derivative of the energy.
+     *
+     * @param index     the index of the parameter derivative, between 0 and getNumEnergyParameterDerivatives()
+     * @return the parameter name
+     */
+    const std::string& getEnergyParameterDerivativeName(int index) const;
     /**
      * Add a particle group.
      *
@@ -351,6 +377,7 @@ private:
     std::vector<GroupInfo> groups;
     std::vector<BondInfo> bonds;
     std::vector<FunctionInfo> functions;
+    std::vector<int> energyParameterDerivatives;
     bool usePeriodic;
 };
 

openmmapi/include/openmm/CustomCompoundBondForce.h

@@ -87,6 +87,10 @@ namespace OpenMM {
  * force->addPerBondParameter("theta0");
  * force->addPerBondParameter("r0");
  * </pre></tt>
+ * 
+ * This class also has the ability to compute derivatives of the potential energy with respect to global parameters.
+ * Call addEnergyParameterDerivative() to request that the derivative with respect to a particular parameter be
+ * computed.  You can then query its value in a Context by calling getState() on it.
  *
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
  * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
@@ -133,6 +137,13 @@ public:
     int getNumGlobalParameters() const {
         return globalParameters.size();
     }
+    /**
+     * Get the number of global parameters with respect to which the derivative of the energy
+     * should be computed.
+     */
+    int getNumEnergyParameterDerivatives() const {
+        return energyParameterDerivatives.size();
+    }
     /**
      * Get the number of tabulated functions that have been defined.
      */
@@ -212,6 +223,21 @@ public:
      * @param defaultValue   the default value of the parameter
      */
     void setGlobalParameterDefaultValue(int index, double defaultValue);
+    /**
+     * Request that this Force compute the derivative of its energy with respect to a global parameter.
+     * The parameter must have already been added with addGlobalParameter().
+     *
+     * @param name             the name of the parameter
+     */
+    void addEnergyParameterDerivative(const std::string& name);
+    /**
+     * Get the name of a global parameter with respect to which this Force should compute the
+     * derivative of the energy.
+     *
+     * @param index     the index of the parameter derivative, between 0 and getNumEnergyParameterDerivatives()
+     * @return the parameter name
+     */
+    const std::string& getEnergyParameterDerivativeName(int index) const;
     /**
      * Add a bond to the force
      *
@@ -323,6 +349,7 @@ private:
     std::vector<GlobalParameterInfo> globalParameters;
     std::vector<BondInfo> bonds;
     std::vector<FunctionInfo> functions;
+    std::vector<int> energyParameterDerivatives;
     bool usePeriodic;
 };
 

openmmapi/src/CustomCentroidBondForce.cpp

@@ -104,6 +104,20 @@ void CustomCentroidBondForce::setGlobalParameterDefaultValue(int index, double d
     globalParameters[index].defaultValue = defaultValue;
 }
 
+void CustomCentroidBondForce::addEnergyParameterDerivative(const string& name) {
+    for (int i = 0; i < globalParameters.size(); i++)
+        if (name == globalParameters[i].name) {
+            energyParameterDerivatives.push_back(i);
+            return;
+        }
+    throw OpenMMException(string("addEnergyParameterDerivative: Unknown global parameter '"+name+"'"));
+}
+
+const string& CustomCentroidBondForce::getEnergyParameterDerivativeName(int index) const {
+    ASSERT_VALID_INDEX(index, energyParameterDerivatives);
+    return globalParameters[energyParameterDerivatives[index]].name;
+}
+
 int CustomCentroidBondForce::addGroup(const vector<int>& particles, const vector<double>& weights) {
     if (particles.size() != weights.size() && weights.size() > 0)
         throw OpenMMException("CustomCentroidBondForce: wrong number of weights specified for a group.");

openmmapi/src/CustomCompoundBondForce.cpp

@@ -105,6 +105,20 @@ void CustomCompoundBondForce::setGlobalParameterDefaultValue(int index, double d
     globalParameters[index].defaultValue = defaultValue;
 }
 
+void CustomCompoundBondForce::addEnergyParameterDerivative(const string& name) {
+    for (int i = 0; i < globalParameters.size(); i++)
+        if (name == globalParameters[i].name) {
+            energyParameterDerivatives.push_back(i);
+            return;
+        }
+    throw OpenMMException(string("addEnergyParameterDerivative: Unknown global parameter '"+name+"'"));
+}
+
+const string& CustomCompoundBondForce::getEnergyParameterDerivativeName(int index) const {
+    ASSERT_VALID_INDEX(index, energyParameterDerivatives);
+    return globalParameters[energyParameterDerivatives[index]].name;
+}
+
 int CustomCompoundBondForce::addBond(const vector<int>& particles, const vector<double>& parameters) {
     if (particles.size() != particlesPerBond)
         throw OpenMMException("CustomCompoundBondForce: wrong number of particles specified for a bond.");

platforms/opencl/src/OpenCLKernels.cpp

@@ -5220,6 +5220,12 @@ void OpenCLCalcCustomCompoundBondForceKernel::initialize(const System& system, c
         compute<<buffer.getType()<<" bondParams"<<(i+1)<<" = "<<argName<<"[index];\n";
     }
     forceExpressions["energy += "] = energyExpression;
+    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
+        string paramName = force.getEnergyParameterDerivativeName(i);
+        string derivVariable = cl.getBondedUtilities().addEnergyParameterDerivative(paramName);
+        Lepton::ParsedExpression derivExpression = energyExpression.differentiate(paramName).optimize();
+        forceExpressions[derivVariable+" += "] = derivExpression;
+    }
     compute << cl.getExpressionUtilities().createExpressions(forceExpressions, variables, functionList, functionDefinitions, "temp");
 
     // Finally, apply forces to atoms.

platforms/reference/include/ReferenceCustomCentroidBondIxn.h

@@ -26,7 +26,7 @@
 #define __ReferenceCustomCentroidBondIxn_H__
 
 #include "ReferenceBondIxn.h"
-#include "lepton/ExpressionProgram.h"
+#include "openmm/internal/CompiledExpressionSet.h"
 #include "lepton/ParsedExpression.h"
 #include <map>
 #include <vector>
@@ -44,12 +44,15 @@ class ReferenceCustomCentroidBondIxn : public ReferenceBondIxn {
       std::vector<std::vector<int> > groupAtoms;
       std::vector<std::vector<double> > normalizedWeights;
       std::vector<std::vector<int> > bondGroups;
-      Lepton::ExpressionProgram energyExpression;
-      std::vector<std::string> bondParamNames;
+      CompiledExpressionSet expressionSet;
+      Lepton::CompiledExpression energyExpression;
+      std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
+      std::vector<int> bondParamIndex;
       std::vector<PositionTermInfo> positionTerms;
       std::vector<DistanceTermInfo> distanceTerms;
       std::vector<AngleTermInfo> angleTerms;
       std::vector<DihedralTermInfo> dihedralTerms;
+      int numParameters;
       bool usePeriodic;
       RealVec boxVectors[3];
 
@@ -60,15 +63,13 @@ class ReferenceCustomCentroidBondIxn : public ReferenceBondIxn {
 
          @param bond             the index of the bond
          @param groupCenters     group center coordinates
-         @param variables        the values of variables that may appear in expressions
          @param forces           force array (forces added)
          @param totalEnergy      total energy
 
          --------------------------------------------------------------------------------------- */
 
       void calculateOneIxn(int bond, std::vector<OpenMM::RealVec>& groupCenters,
-                           std::map<std::string, double>& variables, std::vector<OpenMM::RealVec>& forces,
-                           RealOpenMM* totalEnergy) const;
+                           std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, double* energyParamDerivs);
 
       void computeDelta(int group1, int group2, RealOpenMM* delta, std::vector<OpenMM::RealVec>& groupCenters) const;
 
@@ -86,7 +87,8 @@ class ReferenceCustomCentroidBondIxn : public ReferenceBondIxn {
        ReferenceCustomCentroidBondIxn(int numGroupsPerBond, const std::vector<std::vector<int> >& groupAtoms,
                                const std::vector<std::vector<double> >& normalizedWeights, const std::vector<std::vector<int> >& bondGroups, const Lepton::ParsedExpression& energyExpression,
                                const std::vector<std::string>& bondParameterNames, const std::map<std::string, std::vector<int> >& distances,
-                               const std::map<std::string, std::vector<int> >& angles, const std::map<std::string, std::vector<int> >& dihedrals);
+                               const std::map<std::string, std::vector<int> >& angles, const std::map<std::string, std::vector<int> >& dihedrals,
+                               const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions);
 
       /**---------------------------------------------------------------------------------------
 
@@ -130,7 +132,7 @@ class ReferenceCustomCentroidBondIxn : public ReferenceBondIxn {
 
       void calculatePairIxn(std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** bondParameters,
                             const std::map<std::string, double>& globalParameters,
-                            std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy) const;
+                            std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, double* energyParamDerivs);
 
 // ---------------------------------------------------------------------------------------
 
@@ -139,9 +141,9 @@ class ReferenceCustomCentroidBondIxn : public ReferenceBondIxn {
 class ReferenceCustomCentroidBondIxn::PositionTermInfo {
 public:
     std::string name;
-    int group, component;
-    Lepton::ExpressionProgram forceExpression;
-    PositionTermInfo(const std::string& name, int group, int component, const Lepton::ExpressionProgram& forceExpression) :
+    int group, component, index;
+    Lepton::CompiledExpression forceExpression;
+    PositionTermInfo(const std::string& name, int group, int component, const Lepton::CompiledExpression& forceExpression) :
             name(name), group(group), component(component), forceExpression(forceExpression) {
     }
 };
@@ -149,10 +151,10 @@ public:
 class ReferenceCustomCentroidBondIxn::DistanceTermInfo {
 public:
     std::string name;
-    int g1, g2;
-    Lepton::ExpressionProgram forceExpression;
+    int g1, g2, index;
+    Lepton::CompiledExpression forceExpression;
     mutable RealOpenMM delta[ReferenceForce::LastDeltaRIndex];
-    DistanceTermInfo(const std::string& name, const std::vector<int>& groups, const Lepton::ExpressionProgram& forceExpression) :
+    DistanceTermInfo(const std::string& name, const std::vector<int>& groups, const Lepton::CompiledExpression& forceExpression) :
             name(name), g1(groups[0]), g2(groups[1]), forceExpression(forceExpression) {
     }
 };
@@ -160,11 +162,11 @@ public:
 class ReferenceCustomCentroidBondIxn::AngleTermInfo {
 public:
     std::string name;
-    int g1, g2, g3;
-    Lepton::ExpressionProgram forceExpression;
+    int g1, g2, g3, index;
+    Lepton::CompiledExpression forceExpression;
     mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
     mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
-    AngleTermInfo(const std::string& name, const std::vector<int>& groups, const Lepton::ExpressionProgram& forceExpression) :
+    AngleTermInfo(const std::string& name, const std::vector<int>& groups, const Lepton::CompiledExpression& forceExpression) :
             name(name), g1(groups[0]), g2(groups[1]), g3(groups[2]), forceExpression(forceExpression) {
     }
 };
@@ -172,14 +174,14 @@ public:
 class ReferenceCustomCentroidBondIxn::DihedralTermInfo {
 public:
     std::string name;
-    int g1, g2, g3, g4;
-    Lepton::ExpressionProgram forceExpression;
+    int g1, g2, g3, g4, index;
+    Lepton::CompiledExpression forceExpression;
     mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
     mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
     mutable RealOpenMM delta3[ReferenceForce::LastDeltaRIndex];
     mutable RealOpenMM cross1[3];
     mutable RealOpenMM cross2[3];
-    DihedralTermInfo(const std::string& name, const std::vector<int>& groups, const Lepton::ExpressionProgram& forceExpression) :
+    DihedralTermInfo(const std::string& name, const std::vector<int>& groups, const Lepton::CompiledExpression& forceExpression) :
             name(name), g1(groups[0]), g2(groups[1]), g3(groups[2]), g4(groups[3]), forceExpression(forceExpression) {
     }
 };

platforms/reference/include/ReferenceCustomCompoundBondIxn.h

@@ -26,7 +26,7 @@
 #define __ReferenceCustomCompoundBondIxn_H__
 
 #include "ReferenceBondIxn.h"
-#include "lepton/ExpressionProgram.h"
+#include "openmm/internal/CompiledExpressionSet.h"
 #include "lepton/ParsedExpression.h"
 #include <map>
 #include <vector>
@@ -42,12 +42,15 @@ class ReferenceCustomCompoundBondIxn : public ReferenceBondIxn {
       class AngleTermInfo;
       class DihedralTermInfo;
       std::vector<std::vector<int> > bondAtoms;
-      Lepton::ExpressionProgram energyExpression;
-      std::vector<std::string> bondParamNames;
+      CompiledExpressionSet expressionSet;
+      Lepton::CompiledExpression energyExpression;
+      std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
+      std::vector<int> bondParamIndex;
       std::vector<ParticleTermInfo> particleTerms;
       std::vector<DistanceTermInfo> distanceTerms;
       std::vector<AngleTermInfo> angleTerms;
       std::vector<DihedralTermInfo> dihedralTerms;
+      int numParameters;
       bool usePeriodic;
       RealVec boxVectors[3];
 
@@ -58,15 +61,13 @@ class ReferenceCustomCompoundBondIxn : public ReferenceBondIxn {
 
          @param bond             the index of the bond
          @param atomCoordinates  atom coordinates
-         @param variables        the values of variables that may appear in expressions
          @param forces           force array (forces added)
          @param totalEnergy      total energy
 
          --------------------------------------------------------------------------------------- */
 
       void calculateOneIxn(int bond, std::vector<OpenMM::RealVec>& atomCoordinates,
-                           std::map<std::string, double>& variables, std::vector<OpenMM::RealVec>& forces,
-                           RealOpenMM* totalEnergy) const;
+                           std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, double* energyParamDerivs);
 
       void computeDelta(int atom1, int atom2, RealOpenMM* delta, std::vector<OpenMM::RealVec>& atomCoordinates) const;
 
@@ -83,7 +84,8 @@ class ReferenceCustomCompoundBondIxn : public ReferenceBondIxn {
 
        ReferenceCustomCompoundBondIxn(int numParticlesPerBond, const std::vector<std::vector<int> >& bondAtoms, const Lepton::ParsedExpression& energyExpression,
                                const std::vector<std::string>& bondParameterNames, const std::map<std::string, std::vector<int> >& distances,
-                               const std::map<std::string, std::vector<int> >& angles, const std::map<std::string, std::vector<int> >& dihedrals);
+                               const std::map<std::string, std::vector<int> >& angles, const std::map<std::string, std::vector<int> >& dihedrals,
+                               const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions);
 
       /**---------------------------------------------------------------------------------------
 
@@ -127,7 +129,7 @@ class ReferenceCustomCompoundBondIxn : public ReferenceBondIxn {
 
       void calculatePairIxn(std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** bondParameters,
                             const std::map<std::string, double>& globalParameters,
-                            std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy) const;
+                            std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, double* energyParamDerivs);
 
 // ---------------------------------------------------------------------------------------
 
@@ -136,9 +138,9 @@ class ReferenceCustomCompoundBondIxn : public ReferenceBondIxn {
 class ReferenceCustomCompoundBondIxn::ParticleTermInfo {
 public:
     std::string name;
-    int atom, component;
-    Lepton::ExpressionProgram forceExpression;
-    ParticleTermInfo(const std::string& name, int atom, int component, const Lepton::ExpressionProgram& forceExpression) :
+    int atom, component, index;
+    Lepton::CompiledExpression forceExpression;
+    ParticleTermInfo(const std::string& name, int atom, int component, const Lepton::CompiledExpression& forceExpression) :
             name(name), atom(atom), component(component), forceExpression(forceExpression) {
     }
 };
@@ -146,10 +148,10 @@ public:
 class ReferenceCustomCompoundBondIxn::DistanceTermInfo {
 public:
     std::string name;
-    int p1, p2;
-    Lepton::ExpressionProgram forceExpression;
+    int p1, p2, index;
+    Lepton::CompiledExpression forceExpression;
     mutable RealOpenMM delta[ReferenceForce::LastDeltaRIndex];
-    DistanceTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
+    DistanceTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression) :
             name(name), p1(atoms[0]), p2(atoms[1]), forceExpression(forceExpression) {
     }
 };
@@ -157,11 +159,11 @@ public:
 class ReferenceCustomCompoundBondIxn::AngleTermInfo {
 public:
     std::string name;
-    int p1, p2, p3;
-    Lepton::ExpressionProgram forceExpression;
+    int p1, p2, p3, index;
+    Lepton::CompiledExpression forceExpression;
     mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
     mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
-    AngleTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
+    AngleTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression) :
             name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), forceExpression(forceExpression) {
     }
 };
@@ -169,14 +171,14 @@ public:
 class ReferenceCustomCompoundBondIxn::DihedralTermInfo {
 public:
     std::string name;
-    int p1, p2, p3, p4;
-    Lepton::ExpressionProgram forceExpression;
+    int p1, p2, p3, p4, index;
+    Lepton::CompiledExpression forceExpression;
     mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
     mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
     mutable RealOpenMM delta3[ReferenceForce::LastDeltaRIndex];
     mutable RealOpenMM cross1[3];
     mutable RealOpenMM cross2[3];
-    DihedralTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
+    DihedralTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression) :
             name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), p4(atoms[3]), forceExpression(forceExpression) {
     }
 };

platforms/reference/include/ReferenceKernels.h

@@ -878,7 +878,7 @@ private:
     int numBonds, numParticles;
     RealOpenMM **bondParamArray;
     ReferenceCustomCentroidBondIxn* ixn;
-    std::vector<std::string> globalParameterNames;
+    std::vector<std::string> globalParameterNames, energyParamDerivNames;
     bool usePeriodic;
 };
 
@@ -917,7 +917,7 @@ private:
     int numBonds;
     RealOpenMM **bondParamArray;
     ReferenceCustomCompoundBondIxn* ixn;
-    std::vector<std::string> globalParameterNames;
+    std::vector<std::string> globalParameterNames, energyParamDerivNames;
     bool usePeriodic;
 };
 

platforms/reference/src/ReferenceKernels.cpp

@@ -1757,7 +1757,13 @@ void ReferenceCalcCustomCentroidBondForceKernel::initialize(const System& system
         bondParameterNames.push_back(force.getPerBondParameterName(i));
     for (int i = 0; i < force.getNumGlobalParameters(); i++)
         globalParameterNames.push_back(force.getGlobalParameterName(i));
-    ixn = new ReferenceCustomCentroidBondIxn(force.getNumGroupsPerBond(), groupAtoms, normalizedWeights, bondGroups, energyExpression, bondParameterNames, distances, angles, dihedrals);
+    vector<Lepton::CompiledExpression> energyParamDerivExpressions;
+    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
+        string param = force.getEnergyParameterDerivativeName(i);
+        energyParamDerivNames.push_back(param);
+        energyParamDerivExpressions.push_back(energyExpression.differentiate(param).createCompiledExpression());
+    }
+    ixn = new ReferenceCustomCentroidBondIxn(force.getNumGroupsPerBond(), groupAtoms, normalizedWeights, bondGroups, energyExpression, bondParameterNames, distances, angles, dihedrals, energyParamDerivExpressions);
 
     // Delete the custom functions.
 
@@ -1774,7 +1780,11 @@ double ReferenceCalcCustomCentroidBondForceKernel::execute(ContextImpl& context,
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
     if (usePeriodic)
         ixn->setPeriodic(extractBoxVectors(context));
-    ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
+    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
+    ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
+    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
+    for (int i = 0; i < energyParamDerivNames.size(); i++)
+        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
     return energy;
 }
 
@@ -1837,7 +1847,13 @@ void ReferenceCalcCustomCompoundBondForceKernel::initialize(const System& system
         bondParameterNames.push_back(force.getPerBondParameterName(i));
     for (int i = 0; i < force.getNumGlobalParameters(); i++)
         globalParameterNames.push_back(force.getGlobalParameterName(i));
-    ixn = new ReferenceCustomCompoundBondIxn(force.getNumParticlesPerBond(), bondParticles, energyExpression, bondParameterNames, distances, angles, dihedrals);
+    vector<Lepton::CompiledExpression> energyParamDerivExpressions;
+    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
+        string param = force.getEnergyParameterDerivativeName(i);
+        energyParamDerivNames.push_back(param);
+        energyParamDerivExpressions.push_back(energyExpression.differentiate(param).createCompiledExpression());
+    }
+    ixn = new ReferenceCustomCompoundBondIxn(force.getNumParticlesPerBond(), bondParticles, energyExpression, bondParameterNames, distances, angles, dihedrals, energyParamDerivExpressions);
 
     // Delete the custom functions.
 
@@ -1854,7 +1870,11 @@ double ReferenceCalcCustomCompoundBondForceKernel::execute(ContextImpl& context,
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
     if (usePeriodic)
         ixn->setPeriodic(extractBoxVectors(context));
-    ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
+    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
+    ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
+    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
+    for (int i = 0; i < energyParamDerivNames.size(); i++)
+        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
     return energy;
 }
 

platforms/reference/src/SimTKReference/ReferenceCustomCentroidBondIxn.cpp

@@ -40,23 +40,47 @@ using namespace OpenMM;
 ReferenceCustomCentroidBondIxn::ReferenceCustomCentroidBondIxn(int numGroupsPerBond, const vector<vector<int> >& groupAtoms,
             const vector<vector<double> >& normalizedWeights, const vector<vector<int> >& bondGroups,
             const Lepton::ParsedExpression& energyExpression, const vector<string>& bondParameterNames,
-            const map<string, vector<int> >& distances, const map<string, vector<int> >& angles, const map<string, vector<int> >& dihedrals) :
-            groupAtoms(groupAtoms), normalizedWeights(normalizedWeights), bondGroups(bondGroups), energyExpression(energyExpression.createProgram()), bondParamNames(bondParameterNames), usePeriodic(false) {
+            const map<string, vector<int> >& distances, const map<string, vector<int> >& angles, const map<string, vector<int> >& dihedrals,
+            const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions) :
+            groupAtoms(groupAtoms), normalizedWeights(normalizedWeights), bondGroups(bondGroups), energyExpression(energyExpression.createCompiledExpression()),
+            usePeriodic(false), energyParamDerivExpressions(energyParamDerivExpressions) {
+    expressionSet.registerExpression(this->energyExpression);
+    for (int i = 0; i < this->energyParamDerivExpressions.size(); i++)
+        expressionSet.registerExpression(this->energyParamDerivExpressions[i]);
     for (int i = 0; i < numGroupsPerBond; i++) {
         stringstream xname, yname, zname;
         xname << 'x' << (i+1);
         yname << 'y' << (i+1);
         zname << 'z' << (i+1);
-        positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(xname.str(), i, 0, energyExpression.differentiate(xname.str()).optimize().createProgram()));
-        positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(yname.str(), i, 1, energyExpression.differentiate(yname.str()).optimize().createProgram()));
-        positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(zname.str(), i, 2, energyExpression.differentiate(zname.str()).optimize().createProgram()));
+        positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(xname.str(), i, 0, energyExpression.differentiate(xname.str()).createCompiledExpression()));
+        positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(yname.str(), i, 1, energyExpression.differentiate(yname.str()).createCompiledExpression()));
+        positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(zname.str(), i, 2, energyExpression.differentiate(zname.str()).createCompiledExpression()));
     }
     for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter)
-        distanceTerms.push_back(ReferenceCustomCentroidBondIxn::DistanceTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
+        distanceTerms.push_back(ReferenceCustomCentroidBondIxn::DistanceTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).createCompiledExpression()));
     for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter)
-        angleTerms.push_back(ReferenceCustomCentroidBondIxn::AngleTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
+        angleTerms.push_back(ReferenceCustomCentroidBondIxn::AngleTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).createCompiledExpression()));
     for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter)
-        dihedralTerms.push_back(ReferenceCustomCentroidBondIxn::DihedralTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
+        dihedralTerms.push_back(ReferenceCustomCentroidBondIxn::DihedralTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).createCompiledExpression()));
+    for (int i = 0; i < positionTerms.size(); i++) {
+        expressionSet.registerExpression(positionTerms[i].forceExpression);
+        positionTerms[i].index = expressionSet.getVariableIndex(positionTerms[i].name);
+    }
+    for (int i = 0; i < distanceTerms.size(); i++) {
+        expressionSet.registerExpression(distanceTerms[i].forceExpression);
+        distanceTerms[i].index = expressionSet.getVariableIndex(distanceTerms[i].name);
+    }
+    for (int i = 0; i < angleTerms.size(); i++) {
+        expressionSet.registerExpression(angleTerms[i].forceExpression);
+        angleTerms[i].index = expressionSet.getVariableIndex(angleTerms[i].name);
+    }
+    for (int i = 0; i < dihedralTerms.size(); i++) {
+        expressionSet.registerExpression(dihedralTerms[i].forceExpression);
+        dihedralTerms[i].index = expressionSet.getVariableIndex(dihedralTerms[i].name);
+    }
+    numParameters = bondParameterNames.size();
+    for (int i = 0; i < numParameters; i++)
+        bondParamIndex.push_back(expressionSet.getVariableIndex(bondParameterNames[i]));
 }
 
 ReferenceCustomCentroidBondIxn::~ReferenceCustomCentroidBondIxn() {
@@ -71,7 +95,7 @@ void ReferenceCustomCentroidBondIxn::setPeriodic(OpenMM::RealVec* vectors) {
 
 void ReferenceCustomCentroidBondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates, RealOpenMM** bondParameters,
                                              const map<string, double>& globalParameters, vector<RealVec>& forces,
-                                             RealOpenMM* totalEnergy) const {
+                                             RealOpenMM* totalEnergy, double* energyParamDerivs) {
 
     // First compute the center of each group.
 
@@ -84,13 +108,14 @@ void ReferenceCustomCentroidBondIxn::calculatePairIxn(vector<RealVec>& atomCoord
 
     // Compute the forces on groups.
 
-    map<string, double> variables = globalParameters;
+    for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter)
+        expressionSet.setVariable(expressionSet.getVariableIndex(iter->first), iter->second);
     vector<RealVec> groupForces(numGroups);
     int numBonds = bondGroups.size();
     for (int bond = 0; bond < numBonds; bond++) {
-        for (int j = 0; j < (int) bondParamNames.size(); j++)
-            variables[bondParamNames[j]] = bondParameters[bond][j];
-        calculateOneIxn(bond, groupCenters, variables, groupForces, totalEnergy);
+        for (int i = 0; i < numParameters; i++)
+            expressionSet.setVariable(bondParamIndex[i], bondParameters[bond][i]);
+        calculateOneIxn(bond, groupCenters, groupForces, totalEnergy, energyParamDerivs);
     }
 
     // Apply the forces to the individual atoms.
@@ -102,31 +127,24 @@ void ReferenceCustomCentroidBondIxn::calculatePairIxn(vector<RealVec>& atomCoord
 }
 
 void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<RealVec>& groupCenters,
-                        map<string, double>& variables, vector<RealVec>& forces, RealOpenMM* totalEnergy) const {
-
-    // ---------------------------------------------------------------------------------------
-
-    static const std::string methodName = "\nReferenceCustomCentroidBondIxn::calculateOneIxn";
-
-    // ---------------------------------------------------------------------------------------
-
+                        vector<RealVec>& forces, RealOpenMM* totalEnergy, double* energyParamDerivs) {
     // Compute all of the variables the energy can depend on.
 
     const vector<int>& groups = bondGroups[bond];
     for (int i = 0; i < (int) positionTerms.size(); i++) {
         const PositionTermInfo& term = positionTerms[i];
-        variables[term.name] = groupCenters[groups[term.group]][term.component];
+        expressionSet.setVariable(term.index, groupCenters[groups[term.group]][term.component]);
     }
     for (int i = 0; i < (int) distanceTerms.size(); i++) {
         const DistanceTermInfo& term = distanceTerms[i];
         computeDelta(groups[term.g1], groups[term.g2], term.delta, groupCenters);
-        variables[term.name] = term.delta[ReferenceForce::RIndex];
+        expressionSet.setVariable(term.index, term.delta[ReferenceForce::RIndex]);
     }
     for (int i = 0; i < (int) angleTerms.size(); i++) {
         const AngleTermInfo& term = angleTerms[i];
         computeDelta(groups[term.g1], groups[term.g2], term.delta1, groupCenters);
         computeDelta(groups[term.g3], groups[term.g2], term.delta2, groupCenters);
-        variables[term.name] = computeAngle(term.delta1, term.delta2);
+        expressionSet.setVariable(term.index, computeAngle(term.delta1, term.delta2));
     }
     for (int i = 0; i < (int) dihedralTerms.size(); i++) {
         const DihedralTermInfo& term = dihedralTerms[i];
@@ -135,21 +153,21 @@ void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<RealVec>&
         computeDelta(groups[term.g4], groups[term.g3], term.delta3, groupCenters);
         RealOpenMM dotDihedral, signOfDihedral;
         RealOpenMM* crossProduct[] = {term.cross1, term.cross2};
-        variables[term.name] = getDihedralAngleBetweenThreeVectors(term.delta1, term.delta2, term.delta3, crossProduct, &dotDihedral, term.delta1, &signOfDihedral, 1);
+        expressionSet.setVariable(term.index, getDihedralAngleBetweenThreeVectors(term.delta1, term.delta2, term.delta3, crossProduct, &dotDihedral, term.delta1, &signOfDihedral, 1));
     }
 
     // Apply forces based on individual particle coordinates.
 
     for (int i = 0; i < (int) positionTerms.size(); i++) {
         const PositionTermInfo& term = positionTerms[i];
-        forces[groups[term.group]][term.component] -= term.forceExpression.evaluate(variables);
+        forces[groups[term.group]][term.component] -= term.forceExpression.evaluate();
     }
 
     // Apply forces based on distances.
 
     for (int i = 0; i < (int) distanceTerms.size(); i++) {
         const DistanceTermInfo& term = distanceTerms[i];
-        RealOpenMM dEdR = (RealOpenMM) (term.forceExpression.evaluate(variables)/(term.delta[ReferenceForce::RIndex]));
+        RealOpenMM dEdR = (RealOpenMM) (term.forceExpression.evaluate()/(term.delta[ReferenceForce::RIndex]));
         for (int i = 0; i < 3; i++) {
            RealOpenMM force  = -dEdR*term.delta[i];
            forces[groups[term.g1]][i] -= force;
@@ -161,7 +179,7 @@ void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<RealVec>&
 
     for (int i = 0; i < (int) angleTerms.size(); i++) {
         const AngleTermInfo& term = angleTerms[i];
-        RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate(variables);
+        RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate();
         RealOpenMM thetaCross[ReferenceForce::LastDeltaRIndex];
         SimTKOpenMMUtilities::crossProductVector3(term.delta1, term.delta2, thetaCross);
         RealOpenMM lengthThetaCross = SQRT(DOT3(thetaCross, thetaCross));
@@ -188,7 +206,7 @@ void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<RealVec>&
 
     for (int i = 0; i < (int) dihedralTerms.size(); i++) {
         const DihedralTermInfo& term = dihedralTerms[i];
-        RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate(variables);
+        RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate();
         RealOpenMM internalF[4][3];
         RealOpenMM forceFactors[4];
         RealOpenMM normCross1 = DOT3(term.cross1, term.cross1);
@@ -218,7 +236,12 @@ void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<RealVec>&
     // Add the energy
 
     if (totalEnergy)
-        *totalEnergy += (RealOpenMM) energyExpression.evaluate(variables);
+        *totalEnergy += (RealOpenMM) energyExpression.evaluate();
+    
+    // Compute derivatives of the energy.
+    
+    for (int i = 0; i < energyParamDerivExpressions.size(); i++)
+        energyParamDerivs[i] += energyParamDerivExpressions[i].evaluate();
 }
 
 void ReferenceCustomCentroidBondIxn::computeDelta(int group1, int group2, RealOpenMM* delta, vector<RealVec>& groupCenters) const {

platforms/reference/src/SimTKReference/ReferenceCustomCompoundBondIxn.cpp

@@ -45,23 +45,47 @@ using namespace OpenMM;
 
 ReferenceCustomCompoundBondIxn::ReferenceCustomCompoundBondIxn(int numParticlesPerBond, const vector<vector<int> >& bondAtoms,
             const Lepton::ParsedExpression& energyExpression, const vector<string>& bondParameterNames,
-            const map<string, vector<int> >& distances, const map<string, vector<int> >& angles, const map<string, vector<int> >& dihedrals) :
-            bondAtoms(bondAtoms), energyExpression(energyExpression.createProgram()), bondParamNames(bondParameterNames), usePeriodic(false) {
+            const map<string, vector<int> >& distances, const map<string, vector<int> >& angles, const map<string, vector<int> >& dihedrals,
+            const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions) :
+            bondAtoms(bondAtoms), energyExpression(energyExpression.createCompiledExpression()), usePeriodic(false),
+            energyParamDerivExpressions(energyParamDerivExpressions) {
+    expressionSet.registerExpression(this->energyExpression);
+    for (int i = 0; i < this->energyParamDerivExpressions.size(); i++)
+        expressionSet.registerExpression(this->energyParamDerivExpressions[i]);
     for (int i = 0; i < numParticlesPerBond; i++) {
         stringstream xname, yname, zname;
         xname << 'x' << (i+1);
         yname << 'y' << (i+1);
         zname << 'z' << (i+1);
-        particleTerms.push_back(ReferenceCustomCompoundBondIxn::ParticleTermInfo(xname.str(), i, 0, energyExpression.differentiate(xname.str()).optimize().createProgram()));
-        particleTerms.push_back(ReferenceCustomCompoundBondIxn::ParticleTermInfo(yname.str(), i, 1, energyExpression.differentiate(yname.str()).optimize().createProgram()));
-        particleTerms.push_back(ReferenceCustomCompoundBondIxn::ParticleTermInfo(zname.str(), i, 2, energyExpression.differentiate(zname.str()).optimize().createProgram()));
+        particleTerms.push_back(ReferenceCustomCompoundBondIxn::ParticleTermInfo(xname.str(), i, 0, energyExpression.differentiate(xname.str()).createCompiledExpression()));
+        particleTerms.push_back(ReferenceCustomCompoundBondIxn::ParticleTermInfo(yname.str(), i, 1, energyExpression.differentiate(yname.str()).createCompiledExpression()));
+        particleTerms.push_back(ReferenceCustomCompoundBondIxn::ParticleTermInfo(zname.str(), i, 2, energyExpression.differentiate(zname.str()).createCompiledExpression()));
     }
     for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter)
-        distanceTerms.push_back(ReferenceCustomCompoundBondIxn::DistanceTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
+        distanceTerms.push_back(ReferenceCustomCompoundBondIxn::DistanceTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).createCompiledExpression()));
     for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter)
-        angleTerms.push_back(ReferenceCustomCompoundBondIxn::AngleTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
+        angleTerms.push_back(ReferenceCustomCompoundBondIxn::AngleTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).createCompiledExpression()));
     for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter)
-        dihedralTerms.push_back(ReferenceCustomCompoundBondIxn::DihedralTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
+        dihedralTerms.push_back(ReferenceCustomCompoundBondIxn::DihedralTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).createCompiledExpression()));
+    for (int i = 0; i < particleTerms.size(); i++) {
+        expressionSet.registerExpression(particleTerms[i].forceExpression);
+        particleTerms[i].index = expressionSet.getVariableIndex(particleTerms[i].name);
+    }
+    for (int i = 0; i < distanceTerms.size(); i++) {
+        expressionSet.registerExpression(distanceTerms[i].forceExpression);
+        distanceTerms[i].index = expressionSet.getVariableIndex(distanceTerms[i].name);
+    }
+    for (int i = 0; i < angleTerms.size(); i++) {
+        expressionSet.registerExpression(angleTerms[i].forceExpression);
+        angleTerms[i].index = expressionSet.getVariableIndex(angleTerms[i].name);
+    }
+    for (int i = 0; i < dihedralTerms.size(); i++) {
+        expressionSet.registerExpression(dihedralTerms[i].forceExpression);
+        dihedralTerms[i].index = expressionSet.getVariableIndex(dihedralTerms[i].name);
+    }
+    numParameters = bondParameterNames.size();
+    for (int i = 0; i < numParameters; i++)
+        bondParamIndex.push_back(expressionSet.getVariableIndex(bondParameterNames[i]));
 }
 
 /**---------------------------------------------------------------------------------------
@@ -94,14 +118,14 @@ void ReferenceCustomCompoundBondIxn::setPeriodic(OpenMM::RealVec* vectors) {
 
 void ReferenceCustomCompoundBondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates, RealOpenMM** bondParameters,
                                              const map<string, double>& globalParameters, vector<RealVec>& forces,
-                                             RealOpenMM* totalEnergy) const {
-
-    map<string, double> variables = globalParameters;
+                                             RealOpenMM* totalEnergy, double* energyParamDerivs) {
+    for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter)
+        expressionSet.setVariable(expressionSet.getVariableIndex(iter->first), iter->second);
     int numBonds = bondAtoms.size();
     for (int bond = 0; bond < numBonds; bond++) {
-        for (int j = 0; j < (int) bondParamNames.size(); j++)
-            variables[bondParamNames[j]] = bondParameters[bond][j];
-        calculateOneIxn(bond, atomCoordinates, variables, forces, totalEnergy);
+        for (int i = 0; i < numParameters; i++)
+            expressionSet.setVariable(bondParamIndex[i], bondParameters[bond][i]);
+        calculateOneIxn(bond, atomCoordinates, forces, totalEnergy, energyParamDerivs);
     }
 }
 
@@ -111,7 +135,6 @@ void ReferenceCustomCompoundBondIxn::calculatePairIxn(vector<RealVec>& atomCoord
 
      @param bond             the index of the bond
      @param atomCoordinates  atom coordinates
-     @param variables        the values of variables that may appear in expressions
      @param forces           force array (forces added)
      @param energyByAtom     atom energy
      @param totalEnergy      total energy
@@ -119,31 +142,24 @@ void ReferenceCustomCompoundBondIxn::calculatePairIxn(vector<RealVec>& atomCoord
      --------------------------------------------------------------------------------------- */
 
 void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<RealVec>& atomCoordinates,
-                        map<string, double>& variables, vector<RealVec>& forces, RealOpenMM* totalEnergy) const {
-
-    // ---------------------------------------------------------------------------------------
-
-    static const std::string methodName = "\nReferenceCustomCompoundBondIxn::calculateOneIxn";
-
-    // ---------------------------------------------------------------------------------------
-
+                        vector<RealVec>& forces, RealOpenMM* totalEnergy, double* energyParamDerivs) {
     // Compute all of the variables the energy can depend on.
 
     const vector<int>& atoms = bondAtoms[bond];
     for (int i = 0; i < (int) particleTerms.size(); i++) {
         const ParticleTermInfo& term = particleTerms[i];
-        variables[term.name] = atomCoordinates[atoms[term.atom]][term.component];
+        expressionSet.setVariable(term.index, atomCoordinates[atoms[term.atom]][term.component]);
     }
     for (int i = 0; i < (int) distanceTerms.size(); i++) {
         const DistanceTermInfo& term = distanceTerms[i];
         computeDelta(atoms[term.p1], atoms[term.p2], term.delta, atomCoordinates);
-        variables[term.name] = term.delta[ReferenceForce::RIndex];
+        expressionSet.setVariable(term.index, term.delta[ReferenceForce::RIndex]);
     }
     for (int i = 0; i < (int) angleTerms.size(); i++) {
         const AngleTermInfo& term = angleTerms[i];
         computeDelta(atoms[term.p1], atoms[term.p2], term.delta1, atomCoordinates);
         computeDelta(atoms[term.p3], atoms[term.p2], term.delta2, atomCoordinates);
-        variables[term.name] = computeAngle(term.delta1, term.delta2);
+        expressionSet.setVariable(term.index, computeAngle(term.delta1, term.delta2));
     }
     for (int i = 0; i < (int) dihedralTerms.size(); i++) {
         const DihedralTermInfo& term = dihedralTerms[i];
@@ -152,21 +168,21 @@ void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<RealVec>&
         computeDelta(atoms[term.p4], atoms[term.p3], term.delta3, atomCoordinates);
         RealOpenMM dotDihedral, signOfDihedral;
         RealOpenMM* crossProduct[] = {term.cross1, term.cross2};
-        variables[term.name] = getDihedralAngleBetweenThreeVectors(term.delta1, term.delta2, term.delta3, crossProduct, &dotDihedral, term.delta1, &signOfDihedral, 1);
+        expressionSet.setVariable(term.index,getDihedralAngleBetweenThreeVectors(term.delta1, term.delta2, term.delta3, crossProduct, &dotDihedral, term.delta1, &signOfDihedral, 1));
     }
     
     // Apply forces based on individual particle coordinates.
     
     for (int i = 0; i < (int) particleTerms.size(); i++) {
         const ParticleTermInfo& term = particleTerms[i];
-        forces[atoms[term.atom]][term.component] -= term.forceExpression.evaluate(variables);
+        forces[atoms[term.atom]][term.component] -= term.forceExpression.evaluate();
     }
 
     // Apply forces based on distances.
 
     for (int i = 0; i < (int) distanceTerms.size(); i++) {
         const DistanceTermInfo& term = distanceTerms[i];
-        RealOpenMM dEdR = (RealOpenMM) (term.forceExpression.evaluate(variables)/(term.delta[ReferenceForce::RIndex]));
+        RealOpenMM dEdR = (RealOpenMM) (term.forceExpression.evaluate()/(term.delta[ReferenceForce::RIndex]));
         for (int i = 0; i < 3; i++) {
            RealOpenMM force  = -dEdR*term.delta[i];
            forces[atoms[term.p1]][i] -= force;
@@ -178,7 +194,7 @@ void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<RealVec>&
 
     for (int i = 0; i < (int) angleTerms.size(); i++) {
         const AngleTermInfo& term = angleTerms[i];
-        RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate(variables);
+        RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate();
         RealOpenMM thetaCross[ReferenceForce::LastDeltaRIndex];
         SimTKOpenMMUtilities::crossProductVector3(term.delta1, term.delta2, thetaCross);
         RealOpenMM lengthThetaCross = SQRT(DOT3(thetaCross, thetaCross));
@@ -205,7 +221,7 @@ void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<RealVec>&
 
     for (int i = 0; i < (int) dihedralTerms.size(); i++) {
         const DihedralTermInfo& term = dihedralTerms[i];
-        RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate(variables);
+        RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate();
         RealOpenMM internalF[4][3];
         RealOpenMM forceFactors[4];
         RealOpenMM normCross1 = DOT3(term.cross1, term.cross1);
@@ -235,7 +251,12 @@ void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<RealVec>&
     // Add the energy
 
     if (totalEnergy)
-        *totalEnergy += (RealOpenMM) energyExpression.evaluate(variables);
+        *totalEnergy += (RealOpenMM) energyExpression.evaluate();
+    
+    // Compute derivatives of the energy.
+    
+    for (int i = 0; i < energyParamDerivExpressions.size(); i++)
+        energyParamDerivs[i] += energyParamDerivExpressions[i].evaluate();
 }
 
 void ReferenceCustomCompoundBondIxn::computeDelta(int atom1, int atom2, RealOpenMM* delta, vector<RealVec>& atomCoordinates) const {

tests/TestCustomCentroidBondForce.h

@@ -333,6 +333,65 @@ void testPeriodic() {
     ASSERT_EQUAL_VEC(Vec3(2*0.5*(5.0/12.0), 0, 0), state.getForces()[4], TOL);
 }
 
+void testEnergyParameterDerivatives() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(2.0);
+    system.addParticle(3.0);
+    system.addParticle(4.0);
+    system.addParticle(5.0);
+    CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "k*(distance(g1,g2)-r0)^2");
+    force->addGlobalParameter("r0", 0.0);
+    force->addGlobalParameter("k", 0.0);
+    force->addEnergyParameterDerivative("r0");
+    force->addEnergyParameterDerivative("k");
+    vector<int> particles1;
+    particles1.push_back(0);
+    particles1.push_back(1);
+    vector<int> particles2;
+    particles2.push_back(2);
+    particles2.push_back(3);
+    particles2.push_back(4);
+    force->addGroup(particles1);
+    force->addGroup(particles2);
+    vector<int> groups;
+    groups.push_back(0);
+    groups.push_back(1);
+    vector<double> parameters;
+    force->addBond(groups, parameters);
+    system.addForce(force);
+
+    // The center of mass of group 0 is (1.5, 0, 0).
+
+    vector<Vec3> positions(5);
+    positions[0] = Vec3(2.5, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+
+    // The center of mass of group 1 is (-1, 0, 0).
+
+    positions[2] = Vec3(-6, 0, 0);
+    positions[3] = Vec3(-1, 0, 0);
+    positions[4] = Vec3(2, 0, 0);
+    
+    // Check the derivatives.
+
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    for (int i = 0; i < 10; i++) {
+        double r0 = 0.1*i;
+        double k = 10-i;
+        context.setParameter("r0", r0);
+        context.setParameter("k", k);
+        State state = context.getState(State::ParameterDerivatives);
+        map<string, double> derivs = state.getEnergyParameterDerivatives();
+        double dEdr0 = -2*k*(2.5-r0);
+        double dEdk = (2.5-r0)*(2.5-r0);
+        ASSERT_EQUAL_TOL(dEdr0, derivs["r0"], 1e-5);
+        ASSERT_EQUAL_TOL(dEdk, derivs["k"], 1e-5);
+    }
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -343,6 +402,7 @@ int main(int argc, char* argv[]) {
         testCustomWeights();
         testIllegalVariable();
         testPeriodic();
+        testEnergyParameterDerivatives();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomCompoundBondForce.h

@@ -406,6 +406,48 @@ void testPeriodic() {
     }
 }
 
+void testEnergyParameterDerivatives() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomCompoundBondForce* custom = new CustomCompoundBondForce(4, "k*(dihedral(p1,p2,p3,p4)-theta0)^2");
+    custom->addGlobalParameter("theta0", 0.0);
+    custom->addGlobalParameter("k", 0.0);
+    custom->addEnergyParameterDerivative("theta0");
+    custom->addEnergyParameterDerivative("k");
+    vector<int> particles(4);
+    particles[0] = 0;
+    particles[1] = 1;
+    particles[2] = 2;
+    particles[3] = 3;
+    vector<double> parameters;
+    custom->addBond(particles, parameters);
+    system.addForce(custom);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 2, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    positions[3] = Vec3(1, 1, 1);
+    context.setPositions(positions);
+    double theta = M_PI/4;
+    for (int i = 0; i < 10; i++) {
+        double theta0 = 0.1*i;
+        double k = 10-i;
+        context.setParameter("theta0", theta0);
+        context.setParameter("k", k);
+        State state = context.getState(State::ParameterDerivatives);
+        map<string, double> derivs = state.getEnergyParameterDerivatives();
+        double dEdtheta0 = -2*k*(theta-theta0);
+        double dEdk = (theta-theta0)*(theta-theta0);
+        ASSERT_EQUAL_TOL(dEdtheta0, derivs["theta0"], 1e-5);
+        ASSERT_EQUAL_TOL(dEdk, derivs["k"], 1e-5);
+    }
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -418,6 +460,7 @@ int main(int argc, char* argv[]) {
         testMultipleBonds();
         testIllegalVariable();
         testPeriodic();
+        testEnergyParameterDerivatives();
         runPlatformTests();
     }
     catch(const exception& e) {