Merge pull request #1547 from peastman/paramderivs

Energy derivatives with respect to global parameters

platforms/reference/include/ReferenceCustomNonbondedIxn.h

 
-/* Portions copyright (c) 2009-2013 Stanford University and Simbios.
+/* Portions copyright (c) 2009-2016 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -27,7 +27,7 @@
 
 #include "ReferencePairIxn.h"
 #include "ReferenceNeighborList.h"
-#include "lepton/CompiledExpression.h"
+#include "openmm/internal/CompiledExpressionSet.h"
 #include <map>
 #include <set>
 #include <utility>
@@ -48,10 +48,10 @@ class ReferenceCustomNonbondedIxn {
       Lepton::CompiledExpression energyExpression;
       Lepton::CompiledExpression forceExpression;
       std::vector<std::string> paramNames;
-      std::vector<double*> energyParticleParams;
-      std::vector<double*> forceParticleParams;
-      double* energyR;
-      double* forceR;
+      std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
+      CompiledExpressionSet expressionSet;
+      std::vector<int> particleParamIndex;
+      int rIndex;
       std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups;
 
       /**---------------------------------------------------------------------------------------
@@ -69,7 +69,7 @@ class ReferenceCustomNonbondedIxn {
          --------------------------------------------------------------------------------------- */
 
       void calculateOneIxn(int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& forces,
-                           RealOpenMM* energyByAtom, RealOpenMM* totalEnergy);
+                           RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, double* energyParamDerivs);
 
 
    public:
@@ -81,7 +81,7 @@ class ReferenceCustomNonbondedIxn {
          --------------------------------------------------------------------------------------- */
 
        ReferenceCustomNonbondedIxn(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpression,
-                                   const std::vector<std::string>& parameterNames);
+                                   const std::vector<std::string>& parameterNames, const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions);
 
       /**---------------------------------------------------------------------------------------
 
@@ -155,7 +155,8 @@ class ReferenceCustomNonbondedIxn {
       void calculatePairIxn(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
                             RealOpenMM** atomParameters, std::vector<std::set<int> >& exclusions,
                             RealOpenMM* fixedParameters, const std::map<std::string, double>& globalParameters,
-                            std::vector<OpenMM::RealVec>& forces, RealOpenMM* energyByAtom, RealOpenMM* totalEnergy);
+                            std::vector<OpenMM::RealVec>& forces, RealOpenMM* energyByAtom, RealOpenMM* totalEnergy,
+                            double* energyParamDerivs);
 
 // ---------------------------------------------------------------------------------------
 

platforms/reference/include/ReferenceCustomTorsionIxn.h

@@ -25,7 +25,7 @@
 #define __ReferenceCustomTorsionIxn_H__
 
 #include "ReferenceBondIxn.h"
-#include "lepton/CompiledExpression.h"
+#include "openmm/internal/CompiledExpressionSet.h"
 
 namespace OpenMM {
 
@@ -34,10 +34,10 @@ class ReferenceCustomTorsionIxn : public ReferenceBondIxn {
    private:
       Lepton::CompiledExpression energyExpression;
       Lepton::CompiledExpression forceExpression;
-      std::vector<double*> energyParams;
-      std::vector<double*> forceParams;
-      double* energyTheta;
-      double* forceTheta;
+      std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
+      CompiledExpressionSet expressionSet;
+      std::vector<int> torsionParamIndex;
+      int thetaIndex;
       int numParameters;
       bool usePeriodic;
       RealVec boxVectors[3];
@@ -51,7 +51,8 @@ class ReferenceCustomTorsionIxn : public ReferenceBondIxn {
          --------------------------------------------------------------------------------------- */
 
        ReferenceCustomTorsionIxn(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpression,
-                              const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters);
+                              const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters,
+                              const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions);
 
       /**---------------------------------------------------------------------------------------
 
@@ -85,7 +86,7 @@ class ReferenceCustomTorsionIxn : public ReferenceBondIxn {
 
       void calculateBondIxn(int* atomIndices, std::vector<OpenMM::RealVec>& atomCoordinates,
                             RealOpenMM* parameters, std::vector<OpenMM::RealVec>& forces,
-                            RealOpenMM* totalEnergy) const;
+                            RealOpenMM* totalEnergy, double* energyParamDerivs);
 
 
 };

platforms/reference/include/ReferenceHarmonicBondIxn.h

@@ -79,7 +79,7 @@ class ReferenceHarmonicBondIxn : public ReferenceBondIxn {
       
       void calculateBondIxn(int* atomIndices, std::vector<OpenMM::RealVec>& atomCoordinates,
                             RealOpenMM* parameters, std::vector<OpenMM::RealVec>& forces,
-                            RealOpenMM* totalEnergy) const;
+                            RealOpenMM* totalEnergy, double* energyParamDerivs);
 
 };
 

platforms/reference/include/ReferenceKernels.h

@@ -38,7 +38,6 @@
 #include "ReferenceNeighborList.h"
 #include "lepton/CompiledExpression.h"
 #include "lepton/CustomFunction.h"
-#include "lepton/ExpressionProgram.h"
 
 namespace OpenMM {
 
@@ -157,6 +156,12 @@ public:
      * @param forces  on exit, this contains the forces
      */
     void getForces(ContextImpl& context, std::vector<Vec3>& forces);
+    /**
+     * Get the current derivatives of the energy with respect to context parameters.
+     *
+     * @param derivs  on exit, this contains the derivatives
+     */
+    void getEnergyParameterDerivatives(ContextImpl& context, std::map<std::string, double>& derivs);
     /**
      * Get the current periodic box vectors.
      *
@@ -319,7 +324,8 @@ private:
     int **bondIndexArray;
     RealOpenMM **bondParamArray;
     Lepton::CompiledExpression energyExpression, forceExpression;
-    std::vector<std::string> parameterNames, globalParameterNames;
+    std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
+    std::vector<std::string> parameterNames, globalParameterNames, energyParamDerivNames;
     bool usePeriodic;
 };
 
@@ -397,7 +403,8 @@ private:
     int **angleIndexArray;
     RealOpenMM **angleParamArray;
     Lepton::CompiledExpression energyExpression, forceExpression;
-    std::vector<std::string> parameterNames, globalParameterNames;
+    std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
+    std::vector<std::string> parameterNames, globalParameterNames, energyParamDerivNames;
     bool usePeriodic;
 };
 
@@ -550,7 +557,8 @@ private:
     int **torsionIndexArray;
     RealOpenMM **torsionParamArray;
     Lepton::CompiledExpression energyExpression, forceExpression;
-    std::vector<std::string> parameterNames, globalParameterNames;
+    std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
+    std::vector<std::string> parameterNames, globalParameterNames, energyParamDerivNames;
     bool usePeriodic;
 };
 
@@ -647,8 +655,10 @@ private:
     std::map<std::string, double> globalParamValues;
     std::vector<std::set<int> > exclusions;
     Lepton::CompiledExpression energyExpression, forceExpression;
-    std::vector<std::string> parameterNames, globalParameterNames;
+    std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
+    std::vector<std::string> parameterNames, globalParameterNames, energyParamDerivNames;
     std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups;
+    std::vector<double> longRangeCoefficientDerivs;
     NonbondedMethod nonbondedMethod;
     NeighborList* neighborList;
 };
@@ -727,14 +737,16 @@ private:
     RealOpenMM **particleParamArray;
     RealOpenMM nonbondedCutoff;
     std::vector<std::set<int> > exclusions;
-    std::vector<std::string> particleParameterNames, globalParameterNames, valueNames;
-    std::vector<Lepton::ExpressionProgram> valueExpressions;
-    std::vector<std::vector<Lepton::ExpressionProgram> > valueDerivExpressions;
-    std::vector<std::vector<Lepton::ExpressionProgram> > valueGradientExpressions;
+    std::vector<std::string> particleParameterNames, globalParameterNames, energyParamDerivNames, valueNames;
+    std::vector<Lepton::CompiledExpression> valueExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > valueDerivExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > valueGradientExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > valueParamDerivExpressions;
     std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
-    std::vector<Lepton::ExpressionProgram> energyExpressions;
-    std::vector<std::vector<Lepton::ExpressionProgram> > energyDerivExpressions;
-    std::vector<std::vector<Lepton::ExpressionProgram> > energyGradientExpressions;
+    std::vector<Lepton::CompiledExpression> energyExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > energyDerivExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > energyGradientExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > energyParamDerivExpressions;
     std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
     NonbondedMethod nonbondedMethod;
     NeighborList* neighborList;
@@ -867,7 +879,7 @@ private:
     int numBonds, numParticles;
     RealOpenMM **bondParamArray;
     ReferenceCustomCentroidBondIxn* ixn;
-    std::vector<std::string> globalParameterNames;
+    std::vector<std::string> globalParameterNames, energyParamDerivNames;
     bool usePeriodic;
 };
 
@@ -906,7 +918,7 @@ private:
     int numBonds;
     RealOpenMM **bondParamArray;
     ReferenceCustomCompoundBondIxn* ixn;
-    std::vector<std::string> globalParameterNames;
+    std::vector<std::string> globalParameterNames, energyParamDerivNames;
     bool usePeriodic;
 };
 

platforms/reference/include/ReferenceLJCoulomb14.h

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -64,7 +64,7 @@ class OPENMM_EXPORT ReferenceLJCoulomb14 : public ReferenceBondIxn {
       
       void calculateBondIxn(int* atomIndices, std::vector<OpenMM::RealVec>& atomCoordinates,
                             RealOpenMM* parameters, std::vector<OpenMM::RealVec>& forces,
-                            RealOpenMM* totalEnergy) const;
+                            RealOpenMM* totalEnergy, double* energyParamDerivs);
 
 };
 

platforms/reference/include/ReferencePlatform.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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -68,6 +68,7 @@ public:
     void* periodicBoxSize;
     void* periodicBoxVectors;
     void* constraints;
+    void* energyParameterDerivatives;
 };
 } // namespace OpenMM
 

platforms/reference/include/ReferenceProperDihedralBond.h

@@ -80,7 +80,7 @@ class OPENMM_EXPORT ReferenceProperDihedralBond : public ReferenceBondIxn {
       
       void calculateBondIxn(int* atomIndices, std::vector<OpenMM::RealVec>& atomCoordinates,
                             RealOpenMM* parameters, std::vector<OpenMM::RealVec>& forces,
-                            RealOpenMM* totalEnergy) const;
+                            RealOpenMM* totalEnergy, double* energyParamDerivs);
 
 };
 

platforms/reference/include/ReferenceRbDihedralBond.h

@@ -78,7 +78,7 @@ class OPENMM_EXPORT ReferenceRbDihedralBond : public ReferenceBondIxn {
       
       void calculateBondIxn(int* atomIndices, std::vector<OpenMM::RealVec>& atomCoordinates,
                             RealOpenMM* parameters, std::vector<OpenMM::RealVec>& forces,
-                            RealOpenMM* totalEnergy) const;
+                            RealOpenMM* totalEnergy, double* energyParamDerivs);
 
 };
 

platforms/reference/src/ReferenceKernels.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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -146,6 +146,11 @@ static ReferenceConstraints& extractConstraints(ContextImpl& context) {
     return *(ReferenceConstraints*) data->constraints;
 }
 
+static map<string, double>& extractEnergyParameterDerivatives(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return *((map<string, double>*) data->energyParameterDerivatives);
+}
+
 /**
  * Make sure an expression doesn't use any undefined variables.
  */
@@ -208,6 +213,8 @@ void ReferenceCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context,
     }
     else
         savedForces = forceData;
+    for (map<string, double>::const_iterator iter = context.getParameters().begin(); iter != context.getParameters().end(); ++iter)
+        extractEnergyParameterDerivatives(context)[iter->first] = 0;
 }
 
 double ReferenceCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForces, bool includeEnergy, int groups, bool& valid) {
@@ -273,6 +280,10 @@ void ReferenceUpdateStateDataKernel::getForces(ContextImpl& context, std::vector
         forces[i] = Vec3(forceData[i][0], forceData[i][1], forceData[i][2]);
 }
 
+void ReferenceUpdateStateDataKernel::getEnergyParameterDerivatives(ContextImpl& context, map<string, double>& derivs) {
+    derivs = extractEnergyParameterDerivatives(context);
+}
+
 void ReferenceUpdateStateDataKernel::getPeriodicBoxVectors(ContextImpl& context, Vec3& a, Vec3& b, Vec3& c) const {
     RealVec* vectors = extractBoxVectors(context);
     a = vectors[0];
@@ -437,6 +448,11 @@ void ReferenceCalcCustomBondForceKernel::initialize(const System& system, const
         parameterNames.push_back(force.getPerBondParameterName(i));
     for (int i = 0; i < force.getNumGlobalParameters(); i++)
         globalParameterNames.push_back(force.getGlobalParameterName(i));
+    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
+        string param = force.getEnergyParameterDerivativeName(i);
+        energyParamDerivNames.push_back(param);
+        energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
+    }
     set<string> variables;
     variables.insert("r");
     variables.insert(parameterNames.begin(), parameterNames.end());
@@ -451,11 +467,15 @@ double ReferenceCalcCustomBondForceKernel::execute(ContextImpl& context, bool in
     map<string, double> globalParameters;
     for (int i = 0; i < (int) globalParameterNames.size(); i++)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
-    ReferenceBondForce refBondForce;
-    ReferenceCustomBondIxn bond(energyExpression, forceExpression, parameterNames, globalParameters);
+    ReferenceCustomBondIxn bond(energyExpression, forceExpression, parameterNames, globalParameters, energyParamDerivExpressions);
     if (usePeriodic)
         bond.setPeriodic(extractBoxVectors(context));
-    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, includeEnergy ? &energy : NULL, bond);
+    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
+    for (int i = 0; i < numBonds; i++)
+        bond.calculateBondIxn(bondIndexArray[i], posData, bondParamArray[i], 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;
 }
 
@@ -562,6 +582,11 @@ void ReferenceCalcCustomAngleForceKernel::initialize(const System& system, const
         parameterNames.push_back(force.getPerAngleParameterName(i));
     for (int i = 0; i < force.getNumGlobalParameters(); i++)
         globalParameterNames.push_back(force.getGlobalParameterName(i));
+    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
+        string param = force.getEnergyParameterDerivativeName(i);
+        energyParamDerivNames.push_back(param);
+        energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
+    }
     set<string> variables;
     variables.insert("theta");
     variables.insert(parameterNames.begin(), parameterNames.end());
@@ -576,11 +601,15 @@ double ReferenceCalcCustomAngleForceKernel::execute(ContextImpl& context, bool i
     map<string, double> globalParameters;
     for (int i = 0; i < (int) globalParameterNames.size(); i++)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
-    ReferenceBondForce refBondForce;
-    ReferenceCustomAngleIxn customAngle(energyExpression, forceExpression, parameterNames, globalParameters);
+    ReferenceCustomAngleIxn customAngle(energyExpression, forceExpression, parameterNames, globalParameters, energyParamDerivExpressions);
     if (usePeriodic)
         customAngle.setPeriodic(extractBoxVectors(context));
-    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, includeEnergy ? &energy : NULL, customAngle);
+    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
+    for (int i = 0; i < numAngles; i++)
+        customAngle.calculateBondIxn(angleIndexArray[i], posData, angleParamArray[i], 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;
 }
 
@@ -823,6 +852,11 @@ void ReferenceCalcCustomTorsionForceKernel::initialize(const System& system, con
         parameterNames.push_back(force.getPerTorsionParameterName(i));
     for (int i = 0; i < force.getNumGlobalParameters(); i++)
         globalParameterNames.push_back(force.getGlobalParameterName(i));
+    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
+        string param = force.getEnergyParameterDerivativeName(i);
+        energyParamDerivNames.push_back(param);
+        energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
+    }
     set<string> variables;
     variables.insert("theta");
     variables.insert(parameterNames.begin(), parameterNames.end());
@@ -837,11 +871,15 @@ double ReferenceCalcCustomTorsionForceKernel::execute(ContextImpl& context, bool
     map<string, double> globalParameters;
     for (int i = 0; i < (int) globalParameterNames.size(); i++)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
-    ReferenceBondForce refBondForce;
-    ReferenceCustomTorsionIxn customTorsion(energyExpression, forceExpression, parameterNames, globalParameters);
+    ReferenceCustomTorsionIxn customTorsion(energyExpression, forceExpression, parameterNames, globalParameters, energyParamDerivExpressions);
     if (usePeriodic)
         customTorsion.setPeriodic(extractBoxVectors(context));
-    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, customTorsion);
+    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
+    for (int i = 0; i < numTorsions; i++)
+        customTorsion.calculateBondIxn(torsionIndexArray[i], posData, torsionParamArray[i], 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;
 }
 
@@ -1088,6 +1126,11 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
         globalParameterNames.push_back(force.getGlobalParameterName(i));
         globalParamValues[force.getGlobalParameterName(i)] = force.getGlobalParameterDefaultValue(i);
     }
+    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
+        string param = force.getEnergyParameterDerivativeName(i);
+        energyParamDerivNames.push_back(param);
+        energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
+    }
     set<string> variables;
     variables.insert("r");
     for (int i = 0; i < numParameters; i++) {
@@ -1127,7 +1170,7 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
     vector<RealVec>& forceData = extractForces(context);
     RealVec* boxVectors = extractBoxVectors(context);
     RealOpenMM energy = 0;
-    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames);
+    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames, energyParamDerivExpressions);
     bool periodic = (nonbondedMethod == CutoffPeriodic);
     if (nonbondedMethod != NoCutoff) {
         computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic, nonbondedCutoff, 0.0);
@@ -1150,15 +1193,22 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
     }
     if (useSwitchingFunction)
         ixn.setUseSwitchingFunction(switchingDistance);
-    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusions, 0, globalParamValues, forceData, 0, includeEnergy ? &energy : NULL);
+    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
+    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusions, 0, globalParamValues, forceData, 0, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
+    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
+    for (int i = 0; i < energyParamDerivNames.size(); i++)
+        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
     
     // Add in the long range correction.
     
     if (!hasInitializedLongRangeCorrection || (globalParamsChanged && forceCopy != NULL)) {
-        longRangeCoefficient = CustomNonbondedForceImpl::calcLongRangeCorrection(*forceCopy, context.getOwner());
+        CustomNonbondedForceImpl::calcLongRangeCorrection(*forceCopy, context.getOwner(), longRangeCoefficient, longRangeCoefficientDerivs);
         hasInitializedLongRangeCorrection = true;
     }
-    energy += longRangeCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
+    double volume = boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2];
+    energy += longRangeCoefficient/volume;
+    for (int i = 0; i < longRangeCoefficientDerivs.size(); i++)
+        energyParamDerivs[energyParamDerivNames[i]] += longRangeCoefficientDerivs[i]/volume;
     return energy;
 }
 
@@ -1180,7 +1230,7 @@ void ReferenceCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImp
     // If necessary, recompute the long range correction.
     
     if (forceCopy != NULL) {
-        longRangeCoefficient = CustomNonbondedForceImpl::calcLongRangeCorrection(force, context.getOwner());
+        CustomNonbondedForceImpl::calcLongRangeCorrection(force, context.getOwner(), longRangeCoefficient, longRangeCoefficientDerivs);
         hasInitializedLongRangeCorrection = true;
         *forceCopy = force;
     }
@@ -1309,6 +1359,7 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
 
     valueDerivExpressions.resize(force.getNumComputedValues());
     valueGradientExpressions.resize(force.getNumComputedValues());
+    valueParamDerivExpressions.resize(force.getNumComputedValues());
     set<string> particleVariables, pairVariables;
     pairVariables.insert("r");
     particleVariables.insert("x");
@@ -1326,21 +1377,26 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
         CustomGBForce::ComputationType type;
         force.getComputedValueParameters(i, name, expression, type);
         Lepton::ParsedExpression ex = Lepton::Parser::parse(expression, functions).optimize();
-        valueExpressions.push_back(ex.createProgram());
+        valueExpressions.push_back(ex.createCompiledExpression());
         valueTypes.push_back(type);
         valueNames.push_back(name);
         if (i == 0) {
-            valueDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
+            valueDerivExpressions[i].push_back(ex.differentiate("r").createCompiledExpression());
             validateVariables(ex.getRootNode(), pairVariables);
         }
         else {
-            valueGradientExpressions[i].push_back(ex.differentiate("x").optimize().createProgram());
-            valueGradientExpressions[i].push_back(ex.differentiate("y").optimize().createProgram());
-            valueGradientExpressions[i].push_back(ex.differentiate("z").optimize().createProgram());
+            valueGradientExpressions[i].push_back(ex.differentiate("x").createCompiledExpression());
+            valueGradientExpressions[i].push_back(ex.differentiate("y").createCompiledExpression());
+            valueGradientExpressions[i].push_back(ex.differentiate("z").createCompiledExpression());
             for (int j = 0; j < i; j++)
-                valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
+                valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).createCompiledExpression());
             validateVariables(ex.getRootNode(), particleVariables);
         }
+        for (int j = 0; j < force.getNumEnergyParameterDerivatives(); j++) {
+            string param = force.getEnergyParameterDerivativeName(j);
+            energyParamDerivNames.push_back(param);
+            valueParamDerivExpressions[i].push_back(ex.differentiate(param).createCompiledExpression());
+        }
         particleVariables.insert(name);
         pairVariables.insert(name+"1");
         pairVariables.insert(name+"2");
@@ -1350,29 +1406,32 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
 
     energyDerivExpressions.resize(force.getNumEnergyTerms());
     energyGradientExpressions.resize(force.getNumEnergyTerms());
+    energyParamDerivExpressions.resize(force.getNumEnergyTerms());
     for (int i = 0; i < force.getNumEnergyTerms(); i++) {
         string expression;
         CustomGBForce::ComputationType type;
         force.getEnergyTermParameters(i, expression, type);
         Lepton::ParsedExpression ex = Lepton::Parser::parse(expression, functions).optimize();
-        energyExpressions.push_back(ex.createProgram());
+        energyExpressions.push_back(ex.createCompiledExpression());
         energyTypes.push_back(type);
         if (type != CustomGBForce::SingleParticle)
-            energyDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
+            energyDerivExpressions[i].push_back(ex.differentiate("r").createCompiledExpression());
         for (int j = 0; j < force.getNumComputedValues(); j++) {
             if (type == CustomGBForce::SingleParticle) {
-                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
-                energyGradientExpressions[i].push_back(ex.differentiate("x").optimize().createProgram());
-                energyGradientExpressions[i].push_back(ex.differentiate("y").optimize().createProgram());
-                energyGradientExpressions[i].push_back(ex.differentiate("z").optimize().createProgram());
+                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).createCompiledExpression());
+                energyGradientExpressions[i].push_back(ex.differentiate("x").createCompiledExpression());
+                energyGradientExpressions[i].push_back(ex.differentiate("y").createCompiledExpression());
+                energyGradientExpressions[i].push_back(ex.differentiate("z").createCompiledExpression());
                 validateVariables(ex.getRootNode(), particleVariables);
             }
             else {
-                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").optimize().createProgram());
-                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").optimize().createProgram());
+                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").createCompiledExpression());
+                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").createCompiledExpression());
                 validateVariables(ex.getRootNode(), pairVariables);
             }
         }
+        for (int j = 0; j < force.getNumEnergyParameterDerivatives(); j++)
+            energyParamDerivExpressions[i].push_back(ex.differentiate(force.getEnergyParameterDerivativeName(j)).createCompiledExpression());
     }
 
     // Delete the custom functions.
@@ -1385,8 +1444,8 @@ double ReferenceCalcCustomGBForceKernel::execute(ContextImpl& context, bool incl
     vector<RealVec>& posData = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
     RealOpenMM energy = 0;
-    ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueNames, valueTypes, energyExpressions,
-        energyDerivExpressions, energyGradientExpressions, energyTypes, particleParameterNames);
+    ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueParamDerivExpressions, valueNames, valueTypes,
+        energyExpressions, energyDerivExpressions, energyGradientExpressions, energyParamDerivExpressions, energyTypes, particleParameterNames);
     bool periodic = (nonbondedMethod == CutoffPeriodic);
     if (periodic)
         ixn.setPeriodic(extractBoxVectors(context));
@@ -1397,7 +1456,11 @@ double ReferenceCalcCustomGBForceKernel::execute(ContextImpl& context, bool incl
     map<string, double> globalParameters;
     for (int i = 0; i < (int) globalParameterNames.size(); i++)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
-    ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, includeEnergy ? &energy : NULL);
+    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
+    ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, 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;
 }
 
@@ -1707,7 +1770,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.
 
@@ -1724,7 +1793,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;
 }
 
@@ -1787,7 +1860,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.
 
@@ -1804,7 +1883,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/ReferencePlatform.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -36,6 +36,7 @@
 #include "openmm/internal/ContextImpl.h"
 #include "SimTKOpenMMRealType.h"
 #include "RealVec.h"
+#include <map>
 #include <vector>
 
 using namespace OpenMM;
@@ -99,6 +100,7 @@ ReferencePlatform::PlatformData::PlatformData(const System& system) : time(0.0),
     periodicBoxSize = new RealVec();
     periodicBoxVectors = new RealVec[3];
     constraints = new ReferenceConstraints(system);
+    energyParameterDerivatives = new map<string, double>();
 }
 
 ReferencePlatform::PlatformData::~PlatformData() {
@@ -108,4 +110,5 @@ ReferencePlatform::PlatformData::~PlatformData() {
     delete (RealVec*) periodicBoxSize;
     delete[] (RealVec*) periodicBoxVectors;
     delete (ReferenceConstraints*) constraints;
+    delete (map<string, double>*) energyParameterDerivatives;
 }

platforms/reference/src/SimTKReference/ReferenceAngleBondIxn.cpp

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -129,7 +129,7 @@ void ReferenceAngleBondIxn::calculateBondIxn(int* atomIndices,
                                              vector<RealVec>& atomCoordinates,
                                              RealOpenMM* parameters,
                                              vector<RealVec>& forces,
-                                             RealOpenMM* totalEnergy) const {
+                                             RealOpenMM* totalEnergy, double* energyParamDerivs) {
 
    // constants -- reduce Visual Studio warnings regarding conversions between float & double
 

platforms/reference/src/SimTKReference/ReferenceBondForce.cpp

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -97,7 +97,7 @@ void ReferenceBondForce::calculateForce(int numberOfBonds, int** atomIndices,
       // calculate bond ixn
 
       referenceBondIxn.calculateBondIxn(atomIndices[ii], atomCoordinates, parameters[ii], 
-                                        forces, totalEnergy);
+                                        forces, totalEnergy, NULL);
    }
 }
 

platforms/reference/src/SimTKReference/ReferenceBondIxn.cpp

 
-/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -78,7 +78,7 @@ ReferenceBondIxn::~ReferenceBondIxn() {
      
    void ReferenceBondIxn::calculateBondIxn(int* atomIndices, vector<RealVec>& atomCoordinates,
                                            RealOpenMM* parameters, vector<RealVec>& forces,
-                                           RealOpenMM* totalEnergy) const {
+                                           RealOpenMM* totalEnergy, double* energyParamDerivs) {
    // ---------------------------------------------------------------------------------------
 
    // static const std::string methodName = "\nReferenceBondIxn::calculateBondIxn";

platforms/reference/src/SimTKReference/ReferenceCMAPTorsionIxn.cpp

@@ -207,5 +207,5 @@ void ReferenceCMAPTorsionIxn::calculateOneIxn(int index, vector<RealVec>& atomCo
    --------------------------------------------------------------------------------------- */
 
 void ReferenceCMAPTorsionIxn::calculateBondIxn(int* atomIndices, vector<RealVec>& atomCoordinates,
-        RealOpenMM* parameters, vector<RealVec>& forces, RealOpenMM* totalEnergy) const {
+        RealOpenMM* parameters, vector<RealVec>& forces, RealOpenMM* totalEnergy, double* energyParamDerivs) {
 }

platforms/reference/src/SimTKReference/ReferenceCustomAngleIxn.cpp

@@ -38,20 +38,19 @@ using namespace std;
    --------------------------------------------------------------------------------------- */
 
 ReferenceCustomAngleIxn::ReferenceCustomAngleIxn(const Lepton::CompiledExpression& energyExpression,
-        const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
-        energyExpression(energyExpression), forceExpression(forceExpression), usePeriodic(false) {
-    
-    energyTheta = ReferenceForce::getVariablePointer(this->energyExpression, "theta");
-    forceTheta = ReferenceForce::getVariablePointer(this->forceExpression, "theta");
+        const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters,
+        const vector<Lepton::CompiledExpression> energyParamDerivExpressions) :
+        energyExpression(energyExpression), forceExpression(forceExpression), usePeriodic(false), energyParamDerivExpressions(energyParamDerivExpressions) {
+    expressionSet.registerExpression(this->energyExpression);
+    expressionSet.registerExpression(this->forceExpression);
+    for (int i = 0; i < this->energyParamDerivExpressions.size(); i++)
+        expressionSet.registerExpression(this->energyParamDerivExpressions[i]);
+    thetaIndex = expressionSet.getVariableIndex("theta");
     numParameters = parameterNames.size();
-    for (int i = 0; i < (int) numParameters; i++) {
-        energyParams.push_back(ReferenceForce::getVariablePointer(this->energyExpression, parameterNames[i]));
-        forceParams.push_back(ReferenceForce::getVariablePointer(this->forceExpression, parameterNames[i]));
-    }
-    for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter) {
-        ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->energyExpression, iter->first), iter->second);
-        ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->forceExpression, iter->first), iter->second);
-    }
+    for (int i = 0; i < (int) numParameters; i++)
+        angleParamIndex.push_back(expressionSet.getVariableIndex(parameterNames[i]));
+    for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter)
+        expressionSet.setVariable(expressionSet.getVariableIndex(iter->first), iter->second);
 }
 
 /**---------------------------------------------------------------------------------------
@@ -86,18 +85,10 @@ void ReferenceCustomAngleIxn::calculateBondIxn(int* atomIndices,
                                                vector<RealVec>& atomCoordinates,
                                                RealOpenMM* parameters,
                                                vector<RealVec>& forces,
-                                               RealOpenMM* totalEnergy) const {
-
-   static const std::string methodName = "\nReferenceCustomAngleIxn::calculateAngleIxn";
-
-   static const RealOpenMM zero        = 0.0;
-   static const RealOpenMM one         = 1.0;
-
+                                               RealOpenMM* totalEnergy, double* energyParamDerivs) {
    RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
-   for (int i = 0; i < numParameters; i++) {
-       ReferenceForce::setVariable(energyParams[i], parameters[i]);
-       ReferenceForce::setVariable(forceParams[i], parameters[i]);
-   }
+   for (int i = 0; i < numParameters; i++)
+       expressionSet.setVariable(angleParamIndex[i], parameters[i]);
 
    // ---------------------------------------------------------------------------------------
 
@@ -122,14 +113,13 @@ void ReferenceCustomAngleIxn::calculateBondIxn(int* atomIndices,
    RealOpenMM dot = DOT3(deltaR[0], deltaR[1]);
    RealOpenMM cosine = dot/SQRT((deltaR[0][ReferenceForce::R2Index]*deltaR[1][ReferenceForce::R2Index]));
    RealOpenMM angle;
-   if (cosine >= one)
-      angle = zero;
-   else if (cosine <= -one)
+   if (cosine >= 1.0)
+      angle = 0.0;
+   else if (cosine <= -1.0)
       angle = PI_M;
    else
       angle = ACOS(cosine);
-   ReferenceForce::setVariable(energyTheta, angle);
-   ReferenceForce::setVariable(forceTheta, angle);
+   expressionSet.setVariable(thetaIndex, angle);
 
    // Compute the force and energy, and apply them to the atoms.
    
@@ -156,6 +146,11 @@ void ReferenceCustomAngleIxn::calculateBondIxn(int* atomIndices,
       }
    }
 
+   // Record parameter derivatives.
+
+   for (int i = 0; i < energyParamDerivExpressions.size(); i++)
+       energyParamDerivs[i] += energyParamDerivExpressions[i].evaluate();
+   
    // accumulate energies
 
    if (totalEnergy != NULL)

platforms/reference/src/SimTKReference/ReferenceCustomBondIxn.cpp

@@ -39,19 +39,19 @@ using namespace OpenMM;
    --------------------------------------------------------------------------------------- */
 
 ReferenceCustomBondIxn::ReferenceCustomBondIxn(const Lepton::CompiledExpression& energyExpression,
-        const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
-        energyExpression(energyExpression), forceExpression(forceExpression), usePeriodic(false) {
-    energyR = ReferenceForce::getVariablePointer(this->energyExpression, "r");
-    forceR = ReferenceForce::getVariablePointer(this->forceExpression, "r");
+        const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters,
+        const vector<Lepton::CompiledExpression> energyParamDerivExpressions) :
+        energyExpression(energyExpression), forceExpression(forceExpression), usePeriodic(false), energyParamDerivExpressions(energyParamDerivExpressions) {
+    expressionSet.registerExpression(this->energyExpression);
+    expressionSet.registerExpression(this->forceExpression);
+    for (int i = 0; i < this->energyParamDerivExpressions.size(); i++)
+        expressionSet.registerExpression(this->energyParamDerivExpressions[i]);
+    rIndex = expressionSet.getVariableIndex("r");
     numParameters = parameterNames.size();
-    for (int i = 0; i < (int) numParameters; i++) {
-        energyParams.push_back(ReferenceForce::getVariablePointer(this->energyExpression, parameterNames[i]));
-        forceParams.push_back(ReferenceForce::getVariablePointer(this->forceExpression, parameterNames[i]));
-    }
-    for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter) {
-        ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->energyExpression, iter->first), iter->second);
-        ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->forceExpression, iter->first), iter->second);
-    }
+    for (int i = 0; i < (int) numParameters; i++)
+        bondParamIndex.push_back(expressionSet.getVariableIndex(parameterNames[i]));
+    for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter)
+        expressionSet.setVariable(expressionSet.getVariableIndex(iter->first), iter->second);
 }
 
 /**---------------------------------------------------------------------------------------
@@ -86,21 +86,10 @@ void ReferenceCustomBondIxn::calculateBondIxn(int* atomIndices,
                                               vector<RealVec>& atomCoordinates,
                                               RealOpenMM* parameters,
                                               vector<RealVec>& forces,
-                                              RealOpenMM* totalEnergy) const {
-
-   static const std::string methodName = "\nReferenceCustomBondIxn::calculateBondIxn";
-
-   static const int twoI               = 2;
-
-   static const RealOpenMM zero        = 0.0;
-   static const RealOpenMM two         = 2.0;
-   static const RealOpenMM half        = 0.5;
-
+                                              RealOpenMM* totalEnergy, double* energyParamDerivs) {
    RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
-   for (int i = 0; i < numParameters; i++) {
-       ReferenceForce::setVariable(energyParams[i], parameters[i]);
-       ReferenceForce::setVariable(forceParams[i], parameters[i]);
-   }
+   for (int i = 0; i < numParameters; i++)
+       expressionSet.setVariable(bondParamIndex[i], parameters[i]);
 
    // ---------------------------------------------------------------------------------------
 
@@ -113,10 +102,9 @@ void ReferenceCustomBondIxn::calculateBondIxn(int* atomIndices,
    else
        ReferenceForce::getDeltaR(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR);
    
-   ReferenceForce::setVariable(energyR, deltaR[ReferenceForce::RIndex]);
-   ReferenceForce::setVariable(forceR, deltaR[ReferenceForce::RIndex]);
+   expressionSet.setVariable(rIndex, deltaR[ReferenceForce::RIndex]);
    RealOpenMM dEdR            = (RealOpenMM) forceExpression.evaluate();
-   dEdR                       = deltaR[ReferenceForce::RIndex] > zero ? (dEdR/deltaR[ReferenceForce::RIndex]) : zero;
+   dEdR                       = deltaR[ReferenceForce::RIndex] > 0 ? (dEdR/deltaR[ReferenceForce::RIndex]) : 0;
 
    forces[atomAIndex][0]     += dEdR*deltaR[ReferenceForce::XIndex];
    forces[atomAIndex][1]     += dEdR*deltaR[ReferenceForce::YIndex];
@@ -126,6 +114,8 @@ void ReferenceCustomBondIxn::calculateBondIxn(int* atomIndices,
    forces[atomBIndex][1]     -= dEdR*deltaR[ReferenceForce::YIndex];
    forces[atomBIndex][2]     -= dEdR*deltaR[ReferenceForce::ZIndex];
 
+   for (int i = 0; i < energyParamDerivExpressions.size(); i++)
+       energyParamDerivs[i] += energyParamDerivExpressions[i].evaluate();
    if (totalEnergy != NULL)
        *totalEnergy += (RealOpenMM) energyExpression.evaluate();
 }

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 {

platforms/reference/src/SimTKReference/ReferenceCustomDynamics.cpp

@@ -36,6 +36,28 @@
 
 using namespace std;
 using namespace OpenMM;
+using namespace Lepton;
+
+class ReferenceCustomDynamics::DerivFunction : public CustomFunction {
+public:
+    DerivFunction(map<string, double>& energyParamDerivs, const string& param) : energyParamDerivs(energyParamDerivs), param(param) {
+    }
+    int getNumArguments() const {
+        return 0;
+    }
+    double evaluate(const double* arguments) const {
+        return energyParamDerivs[param];
+    }
+    double evaluateDerivative(const double* arguments, const int* derivOrder) const {
+        return 0;
+    }
+    CustomFunction* clone() const {
+        return new DerivFunction(energyParamDerivs, param);
+    }
+private:
+    map<string, double>& energyParamDerivs;
+    string param;
+};
 
 /**---------------------------------------------------------------------------------------
 
@@ -56,7 +78,7 @@ ReferenceCustomDynamics::ReferenceCustomDynamics(int numberOfAtoms, const Custom
         string expression;
         integrator.getComputationStep(i, stepType[i], stepVariable[i], expression);
     }
-    kineticEnergyExpression = Lepton::Parser::parse(integrator.getKineticEnergyExpression()).optimize().createCompiledExpression();
+    kineticEnergyExpression = Parser::parse(integrator.getKineticEnergyExpression()).optimize().createCompiledExpression();
     expressionSet.registerExpression(kineticEnergyExpression);
     kineticEnergyNeedsForce = false;
     if (kineticEnergyExpression.getVariables().find("f") != kineticEnergyExpression.getVariables().end())
@@ -78,13 +100,13 @@ void ReferenceCustomDynamics::initialize(ContextImpl& context, vector<RealOpenMM
 
     int numSteps = stepType.size();
     vector<int> forceGroup;
-    vector<vector<Lepton::ParsedExpression> > expressions;
+    vector<vector<ParsedExpression> > expressions;
     CustomIntegratorUtilities::analyzeComputations(context, integrator, expressions, comparisons, blockEnd, invalidatesForces, needsForces, needsEnergy, computeBothForceAndEnergy, forceGroup);
     stepExpressions.resize(expressions.size());
     for (int i = 0; i < numSteps; i++) {
         stepExpressions[i].resize(expressions[i].size());
         for (int j = 0; j < (int) expressions[i].size(); j++) {
-            stepExpressions[i][j] = expressions[i][j].createCompiledExpression();
+            stepExpressions[i][j] = ParsedExpression(replaceDerivFunctions(expressions[i][j].getRootNode(), context)).createCompiledExpression();
             expressionSet.registerExpression(stepExpressions[i][j]);
         }
         if (stepType[i] == CustomIntegrator::WhileBlockStart)
@@ -141,6 +163,22 @@ void ReferenceCustomDynamics::initialize(ContextImpl& context, vector<RealOpenMM
         stepVariableIndex.push_back(expressionSet.getVariableIndex(stepVariable[i]));
 }
 
+ExpressionTreeNode ReferenceCustomDynamics::replaceDerivFunctions(const ExpressionTreeNode& node, ContextImpl& context) {
+    const Operation& op = node.getOperation();
+    if (op.getId() == Operation::CUSTOM && op.getName() == "deriv") {
+        string param = node.getChildren()[1].getOperation().getName();
+        if (context.getParameters().find(param) == context.getParameters().end())
+            throw OpenMMException("The second argument to deriv() must be a context parameter");
+        return ExpressionTreeNode(new Operation::Custom("deriv", new DerivFunction(energyParamDerivs, param)));
+    }
+    else {
+        vector<ExpressionTreeNode> children;
+        for (int i = 0; i < (int) node.getChildren().size(); i++)
+            children.push_back(replaceDerivFunctions(node.getChildren()[i], context));
+        return ExpressionTreeNode(op.clone(), children);
+    }
+}
+
 /**---------------------------------------------------------------------------------------
 
    Update -- driver routine for performing Custom dynamics update of coordinates
@@ -178,8 +216,10 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
             bool computeEnergy = needsEnergy[step] || computeBothForceAndEnergy[step];
             recordChangedParameters(context, globals);
             RealOpenMM e = context.calcForcesAndEnergy(computeForce, computeEnergy, forceGroupFlags[step]);
-            if (computeEnergy)
+            if (computeEnergy) {
                 energy = e;
+                context.getEnergyParameterDerivatives(energyParamDerivs);
+            }
             forcesAreValid = true;
         }
         expressionSet.setVariable(energyVariableIndex[step], energy);
@@ -266,7 +306,7 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
 
 void ReferenceCustomDynamics::computePerDof(int numberOfAtoms, vector<RealVec>& results, const vector<RealVec>& atomCoordinates,
               const vector<RealVec>& velocities, const vector<RealVec>& forces, const vector<RealOpenMM>& masses,
-              const vector<vector<RealVec> >& perDof, const Lepton::CompiledExpression& expression, int forceIndex) {
+              const vector<vector<RealVec> >& perDof, const CompiledExpression& expression, int forceIndex) {
     // Loop over all degrees of freedom.
 
     for (int i = 0; i < numberOfAtoms; i++) {

platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.cpp

 
-/* Portions copyright (c) 2009 Stanford University and Simbios.
+/* Portions copyright (c) 2009-2016 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -42,38 +42,62 @@ using namespace OpenMM;
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceCustomGBIxn::ReferenceCustomGBIxn(const vector<Lepton::ExpressionProgram>& valueExpressions,
-                     const vector<vector<Lepton::ExpressionProgram> > valueDerivExpressions,
-                     const vector<vector<Lepton::ExpressionProgram> > valueGradientExpressions,
+ReferenceCustomGBIxn::ReferenceCustomGBIxn(const vector<Lepton::CompiledExpression>& valueExpressions,
+                     const vector<vector<Lepton::CompiledExpression> > valueDerivExpressions,
+                     const vector<vector<Lepton::CompiledExpression> > valueGradientExpressions,
+                     const vector<vector<Lepton::CompiledExpression> > valueParamDerivExpressions,
                      const vector<string>& valueNames,
                      const vector<OpenMM::CustomGBForce::ComputationType>& valueTypes,
-                     const vector<Lepton::ExpressionProgram>& energyExpressions,
-                     const vector<vector<Lepton::ExpressionProgram> > energyDerivExpressions,
-                     const vector<vector<Lepton::ExpressionProgram> > energyGradientExpressions,
+                     const vector<Lepton::CompiledExpression>& energyExpressions,
+                     const vector<vector<Lepton::CompiledExpression> > energyDerivExpressions,
+                     const vector<vector<Lepton::CompiledExpression> > energyGradientExpressions,
+                     const vector<vector<Lepton::CompiledExpression> > energyParamDerivExpressions,
                      const vector<OpenMM::CustomGBForce::ComputationType>& energyTypes,
                      const vector<string>& parameterNames) :
-            cutoff(false), periodic(false), valueExpressions(valueExpressions), valueDerivExpressions(valueDerivExpressions), valueGradientExpressions(valueGradientExpressions),
-            valueNames(valueNames), valueTypes(valueTypes), energyExpressions(energyExpressions), energyDerivExpressions(energyDerivExpressions), energyGradientExpressions(energyGradientExpressions),
-            energyTypes(energyTypes), paramNames(parameterNames) {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceCustomGBIxn::ReferenceCustomGBIxn";
-
-   // ---------------------------------------------------------------------------------------
-
-    for (int i = 0; i < (int) paramNames.size(); i++) {
+            cutoff(false), periodic(false), valueExpressions(valueExpressions), valueDerivExpressions(valueDerivExpressions), valueGradientExpressions(valueGradientExpressions), valueParamDerivExpressions(valueParamDerivExpressions),
+            valueTypes(valueTypes), energyExpressions(energyExpressions), energyDerivExpressions(energyDerivExpressions), energyGradientExpressions(energyGradientExpressions), energyParamDerivExpressions(energyParamDerivExpressions),
+            energyTypes(energyTypes) {
+
+    for (int i = 0; i < this->valueExpressions.size(); i++)
+        expressionSet.registerExpression(this->valueExpressions[i]);
+    for (int i = 0; i < this->valueDerivExpressions.size(); i++)
+        for (int j = 0; j < this->valueDerivExpressions[i].size(); j++)
+            expressionSet.registerExpression(this->valueDerivExpressions[i][j]);
+    for (int i = 0; i < this->valueGradientExpressions.size(); i++)
+        for (int j = 0; j < this->valueGradientExpressions[i].size(); j++)
+            expressionSet.registerExpression(this->valueGradientExpressions[i][j]);
+    for (int i = 0; i < this->valueParamDerivExpressions.size(); i++)
+        for (int j = 0; j < this->valueParamDerivExpressions[i].size(); j++)
+            expressionSet.registerExpression(this->valueParamDerivExpressions[i][j]);
+    for (int i = 0; i < this->energyExpressions.size(); i++)
+        expressionSet.registerExpression(this->energyExpressions[i]);
+    for (int i = 0; i < this->energyDerivExpressions.size(); i++)
+        for (int j = 0; j < this->energyDerivExpressions[i].size(); j++)
+            expressionSet.registerExpression(this->energyDerivExpressions[i][j]);
+    for (int i = 0; i < this->energyGradientExpressions.size(); i++)
+        for (int j = 0; j < this->energyGradientExpressions[i].size(); j++)
+            expressionSet.registerExpression(this->energyGradientExpressions[i][j]);
+    for (int i = 0; i < this->energyParamDerivExpressions.size(); i++)
+        for (int j = 0; j < this->energyParamDerivExpressions[i].size(); j++)
+            expressionSet.registerExpression(this->energyParamDerivExpressions[i][j]);
+    rIndex = expressionSet.getVariableIndex("r");
+    xIndex = expressionSet.getVariableIndex("x");
+    yIndex = expressionSet.getVariableIndex("y");
+    zIndex = expressionSet.getVariableIndex("z");
+    for (int i = 0; i < (int) parameterNames.size(); i++) {
+        paramIndex.push_back(expressionSet.getVariableIndex(parameterNames[i]));
         for (int j = 1; j < 3; j++) {
             stringstream name;
-            name << paramNames[i] << j;
-            particleParamNames.push_back(name.str());
+            name << parameterNames[i] << j;
+            particleParamIndex.push_back(expressionSet.getVariableIndex(name.str()));
         }
     }
     for (int i = 0; i < (int) valueNames.size(); i++) {
+        valueIndex.push_back(expressionSet.getVariableIndex(valueNames[i]));
         for (int j = 1; j < 3; j++) {
             stringstream name;
             name << valueNames[i] << j;
-            particleValueNames.push_back(name.str());
+            particleValueIndex.push_back(expressionSet.getVariableIndex(name.str()));
         }
     }
 }
@@ -85,13 +109,6 @@ ReferenceCustomGBIxn::ReferenceCustomGBIxn(const vector<Lepton::ExpressionProgra
    --------------------------------------------------------------------------------------- */
 
 ReferenceCustomGBIxn::~ReferenceCustomGBIxn() {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceCustomGBIxn::~ReferenceCustomGBIxn";
-
-   // ---------------------------------------------------------------------------------------
-
 }
 
   /**---------------------------------------------------------------------------------------
@@ -135,55 +152,73 @@ ReferenceCustomGBIxn::~ReferenceCustomGBIxn() {
 
 void ReferenceCustomGBIxn::calculateIxn(int numberOfAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
                                            const vector<set<int> >& exclusions, map<string, double>& globalParameters, vector<RealVec>& forces,
-                                           RealOpenMM* totalEnergy) const {
+                                           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);
+    
+    // Initialize arrays for storing values.
+    
+    int numValues = valueTypes.size();
+    int numDerivs = valueParamDerivExpressions[0].size();
+    values.resize(numValues);
+    dEdV.resize(numValues, vector<RealOpenMM>(numberOfAtoms, 0.0));
+    dValuedParam.resize(numValues);
+    for (int i = 0; i < numValues; i++)
+        dValuedParam[i].resize(numDerivs, vector<RealOpenMM>(numberOfAtoms, 0.0));
+
     // First calculate the computed values.
 
-    int numValues = valueTypes.size();
-    vector<vector<RealOpenMM> > values(numValues);
     for (int valueIndex = 0; valueIndex < numValues; valueIndex++) {
         if (valueTypes[valueIndex] == OpenMM::CustomGBForce::SingleParticle)
-            calculateSingleParticleValue(valueIndex, numberOfAtoms, atomCoordinates, values, globalParameters, atomParameters);
+            calculateSingleParticleValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters);
         else if (valueTypes[valueIndex] == OpenMM::CustomGBForce::ParticlePair)
-            calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true);
+            calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, exclusions, true);
         else
-            calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, false);
+            calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, exclusions, false);
     }
 
     // Now calculate the energy and its derivates.
 
-    vector<vector<RealOpenMM> > dEdV(numValues, vector<RealOpenMM>(numberOfAtoms, (RealOpenMM) 0));
     for (int termIndex = 0; termIndex < (int) energyExpressions.size(); termIndex++) {
         if (energyTypes[termIndex] == OpenMM::CustomGBForce::SingleParticle)
-            calculateSingleParticleEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, values, globalParameters, atomParameters, forces, totalEnergy, dEdV);
+            calculateSingleParticleEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, forces, totalEnergy, energyParamDerivs);
         else if (energyTypes[termIndex] == OpenMM::CustomGBForce::ParticlePair)
-            calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true, forces, totalEnergy, dEdV);
+            calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, exclusions, true, forces, totalEnergy, energyParamDerivs);
         else
-            calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, false, forces, totalEnergy, dEdV);
+            calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, exclusions, false, forces, totalEnergy, energyParamDerivs);
     }
 
     // Apply the chain rule to evaluate forces.
 
-    calculateChainRuleForces(numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, forces, dEdV);
+    calculateChainRuleForces(numberOfAtoms, atomCoordinates, atomParameters, exclusions, forces, energyParamDerivs);
 }
 
-void ReferenceCustomGBIxn::calculateSingleParticleValue(int index, int numAtoms, vector<RealVec>& atomCoordinates, vector<vector<RealOpenMM> >& values,
-        const map<string, double>& globalParameters, RealOpenMM** atomParameters) const {
+void ReferenceCustomGBIxn::calculateSingleParticleValue(int index, int numAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters) {
     values[index].resize(numAtoms);
-    map<string, double> variables = globalParameters;
     for (int i = 0; i < numAtoms; i++) {
-        variables["x"] = atomCoordinates[i][0];
-        variables["y"] = atomCoordinates[i][1];
-        variables["z"] = atomCoordinates[i][2];
-        for (int j = 0; j < (int) paramNames.size(); j++)
-            variables[paramNames[j]] = atomParameters[i][j];
+        expressionSet.setVariable(xIndex, atomCoordinates[i][0]);
+        expressionSet.setVariable(yIndex, atomCoordinates[i][1]);
+        expressionSet.setVariable(zIndex, atomCoordinates[i][2]);
+        for (int j = 0; j < (int) paramIndex.size(); j++)
+            expressionSet.setVariable(paramIndex[j], atomParameters[i][j]);
         for (int j = 0; j < index; j++)
-            variables[valueNames[j]] = values[j][i];
-        values[index][i] = (RealOpenMM) valueExpressions[index].evaluate(variables);
+            expressionSet.setVariable(valueIndex[j], values[j][i]);
+        values[index][i] = (RealOpenMM) valueExpressions[index].evaluate();
+
+        // Calculate derivatives with respect to parameters.
+
+        for (int j = 0; j < valueParamDerivExpressions[index].size(); j++)
+            dValuedParam[index][j][i] += valueParamDerivExpressions[index][j].evaluate();
+        for (int j = 0; j < index; j++) {
+            RealOpenMM dVdV = valueDerivExpressions[index][j].evaluate();
+            for (int k = 0; k < valueParamDerivExpressions[index].size(); k++)
+                dValuedParam[index][k][i] += dVdV*dValuedParam[j][k][i];
+        }
     }
 }
 
 void ReferenceCustomGBIxn::calculateParticlePairValue(int index, int numAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
-        vector<vector<RealOpenMM> >& values, const map<string, double>& globalParameters, const vector<set<int> >& exclusions, bool useExclusions) const {
+        const vector<set<int> >& exclusions, bool useExclusions) {
     values[index].resize(numAtoms);
     for (int i = 0; i < numAtoms; i++)
         values[index][i] = (RealOpenMM) 0.0;
@@ -194,8 +229,8 @@ void ReferenceCustomGBIxn::calculateParticlePairValue(int index, int numAtoms, v
             OpenMM::AtomPair pair = (*neighborList)[i];
             if (useExclusions && exclusions[pair.first].find(pair.second) != exclusions[pair.first].end())
                 continue;
-            calculateOnePairValue(index, pair.first, pair.second, atomCoordinates, atomParameters, globalParameters, values);
-            calculateOnePairValue(index, pair.second, pair.first, atomCoordinates, atomParameters, globalParameters, values);
+            calculateOnePairValue(index, pair.first, pair.second, atomCoordinates, atomParameters);
+            calculateOnePairValue(index, pair.second, pair.first, atomCoordinates, atomParameters);
         }
     }
     else {
@@ -205,15 +240,14 @@ void ReferenceCustomGBIxn::calculateParticlePairValue(int index, int numAtoms, v
             for (int j = i+1; j < numAtoms; j++) {
                 if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
                     continue;
-                calculateOnePairValue(index, i, j, atomCoordinates, atomParameters, globalParameters, values);
-                calculateOnePairValue(index, j, i, atomCoordinates, atomParameters, globalParameters, values);
+                calculateOnePairValue(index, i, j, atomCoordinates, atomParameters);
+                calculateOnePairValue(index, j, i, atomCoordinates, atomParameters);
            }
         }
     }
 }
 
-void ReferenceCustomGBIxn::calculateOnePairValue(int index, int atom1, int atom2, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
-        const map<string, double>& globalParameters, vector<vector<RealOpenMM> >& values) const {
+void ReferenceCustomGBIxn::calculateOnePairValue(int index, int atom1, int atom2, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters) {
     RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
     if (periodic)
         ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxVectors, deltaR);
@@ -222,44 +256,53 @@ void ReferenceCustomGBIxn::calculateOnePairValue(int index, int atom1, int atom2
     RealOpenMM r = deltaR[ReferenceForce::RIndex];
     if (cutoff && r >= cutoffDistance)
         return;
-    map<string, double> variables = globalParameters;
-    for (int i = 0; i < (int) paramNames.size(); i++) {
-        variables[particleParamNames[i*2]] = atomParameters[atom1][i];
-        variables[particleParamNames[i*2+1]] = atomParameters[atom2][i];
+    for (int i = 0; i < (int) paramIndex.size(); i++) {
+        expressionSet.setVariable(particleParamIndex[i*2], atomParameters[atom1][i]);
+        expressionSet.setVariable(particleParamIndex[i*2+1], atomParameters[atom2][i]);
     }
-    variables["r"] = r;
+    expressionSet.setVariable(rIndex, r);
     for (int i = 0; i < index; i++) {
-        variables[particleValueNames[i*2]] = values[i][atom1];
-        variables[particleValueNames[i*2+1]] = values[i][atom2];
+        expressionSet.setVariable(particleValueIndex[i*2], values[i][atom1]);
+        expressionSet.setVariable(particleValueIndex[i*2+1], values[i][atom2]);
     }
-    values[index][atom1] += (RealOpenMM) valueExpressions[index].evaluate(variables);
+    values[index][atom1] += (RealOpenMM) valueExpressions[index].evaluate();
+    
+    // Calculate derivatives with respect to parameters.
+    
+    for (int i = 0; i < valueParamDerivExpressions[index].size(); i++)
+        dValuedParam[index][i][atom1] += valueParamDerivExpressions[index][i].evaluate();
 }
 
-void ReferenceCustomGBIxn::calculateSingleParticleEnergyTerm(int index, int numAtoms, vector<RealVec>& atomCoordinates, const vector<vector<RealOpenMM> >& values,
-        const map<string, double>& globalParameters, RealOpenMM** atomParameters, vector<RealVec>& forces, RealOpenMM* totalEnergy,
-        vector<vector<RealOpenMM> >& dEdV) const {
-    map<string, double> variables = globalParameters;
+void ReferenceCustomGBIxn::calculateSingleParticleEnergyTerm(int index, int numAtoms, vector<RealVec>& atomCoordinates,
+        RealOpenMM** atomParameters, vector<RealVec>& forces, RealOpenMM* totalEnergy, double* energyParamDerivs) {
     for (int i = 0; i < numAtoms; i++) {
-        variables["x"] = atomCoordinates[i][0];
-        variables["y"] = atomCoordinates[i][1];
-        variables["z"] = atomCoordinates[i][2];
-        for (int j = 0; j < (int) paramNames.size(); j++)
-            variables[paramNames[j]] = atomParameters[i][j];
-        for (int j = 0; j < (int) valueNames.size(); j++)
-            variables[valueNames[j]] = values[j][i];
+        expressionSet.setVariable(xIndex, atomCoordinates[i][0]);
+        expressionSet.setVariable(yIndex, atomCoordinates[i][1]);
+        expressionSet.setVariable(zIndex, atomCoordinates[i][2]);
+        for (int j = 0; j < (int) paramIndex.size(); j++)
+            expressionSet.setVariable(paramIndex[j], atomParameters[i][j]);
+        for (int j = 0; j < valueIndex.size(); j++)
+            expressionSet.setVariable(valueIndex[j], values[j][i]);
+        
+        // Compute energy and force.
+        
         if (totalEnergy != NULL)
-            *totalEnergy += (RealOpenMM) energyExpressions[index].evaluate(variables);
-        for (int j = 0; j < (int) valueNames.size(); j++)
-            dEdV[j][i] += (RealOpenMM) energyDerivExpressions[index][j].evaluate(variables);
-        forces[i][0] -= (RealOpenMM) energyGradientExpressions[index][0].evaluate(variables);
-        forces[i][1] -= (RealOpenMM) energyGradientExpressions[index][1].evaluate(variables);
-        forces[i][2] -= (RealOpenMM) energyGradientExpressions[index][2].evaluate(variables);
+            *totalEnergy += (RealOpenMM) energyExpressions[index].evaluate();
+        for (int j = 0; j < (int) valueIndex.size(); j++)
+            dEdV[j][i] += (RealOpenMM) energyDerivExpressions[index][j].evaluate();
+        forces[i][0] -= (RealOpenMM) energyGradientExpressions[index][0].evaluate();
+        forces[i][1] -= (RealOpenMM) energyGradientExpressions[index][1].evaluate();
+        forces[i][2] -= (RealOpenMM) energyGradientExpressions[index][2].evaluate();
+        
+        // Compute derivatives with respect to parameters.
+        
+        for (int k = 0; k < energyParamDerivExpressions[index].size(); k++)
+            energyParamDerivs[k] += energyParamDerivExpressions[index][k].evaluate();
     }
 }
 
 void ReferenceCustomGBIxn::calculateParticlePairEnergyTerm(int index, int numAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
-        const vector<vector<RealOpenMM> >& values, const map<string, double>& globalParameters, const vector<set<int> >& exclusions, bool useExclusions,
-        vector<RealVec>& forces, RealOpenMM* totalEnergy, vector<vector<RealOpenMM> >& dEdV) const {
+        const vector<set<int> >& exclusions, bool useExclusions, vector<RealVec>& forces, RealOpenMM* totalEnergy, double* energyParamDerivs) {
     if (cutoff) {
         // Loop over all pairs in the neighbor list.
 
@@ -267,7 +310,7 @@ void ReferenceCustomGBIxn::calculateParticlePairEnergyTerm(int index, int numAto
             OpenMM::AtomPair pair = (*neighborList)[i];
             if (useExclusions && exclusions[pair.first].find(pair.second) != exclusions[pair.first].end())
                 continue;
-            calculateOnePairEnergyTerm(index, pair.first, pair.second, atomCoordinates, atomParameters, globalParameters, values, forces, totalEnergy, dEdV);
+            calculateOnePairEnergyTerm(index, pair.first, pair.second, atomCoordinates, atomParameters, forces, totalEnergy, energyParamDerivs);
         }
     }
     else {
@@ -277,15 +320,14 @@ void ReferenceCustomGBIxn::calculateParticlePairEnergyTerm(int index, int numAto
             for (int j = i+1; j < numAtoms; j++) {
                 if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
                     continue;
-                calculateOnePairEnergyTerm(index, i, j, atomCoordinates, atomParameters, globalParameters, values, forces, totalEnergy, dEdV);
+                calculateOnePairEnergyTerm(index, i, j, atomCoordinates, atomParameters, forces, totalEnergy, energyParamDerivs);
            }
         }
     }
 }
 
 void ReferenceCustomGBIxn::calculateOnePairEnergyTerm(int index, int atom1, int atom2, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
-        const map<string, double>& globalParameters, const vector<vector<RealOpenMM> >& values, vector<RealVec>& forces, RealOpenMM* totalEnergy,
-        vector<vector<RealOpenMM> >& dEdV) const {
+        vector<RealVec>& forces, RealOpenMM* totalEnergy, double* energyParamDerivs) {
     // Compute the displacement.
 
     RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
@@ -299,44 +341,47 @@ void ReferenceCustomGBIxn::calculateOnePairEnergyTerm(int index, int atom1, int
 
     // Record variables for evaluating expressions.
 
-    map<string, double> variables = globalParameters;
-    for (int i = 0; i < (int) paramNames.size(); i++) {
-        variables[particleParamNames[i*2]] = atomParameters[atom1][i];
-        variables[particleParamNames[i*2+1]] = atomParameters[atom2][i];
+    for (int i = 0; i < (int) paramIndex.size(); i++) {
+        expressionSet.setVariable(particleParamIndex[i*2], atomParameters[atom1][i]);
+        expressionSet.setVariable(particleParamIndex[i*2+1], atomParameters[atom2][i]);
     }
-    variables["r"] = r;
-    for (int i = 0; i < (int) valueNames.size(); i++) {
-        variables[particleValueNames[i*2]] = values[i][atom1];
-        variables[particleValueNames[i*2+1]] = values[i][atom2];
+    expressionSet.setVariable(rIndex, r);
+    for (int i = 0; i < (int) valueIndex.size(); i++) {
+        expressionSet.setVariable(particleValueIndex[i*2], values[i][atom1]);
+        expressionSet.setVariable(particleValueIndex[i*2+1], values[i][atom2]);
     }
 
     // Evaluate the energy and its derivatives.
 
     if (totalEnergy != NULL)
-        *totalEnergy += (RealOpenMM) energyExpressions[index].evaluate(variables);
-    RealOpenMM dEdR = (RealOpenMM) energyDerivExpressions[index][0].evaluate(variables);
+        *totalEnergy += (RealOpenMM) energyExpressions[index].evaluate();
+    RealOpenMM dEdR = (RealOpenMM) energyDerivExpressions[index][0].evaluate();
     dEdR *= 1/r;
     for (int i = 0; i < 3; i++) {
        forces[atom1][i] -= dEdR*deltaR[i];
        forces[atom2][i] += dEdR*deltaR[i];
     }
-    for (int i = 0; i < (int) valueNames.size(); i++) {
-        dEdV[i][atom1] += (RealOpenMM) energyDerivExpressions[index][2*i+1].evaluate(variables);
-        dEdV[i][atom2] += (RealOpenMM) energyDerivExpressions[index][2*i+2].evaluate(variables);
+    for (int i = 0; i < (int) valueIndex.size(); i++) {
+        dEdV[i][atom1] += (RealOpenMM) energyDerivExpressions[index][2*i+1].evaluate();
+        dEdV[i][atom2] += (RealOpenMM) energyDerivExpressions[index][2*i+2].evaluate();
     }
+        
+    // Compute derivatives with respect to parameters.
+
+    for (int i = 0; i < energyParamDerivExpressions[index].size(); i++)
+        energyParamDerivs[i] += energyParamDerivExpressions[index][i].evaluate();
 }
 
 void ReferenceCustomGBIxn::calculateChainRuleForces(int numAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
-        const vector<vector<RealOpenMM> >& values, const map<string, double>& globalParameters,
-        const vector<set<int> >& exclusions, vector<RealVec>& forces, vector<vector<RealOpenMM> >& dEdV) const {
+        const vector<set<int> >& exclusions, vector<RealVec>& forces, double* energyParamDerivs) {
     if (cutoff) {
         // Loop over all pairs in the neighbor list.
 
         for (int i = 0; i < (int) neighborList->size(); i++) {
             OpenMM::AtomPair pair = (*neighborList)[i];
             bool isExcluded = (exclusions[pair.first].find(pair.second) != exclusions[pair.first].end());
-            calculateOnePairChainRule(pair.first, pair.second, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
-            calculateOnePairChainRule(pair.second, pair.first, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
+            calculateOnePairChainRule(pair.first, pair.second, atomCoordinates, atomParameters, forces, isExcluded);
+            calculateOnePairChainRule(pair.second, pair.first, atomCoordinates, atomParameters, forces, isExcluded);
         }
     }
     else {
@@ -345,45 +390,50 @@ void ReferenceCustomGBIxn::calculateChainRuleForces(int numAtoms, vector<RealVec
         for (int i = 0; i < numAtoms; i++) {
             for (int j = i+1; j < numAtoms; j++) {
                 bool isExcluded = (exclusions[i].find(j) != exclusions[i].end());
-                calculateOnePairChainRule(i, j, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
-                calculateOnePairChainRule(j, i, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
+                calculateOnePairChainRule(i, j, atomCoordinates, atomParameters, forces, isExcluded);
+                calculateOnePairChainRule(j, i, atomCoordinates, atomParameters, forces, isExcluded);
            }
         }
     }
 
     // Compute chain rule terms for computed values that depend explicitly on particle coordinates.
 
-    map<string, double> variables = globalParameters;
     for (int i = 0; i < numAtoms; i++) {
-        variables["x"] = atomCoordinates[i][0];
-        variables["y"] = atomCoordinates[i][1];
-        variables["z"] = atomCoordinates[i][2];
+        expressionSet.setVariable(xIndex, atomCoordinates[i][0]);
+        expressionSet.setVariable(yIndex, atomCoordinates[i][1]);
+        expressionSet.setVariable(zIndex, atomCoordinates[i][2]);
         vector<RealOpenMM> dVdX(valueDerivExpressions.size(), 0.0);
         vector<RealOpenMM> dVdY(valueDerivExpressions.size(), 0.0);
         vector<RealOpenMM> dVdZ(valueDerivExpressions.size(), 0.0);
-        for (int j = 0; j < (int) paramNames.size(); j++)
-            variables[paramNames[j]] = atomParameters[i][j];
-        for (int j = 1; j < (int) valueNames.size(); j++) {
-            variables[valueNames[j-1]] = values[j-1][i];
+        for (int j = 0; j < (int) paramIndex.size(); j++)
+            expressionSet.setVariable(paramIndex[j], atomParameters[i][j]);
+        for (int j = 1; j < (int) valueIndex.size(); j++) {
+            expressionSet.setVariable(valueIndex[j-1], values[j-1][i]);
             for (int k = 1; k < j; k++) {
-                RealOpenMM dVdV = (RealOpenMM) valueDerivExpressions[j][k].evaluate(variables);
+                RealOpenMM dVdV = (RealOpenMM) valueDerivExpressions[j][k].evaluate();
                 dVdX[j] += dVdV*dVdX[k];
                 dVdY[j] += dVdV*dVdY[k];
                 dVdZ[j] += dVdV*dVdZ[k];
             }
-            dVdX[j] += (RealOpenMM) valueGradientExpressions[j][0].evaluate(variables);
-            dVdY[j] += (RealOpenMM) valueGradientExpressions[j][1].evaluate(variables);
-            dVdZ[j] += (RealOpenMM) valueGradientExpressions[j][2].evaluate(variables);
+            dVdX[j] += (RealOpenMM) valueGradientExpressions[j][0].evaluate();
+            dVdY[j] += (RealOpenMM) valueGradientExpressions[j][1].evaluate();
+            dVdZ[j] += (RealOpenMM) valueGradientExpressions[j][2].evaluate();
             forces[i][0] -= dEdV[j][i]*dVdX[j];
             forces[i][1] -= dEdV[j][i]*dVdY[j];
             forces[i][2] -= dEdV[j][i]*dVdZ[j];
         }
     }
+        
+    // Compute chain rule terms for derivatives with respect to parameters.
+
+    for (int i = 0; i < numAtoms; i++)
+        for (int j = 0; j < (int) valueIndex.size(); j++)
+            for (int k = 0; k < dValuedParam[j].size(); k++)
+                energyParamDerivs[k] += dEdV[j][i]*dValuedParam[j][k][i];
 }
 
 void ReferenceCustomGBIxn::calculateOnePairChainRule(int atom1, int atom2, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
-        const map<string, double>& globalParameters, const vector<vector<RealOpenMM> >& values, vector<RealVec>& forces,
-        vector<vector<RealOpenMM> >& dEdV, bool isExcluded) const {
+        vector<RealVec>& forces, bool isExcluded) {
     // Compute the displacement.
 
     RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
@@ -397,14 +447,13 @@ void ReferenceCustomGBIxn::calculateOnePairChainRule(int atom1, int atom2, vecto
 
     // Record variables for evaluating expressions.
 
-    map<string, double> variables = globalParameters;
-    for (int i = 0; i < (int) paramNames.size(); i++) {
-        variables[particleParamNames[i*2]] = atomParameters[atom1][i];
-        variables[particleParamNames[i*2+1]] = atomParameters[atom2][i];
+    for (int i = 0; i < (int) paramIndex.size(); i++) {
+        expressionSet.setVariable(particleParamIndex[i*2], atomParameters[atom1][i]);
+        expressionSet.setVariable(particleParamIndex[i*2+1], atomParameters[atom2][i]);
     }
-    variables["r"] = r;
-    variables[particleValueNames[0]] = values[0][atom1];
-    variables[particleValueNames[1]] = values[0][atom2];
+    expressionSet.setVariable(rIndex, r);
+    expressionSet.setVariable(particleValueIndex[0], values[0][atom1]);
+    expressionSet.setVariable(particleValueIndex[1], values[0][atom2]);
 
     // Evaluate the derivative of each parameter with respect to position and apply forces.
 
@@ -415,24 +464,23 @@ void ReferenceCustomGBIxn::calculateOnePairChainRule(int atom1, int atom2, vecto
     vector<RealOpenMM> dVdR1(valueDerivExpressions.size(), 0.0);
     vector<RealOpenMM> dVdR2(valueDerivExpressions.size(), 0.0);
     if (!isExcluded || valueTypes[0] != OpenMM::CustomGBForce::ParticlePair) {
-        dVdR1[0] = (RealOpenMM) valueDerivExpressions[0][0].evaluate(variables);;
+        dVdR1[0] = (RealOpenMM) valueDerivExpressions[0][0].evaluate();
         dVdR2[0] = -dVdR1[0];
         for (int i = 0; i < 3; i++) {
             forces[atom1][i] -= dEdV[0][atom1]*dVdR1[0]*deltaR[i];
             forces[atom2][i] -= dEdV[0][atom1]*dVdR2[0]*deltaR[i];
         }
     }
-    variables = globalParameters;
-    for (int i = 0; i < (int) paramNames.size(); i++)
-        variables[paramNames[i]] = atomParameters[atom1][i];
-    variables[valueNames[0]] = values[0][atom1];
-    for (int i = 1; i < (int) valueNames.size(); i++) {
-        variables[valueNames[i]] = values[i][atom1];
-        variables["x"] = atomCoordinates[atom1][0];
-        variables["y"] = atomCoordinates[atom1][1];
-        variables["z"] = atomCoordinates[atom1][2];
+    for (int i = 0; i < (int) paramIndex.size(); i++)
+        expressionSet.setVariable(paramIndex[i], atomParameters[atom1][i]);
+    expressionSet.setVariable(valueIndex[0], values[0][atom1]);
+    for (int i = 1; i < (int) valueIndex.size(); i++) {
+        expressionSet.setVariable(valueIndex[i], values[i][atom1]);
+        expressionSet.setVariable(xIndex, atomCoordinates[atom1][0]);
+        expressionSet.setVariable(yIndex, atomCoordinates[atom1][1]);
+        expressionSet.setVariable(zIndex, atomCoordinates[atom1][2]);
         for (int j = 0; j < i; j++) {
-            RealOpenMM dVdV = (RealOpenMM) valueDerivExpressions[i][j].evaluate(variables);
+            RealOpenMM dVdV = (RealOpenMM) valueDerivExpressions[i][j].evaluate();
             dVdR1[i] += dVdV*dVdR1[j];
             dVdR2[i] += dVdV*dVdR2[j];
         }