Began creating CPU version of CustomGBForce

platforms/cpu/include/CpuCustomGBForce.h

+
+/* Portions copyright (c) 2009-2014 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef OPENMM_CPU_CUSTOM_GB_FORCE_H__
+#define OPENMM_CPU_CUSTOM_GB_FORCE_H__
+
+#include "ReferenceNeighborList.h"
+#include "CompiledExpressionSet.h"
+#include "lepton/CompiledExpression.h"
+#include "openmm/CustomGBForce.h"
+#include <map>
+#include <set>
+#include <vector>
+
+class CpuCustomGBForce {
+private:
+    bool cutoff;
+    bool periodic;
+    const OpenMM::NeighborList* neighborList;
+    RealOpenMM periodicBoxSize[3];
+    RealOpenMM cutoffDistance;
+    OpenMM::CompiledExpressionSet expressionSet;
+    std::vector<Lepton::CompiledExpression> valueExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > valueDerivExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > valueGradientExpressions;
+    std::vector<std::string> valueNames;
+    std::vector<int> valueIndex;
+    std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
+    std::vector<Lepton::CompiledExpression> energyExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > energyDerivExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > energyGradientExpressions;
+    std::vector<std::string> paramNames;
+    std::vector<int> paramIndex;
+    std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
+    std::vector<std::string> particleParamNames;
+    std::vector<std::string> particleValueNames;
+    std::vector<int> particleParamIndex;
+    std::vector<int> particleValueIndex;
+    int xindex, yindex, zindex, rindex;
+
+    /**
+     * Calculate a computed value of type SingleParticle
+     * 
+     * @param index            the index of the value to compute
+     * @param numAtoms         number of atoms
+     * @param atomCoordinates  atom coordinates
+     * @param values           the vector to store computed values into
+     * @param globalParameters the values of global parameters
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     */
+
+    void calculateSingleParticleValue(int index, int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<std::vector<RealOpenMM> >& values,
+                                      const std::map<std::string, double>& globalParameters, RealOpenMM** atomParameters);
+
+    /**
+     * Calculate a computed value that is based on particle pairs
+     * 
+     * @param index            the index of the value to compute
+     * @param numAtoms         number of atoms
+     * @param atomCoordinates  atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param values           the vector to store computed values into
+     * @param globalParameters the values of global parameters
+     * @param exclusions       exclusions[i] is the set of excluded indices for atom i
+     * @param useExclusions    specifies whether to use exclusions
+     */
+
+    void calculateParticlePairValue(int index, int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
+                                    std::vector<std::vector<RealOpenMM> >& values,
+                                    const std::map<std::string, double>& globalParameters,
+                                    const std::vector<std::set<int> >& exclusions, bool useExclusions);
+
+    /**
+     * Evaluate a single atom pair as part of calculating a computed value
+     * 
+     * @param index            the index of the value to compute
+     * @param atom1            the index of the first atom in the pair
+     * @param atom2            the index of the second atom in the pair
+     * @param atomCoordinates  atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param globalParameters the values of global parameters
+     * @param values           the vector to store computed values into
+     */
+
+    void calculateOnePairValue(int index, int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
+                               const std::map<std::string, double>& globalParameters,
+                               std::vector<std::vector<RealOpenMM> >& values);
+
+    /**
+     * Calculate an energy term of type SingleParticle
+     * 
+     * @param index            the index of the value to compute
+     * @param numAtoms         number of atoms
+     * @param atomCoordinates  atom coordinates
+     * @param values           the vector containing computed values
+     * @param globalParameters the values of global parameters
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param forces           forces on atoms are added to this
+     * @param totalEnergy      the energy contribution is added to this
+     * @param dEdV             the derivative of energy with respect to computed values is stored in this
+     */
+
+    void calculateSingleParticleEnergyTerm(int index, int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, const std::vector<std::vector<RealOpenMM> >& values,
+                                      const std::map<std::string, double>& globalParameters, RealOpenMM** atomParameters, std::vector<OpenMM::RealVec>& forces,
+                                      RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV);
+
+    /**
+     * Calculate an energy term that is based on particle pairs
+     * 
+     * @param index            the index of the term to compute
+     * @param numAtoms         number of atoms
+     * @param atomCoordinates  atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param values           the vector containing computed values
+     * @param globalParameters the values of global parameters
+     * @param exclusions       exclusions[i] is the set of excluded indices for atom i
+     * @param useExclusions    specifies whether to use exclusions
+     * @param forces           forces on atoms are added to this
+     * @param totalEnergy      the energy contribution is added to this
+     * @param dEdV             the derivative of energy with respect to computed values is stored in this
+     */
+
+    void calculateParticlePairEnergyTerm(int index, int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
+                                    const std::vector<std::vector<RealOpenMM> >& values,
+                                    const std::map<std::string, double>& globalParameters,
+                                    const std::vector<std::set<int> >& exclusions, bool useExclusions,
+                                    std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV);
+
+    /**
+     * Evaluate a single atom pair as part of calculating an energy term
+     * 
+     * @param index            the index of the term to compute
+     * @param atom1            the index of the first atom in the pair
+     * @param atom2            the index of the second atom in the pair
+     * @param atomCoordinates  atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param globalParameters the values of global parameters
+     * @param values           the vector containing computed values
+     * @param forces           forces on atoms are added to this
+     * @param totalEnergy      the energy contribution is added to this
+     * @param dEdV             the derivative of energy with respect to computed values is stored in this
+     */
+
+    void calculateOnePairEnergyTerm(int index, int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
+                               const std::map<std::string, double>& globalParameters,
+                               const std::vector<std::vector<RealOpenMM> >& values,
+                               std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV);
+
+    /**
+     * Apply the chain rule to compute forces on atoms
+     * 
+     * @param numAtoms         number of atoms
+     * @param atomCoordinates  atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param values           the vector containing computed values
+     * @param globalParameters the values of global parameters
+     * @param exclusions       exclusions[i] is the set of excluded indices for atom i
+     * @param forces           forces on atoms are added to this
+     * @param dEdV             the derivative of energy with respect to computed values is stored in this
+     */
+
+    void calculateChainRuleForces(int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
+                                    const std::vector<std::vector<RealOpenMM> >& values,
+                                    const std::map<std::string, double>& globalParameters,
+                                    const std::vector<std::set<int> >& exclusions,
+                                    std::vector<OpenMM::RealVec>& forces, std::vector<std::vector<RealOpenMM> >& dEdV);
+
+    /**
+     * Evaluate a single atom pair as part of applying the chain rule
+     * 
+     * @param atom1            the index of the first atom in the pair
+     * @param atom2            the index of the second atom in the pair
+     * @param atomCoordinates  atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param globalParameters the values of global parameters
+     * @param values           the vector containing computed values
+     * @param forces           forces on atoms are added to this
+     * @param dEdV             the derivative of energy with respect to computed values is stored in this
+     * @param isExcluded       specifies whether this is an excluded pair
+     */
+
+    void calculateOnePairChainRule(int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
+                               const std::map<std::string, double>& globalParameters,
+                               const std::vector<std::vector<RealOpenMM> >& values,
+                               std::vector<OpenMM::RealVec>& forces, std::vector<std::vector<RealOpenMM> >& dEdV,
+                               bool isExcluded);
+
+public:
+
+    /**
+     * Construct a new CpuCustomGBForce.
+     */
+
+     CpuCustomGBForce(const std::vector<Lepton::CompiledExpression>& valueExpressions,
+                          const std::vector<std::vector<Lepton::CompiledExpression> > valueDerivExpressions,
+                          const std::vector<std::vector<Lepton::CompiledExpression> > valueGradientExpressions,
+                          const std::vector<std::string>& valueNames,
+                          const std::vector<OpenMM::CustomGBForce::ComputationType>& valueTypes,
+                          const std::vector<Lepton::CompiledExpression>& energyExpressions,
+                          const std::vector<std::vector<Lepton::CompiledExpression> > energyDerivExpressions,
+                          const std::vector<std::vector<Lepton::CompiledExpression> > energyGradientExpressions,
+                          const std::vector<OpenMM::CustomGBForce::ComputationType>& energyTypes,
+                          const std::vector<std::string>& parameterNames);
+
+     ~CpuCustomGBForce();
+
+    /**
+     * Set the force to use a cutoff.
+     * 
+     * @param distance            the cutoff distance
+     * @param neighbors           the neighbor list to use
+     */
+
+    void setUseCutoff(RealOpenMM distance, const OpenMM::NeighborList& neighbors);
+
+    /**
+     * Set the force to use periodic boundary conditions.  This requires that a cutoff has
+     * already been set, and the smallest side of the periodic box is at least twice the cutoff
+     * distance.
+     * 
+     * @param boxSize             the X, Y, and Z widths of the periodic box
+     */
+
+    void setPeriodic(OpenMM::RealVec& boxSize);
+
+    /**
+     * Calculate custom GB ixn
+     * 
+     * @param numberOfAtoms    number of atoms
+     * @param atomCoordinates  atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param exclusions       exclusions[i] is the set of excluded indices for atom i
+     * @param globalParameters the values of global parameters
+     * @param forces           force array (forces added)
+     * @param totalEnergy      total energy
+     */
+
+    void calculateIxn(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters, const std::vector<std::set<int> >& exclusions,
+                     std::map<std::string, double>& globalParameters, std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy);
+};
+
+#endif // OPENMM_CPU_CUSTOM_GB_FORCE_H__

platforms/cpu/include/CpuKernels.h

@@ -33,6 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "CpuBondForce.h"
+#include "CpuCustomGBForce.h"
 #include "CpuCustomManyParticleForce.h"
 #include "CpuCustomNonbondedForce.h"
 #include "CpuGBSAOBCForce.h"
@@ -303,6 +304,58 @@ private:
     CpuGBSAOBCForce obc;
 };
 
+/**
+ * This kernel is invoked by CustomGBForce to calculate the forces acting on the system.
+ */
+class CpuCalcCustomGBForceKernel : public CalcCustomGBForceKernel {
+public:
+    CpuCalcCustomGBForceKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) :
+            CalcCustomGBForceKernel(name, platform), data(data) {
+    }
+    ~CpuCalcCustomGBForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomGBForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomGBForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomGBForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomGBForce& force);
+private:
+    CpuPlatform::PlatformData& data;
+    int numParticles;
+    bool isPeriodic;
+    RealOpenMM **particleParamArray;
+    RealOpenMM nonbondedCutoff;
+    std::vector<std::set<int> > exclusions;
+    std::vector<std::string> particleParameterNames, globalParameterNames, valueNames;
+    std::vector<Lepton::CompiledExpression> valueExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > valueDerivExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > valueGradientExpressions;
+    std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
+    std::vector<Lepton::CompiledExpression> energyExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > energyDerivExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > energyGradientExpressions;
+    std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
+    NonbondedMethod nonbondedMethod;
+    NeighborList* neighborList;
+};
+
 /**
  * This kernel is invoked by CustomManyParticleForce to calculate the forces acting on the system and the energy of the system.
  */

platforms/cpu/src/CpuCustomGBForce.cpp

+
+/* Portions copyright (c) 2009-2014 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include <string.h>
+#include <sstream>
+
+#include "SimTKOpenMMCommon.h"
+#include "SimTKOpenMMLog.h"
+#include "SimTKOpenMMUtilities.h"
+#include "ReferenceForce.h"
+#include "CpuCustomGBForce.h"
+
+using std::map;
+using std::set;
+using std::string;
+using std::stringstream;
+using std::vector;
+using OpenMM::RealVec;
+
+CpuCustomGBForce::CpuCustomGBForce(const vector<Lepton::CompiledExpression>& valueExpressions,
+                     const vector<vector<Lepton::CompiledExpression> > valueDerivExpressions,
+                     const vector<vector<Lepton::CompiledExpression> > valueGradientExpressions,
+                     const vector<string>& valueNames,
+                     const vector<OpenMM::CustomGBForce::ComputationType>& valueTypes,
+                     const vector<Lepton::CompiledExpression>& energyExpressions,
+                     const vector<vector<Lepton::CompiledExpression> > energyDerivExpressions,
+                     const vector<vector<Lepton::CompiledExpression> > energyGradientExpressions,
+                     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) {
+    for (int i = 0; i < (int) valueExpressions.size(); i++)
+        expressionSet.registerExpression(this->valueExpressions[i]);
+    for (int i = 0; i < (int) valueDerivExpressions.size(); i++)
+        for (int j = 0; j < (int) valueDerivExpressions[i].size(); j++)
+            expressionSet.registerExpression(this->valueDerivExpressions[i][j]);
+    for (int i = 0; i < (int) valueGradientExpressions.size(); i++)
+        for (int j = 0; j < (int) valueGradientExpressions[i].size(); j++)
+            expressionSet.registerExpression(this->valueGradientExpressions[i][j]);
+    for (int i = 0; i < (int) energyExpressions.size(); i++)
+        expressionSet.registerExpression(this->energyExpressions[i]);
+    for (int i = 0; i < (int) energyDerivExpressions.size(); i++)
+        for (int j = 0; j < (int) energyDerivExpressions[i].size(); j++)
+            expressionSet.registerExpression(this->energyDerivExpressions[i][j]);
+    for (int i = 0; i < (int) energyGradientExpressions.size(); i++)
+        for (int j = 0; j < (int) energyGradientExpressions[i].size(); j++)
+            expressionSet.registerExpression(this->energyGradientExpressions[i][j]);
+    xindex = expressionSet.getVariableIndex("x");
+    yindex = expressionSet.getVariableIndex("y");
+    zindex = expressionSet.getVariableIndex("z");
+    rindex = expressionSet.getVariableIndex("r");
+    for (int i = 0; i < (int) paramNames.size(); i++) {
+        paramIndex.push_back(expressionSet.getVariableIndex(paramNames[i]));
+        for (int j = 1; j < 3; j++) {
+            stringstream name;
+            name << paramNames[i] << j;
+            particleParamNames.push_back(name.str());
+            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()));
+        }
+    }
+}
+
+CpuCustomGBForce::~CpuCustomGBForce() {
+}
+
+void CpuCustomGBForce::setUseCutoff(RealOpenMM distance, const OpenMM::NeighborList& neighbors) {
+    cutoff = true;
+    cutoffDistance = distance;
+    neighborList = &neighbors;
+}
+
+void CpuCustomGBForce::setPeriodic(RealVec& boxSize) {
+
+    if (cutoff) {
+        assert(boxSize[0] >= 2.0*cutoffDistance);
+        assert(boxSize[1] >= 2.0*cutoffDistance);
+        assert(boxSize[2] >= 2.0*cutoffDistance);
+    }
+    periodic = true;
+    periodicBoxSize[0] = boxSize[0];
+    periodicBoxSize[1] = boxSize[1];
+    periodicBoxSize[2] = boxSize[2];
+  }
+
+void CpuCustomGBForce::calculateIxn(int numberOfAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
+                                           const vector<set<int> >& exclusions, map<string, double>& globalParameters, vector<RealVec>& forces,
+                                           RealOpenMM* totalEnergy) {
+    for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter)
+        expressionSet.setVariable(expressionSet.getVariableIndex(iter->first), iter->second);
+
+    // 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);
+        else if (valueTypes[valueIndex] == OpenMM::CustomGBForce::ParticlePair)
+            calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true);
+        else
+            calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, false);
+    }
+
+    // Now calculate the energy and its derivatives.
+
+    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);
+        else if (energyTypes[termIndex] == OpenMM::CustomGBForce::ParticlePair)
+            calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true, forces, totalEnergy, dEdV);
+        else
+            calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, false, forces, totalEnergy, dEdV);
+    }
+
+    // Apply the chain rule to evaluate forces.
+
+    calculateChainRuleForces(numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, forces, dEdV);
+}
+
+void CpuCustomGBForce::calculateSingleParticleValue(int index, int numAtoms, vector<RealVec>& atomCoordinates, vector<vector<RealOpenMM> >& values,
+        const map<string, double>& globalParameters, RealOpenMM** atomParameters) {
+    values[index].resize(numAtoms);
+    for (int i = 0; i < numAtoms; 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) paramNames.size(); j++)
+            expressionSet.setVariable(paramIndex[j], atomParameters[i][j]);
+        for (int j = 0; j < index; j++)
+            expressionSet.setVariable(valueIndex[j], values[j][i]);
+        values[index][i] = (RealOpenMM) valueExpressions[index].evaluate();
+    }
+}
+
+void CpuCustomGBForce::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) {
+    values[index].resize(numAtoms);
+    for (int i = 0; i < numAtoms; i++)
+        values[index][i] = (RealOpenMM) 0.0;
+    if (cutoff) {
+        // Loop over all pairs in the neighbor list.
+
+        for (int i = 0; i < (int) neighborList->size(); i++) {
+            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);
+        }
+    }
+    else {
+        // Perform an O(N^2) loop over all atom pairs.
+
+        for (int i = 0; i < numAtoms; i++) {
+            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);
+           }
+        }
+    }
+}
+
+void CpuCustomGBForce::calculateOnePairValue(int index, int atom1, int atom2, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
+        const map<string, double>& globalParameters, vector<vector<RealOpenMM> >& values) {
+    RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+    if (periodic)
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxSize, deltaR);
+    else
+        ReferenceForce::getDeltaR(atomCoordinates[atom2], atomCoordinates[atom1], deltaR);
+    RealOpenMM r = deltaR[ReferenceForce::RIndex];
+    if (cutoff && r >= cutoffDistance)
+        return;
+    for (int i = 0; i < (int) paramNames.size(); i++) {
+        expressionSet.setVariable(particleParamIndex[i*2], atomParameters[atom1][i]);
+        expressionSet.setVariable(particleParamIndex[i*2+1], atomParameters[atom2][i]);
+    }
+    expressionSet.setVariable(rindex, r);
+    for (int i = 0; i < index; i++) {
+        expressionSet.setVariable(particleValueIndex[i*2], values[i][atom1]);
+        expressionSet.setVariable(particleValueIndex[i*2+1], values[i][atom2]);
+    }
+    values[index][atom1] += (RealOpenMM) valueExpressions[index].evaluate();
+}
+
+void CpuCustomGBForce::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) {
+    for (int i = 0; i < numAtoms; 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) paramNames.size(); j++)
+            expressionSet.setVariable(paramIndex[j], atomParameters[i][j]);
+        for (int j = 0; j < (int) valueNames.size(); j++)
+            expressionSet.setVariable(valueIndex[j], values[j][i]);
+        if (totalEnergy != NULL)
+            *totalEnergy += (RealOpenMM) energyExpressions[index].evaluate();
+        for (int j = 0; j < (int) valueNames.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();
+    }
+}
+
+void CpuCustomGBForce::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) {
+    if (cutoff) {
+        // Loop over all pairs in the neighbor list.
+
+        for (int i = 0; i < (int) neighborList->size(); i++) {
+            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);
+        }
+    }
+    else {
+        // Perform an O(N^2) loop over all atom pairs.
+
+        for (int i = 0; i < numAtoms; i++) {
+            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);
+           }
+        }
+    }
+}
+
+void CpuCustomGBForce::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) {
+    // Compute the displacement.
+
+    RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+    if (periodic)
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxSize, deltaR);
+    else
+        ReferenceForce::getDeltaR(atomCoordinates[atom2], atomCoordinates[atom1], deltaR);
+    RealOpenMM r = deltaR[ReferenceForce::RIndex];
+    if (cutoff && r >= cutoffDistance)
+        return;
+
+    // Record variables for evaluating expressions.
+
+    for (int i = 0; i < (int) paramNames.size(); i++) {
+        expressionSet.setVariable(particleParamIndex[i*2], atomParameters[atom1][i]);
+        expressionSet.setVariable(particleParamIndex[i*2+1], atomParameters[atom2][i]);
+    }
+    expressionSet.setVariable(rindex, r);
+    for (int i = 0; i < (int) valueNames.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();
+    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();
+        dEdV[i][atom2] += (RealOpenMM) energyDerivExpressions[index][2*i+2].evaluate();
+    }
+}
+
+void CpuCustomGBForce::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) {
+    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);
+        }
+    }
+    else {
+        // Perform an O(N^2) loop over all atom pairs.
+
+        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);
+           }
+        }
+    }
+
+    // Compute chain rule terms for computed values that depend explicitly on particle coordinates.
+
+    for (int i = 0; i < numAtoms; i++) {
+        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++)
+            expressionSet.setVariable(paramIndex[j], atomParameters[i][j]);
+        for (int j = 1; j < (int) valueNames.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();
+                dVdX[j] += dVdV*dVdX[k];
+                dVdY[j] += dVdV*dVdY[k];
+                dVdZ[j] += dVdV*dVdZ[k];
+            }
+            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];
+        }
+    }
+}
+
+void CpuCustomGBForce::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) {
+    // Compute the displacement.
+
+    RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+    if (periodic)
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxSize, deltaR);
+    else
+        ReferenceForce::getDeltaR(atomCoordinates[atom2], atomCoordinates[atom1], deltaR);
+    RealOpenMM r = deltaR[ReferenceForce::RIndex];
+    if (cutoff && r >= cutoffDistance)
+        return;
+
+    // Record variables for evaluating expressions.
+
+    for (int i = 0; i < (int) paramNames.size(); i++) {
+        expressionSet.setVariable(particleParamIndex[i*2], atomParameters[atom1][i]);
+        expressionSet.setVariable(particleParamIndex[i*2+1], atomParameters[atom2][i]);
+    }
+    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.
+
+    RealOpenMM rinv = 1/r;
+    deltaR[0] *= rinv;
+    deltaR[1] *= rinv;
+    deltaR[2] *= rinv;
+    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();
+        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];
+        }
+    }
+    for (int i = 0; i < (int) paramNames.size(); i++)
+        expressionSet.setVariable(paramIndex[i], atomParameters[atom1][i]);
+    expressionSet.setVariable(valueIndex[0], values[0][atom1]);
+    for (int i = 1; i < (int) valueNames.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();
+            dVdR1[i] += dVdV*dVdR1[j];
+            dVdR2[i] += dVdV*dVdR2[j];
+        }
+        for (int k = 0; k < 3; k++) {
+            forces[atom1][k] -= dEdV[i][atom1]*dVdR1[i]*deltaR[k];
+            forces[atom2][k] -= dEdV[i][atom1]*dVdR2[i]*deltaR[k];
+        }
+    }
+}

platforms/cpu/src/CpuKernelFactory.cpp

@@ -53,6 +53,8 @@ KernelImpl* CpuKernelFactory::createKernelImpl(std::string name, const Platform&
         return new CpuCalcCustomManyParticleForceKernel(name, platform, data);
     if (name == CalcGBSAOBCForceKernel::Name())
         return new CpuCalcGBSAOBCForceKernel(name, platform, data);
+    if (name == CalcCustomGBForceKernel::Name())
+        return new CpuCalcCustomGBForceKernel(name, platform, data);
     if (name == IntegrateLangevinStepKernel::Name())
         return new CpuIntegrateLangevinStepKernel(name, platform, data);
     throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '") + name + "'").c_str());

platforms/cpu/src/CpuKernels.cpp

@@ -835,6 +835,162 @@ void CpuCalcGBSAOBCForceKernel::copyParametersToContext(ContextImpl& context, co
     obc.setParticleParameters(particleParams);
 }
 
+CpuCalcCustomGBForceKernel::~CpuCalcCustomGBForceKernel() {
+    if (particleParamArray != NULL) {
+        for (int i = 0; i < numParticles; i++)
+            delete[] particleParamArray[i];
+        delete[] particleParamArray;
+    }
+    if (neighborList != NULL)
+        delete neighborList;
+}
+
+void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGBForce& force) {
+    if (force.getNumComputedValues() > 0) {
+        string name, expression;
+        CustomGBForce::ComputationType type;
+        force.getComputedValueParameters(0, name, expression, type);
+        if (type == CustomGBForce::SingleParticle)
+            throw OpenMMException("CpuPlatform requires that the first computed value for a CustomGBForce be of type ParticlePair or ParticlePairNoExclusions.");
+        for (int i = 1; i < force.getNumComputedValues(); i++) {
+            force.getComputedValueParameters(i, name, expression, type);
+            if (type != CustomGBForce::SingleParticle)
+                throw OpenMMException("CpuPlatform requires that a CustomGBForce only have one computed value of type ParticlePair or ParticlePairNoExclusions.");
+        }
+    }
+
+    // Record the exclusions.
+
+    numParticles = force.getNumParticles();
+    exclusions.resize(numParticles);
+    for (int i = 0; i < force.getNumExclusions(); i++) {
+        int particle1, particle2;
+        force.getExclusionParticles(i, particle1, particle2);
+        exclusions[particle1].insert(particle2);
+        exclusions[particle2].insert(particle1);
+    }
+
+    // Build the arrays.
+
+    int numPerParticleParameters = force.getNumPerParticleParameters();
+    particleParamArray = new double*[numParticles];
+    for (int i = 0; i < numParticles; i++)
+        particleParamArray[i] = new double[numPerParticleParameters];
+    for (int i = 0; i < numParticles; ++i) {
+        vector<double> parameters;
+        force.getParticleParameters(i, parameters);
+        for (int j = 0; j < numPerParticleParameters; j++)
+            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
+    }
+    for (int i = 0; i < numPerParticleParameters; i++)
+        particleParameterNames.push_back(force.getPerParticleParameterName(i));
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        globalParameterNames.push_back(force.getGlobalParameterName(i));
+    nonbondedMethod = CalcCustomGBForceKernel::NonbondedMethod(force.getNonbondedMethod());
+    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
+    if (nonbondedMethod == NoCutoff)
+        neighborList = NULL;
+    else
+        neighborList = new NeighborList();
+
+    // Create custom functions for the tabulated functions.
+
+    map<string, Lepton::CustomFunction*> functions;
+    for (int i = 0; i < force.getNumFunctions(); i++)
+        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
+
+    // Parse the expressions for computed values.
+
+    valueDerivExpressions.resize(force.getNumComputedValues());
+    valueGradientExpressions.resize(force.getNumComputedValues());
+    for (int i = 0; i < force.getNumComputedValues(); i++) {
+        string name, expression;
+        CustomGBForce::ComputationType type;
+        force.getComputedValueParameters(i, name, expression, type);
+        Lepton::ParsedExpression ex = Lepton::Parser::parse(expression, functions).optimize();
+        valueExpressions.push_back(ex.createCompiledExpression());
+        valueTypes.push_back(type);
+        valueNames.push_back(name);
+        if (i == 0)
+            valueDerivExpressions[i].push_back(ex.differentiate("r").createCompiledExpression());
+        else {
+            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]).createCompiledExpression());
+        }
+    }
+
+    // Parse the expressions for energy terms.
+
+    energyDerivExpressions.resize(force.getNumEnergyTerms());
+    energyGradientExpressions.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.createCompiledExpression());
+        energyTypes.push_back(type);
+        if (type != CustomGBForce::SingleParticle)
+            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]).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());
+            }
+            else {
+                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").createCompiledExpression());
+                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").createCompiledExpression());
+            }
+        }
+    }
+
+    // Delete the custom functions.
+
+    for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
+        delete iter->second;
+}
+
+double CpuCalcCustomGBForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    vector<RealVec>& posData = extractPositions(context);
+    vector<RealVec>& forceData = extractForces(context);
+    RealOpenMM energy = 0;
+    CpuCustomGBForce ixn(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueNames, valueTypes, energyExpressions,
+        energyDerivExpressions, energyGradientExpressions, energyTypes, particleParameterNames);
+    bool periodic = (nonbondedMethod == CutoffPeriodic);
+    if (periodic)
+        ixn.setPeriodic(extractBoxSize(context));
+    if (nonbondedMethod != NoCutoff) {
+        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxSize(context), periodic, nonbondedCutoff, 0.0);
+        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
+    }
+    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);
+    return energy;
+}
+
+void CpuCalcCustomGBForceKernel::copyParametersToContext(ContextImpl& context, const CustomGBForce& force) {
+    if (numParticles != force.getNumParticles())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+
+    // Record the values.
+
+    int numParameters = force.getNumPerParticleParameters();
+    vector<double> params;
+    for (int i = 0; i < numParticles; ++i) {
+        vector<double> parameters;
+        force.getParticleParameters(i, parameters);
+        for (int j = 0; j < numParameters; j++)
+            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
+    }
+}
+
 CpuCalcCustomManyParticleForceKernel::~CpuCalcCustomManyParticleForceKernel() {
     if (particleParamArray != NULL) {
         for (int i = 0; i < numParticles; i++)

platforms/cpu/src/CpuPlatform.cpp

@@ -67,6 +67,7 @@ CpuPlatform::CpuPlatform() {
     registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomGBForceKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     platformProperties.push_back(CpuThreads());
     int threads = getNumProcessors();

platforms/cpu/tests/TestCpuCustomGBForce.cpp

+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2014 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests all the different force terms in the reference implementation of CustomGBForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "sfmt/SFMT.h"
+#include "openmm/Context.h"
+#include "CpuPlatform.h"
+#include "openmm/CustomGBForce.h"
+#include "openmm/GBSAOBCForce.h"
+#include "openmm/GBVIForce.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include <iostream>
+#include <vector>
+#include <algorithm>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMethod customMethod) {
+    const int numMolecules = 70;
+    const int numParticles = numMolecules*2;
+    const double boxSize = 10.0;
+    const double cutoff = 2.0;
+    CpuPlatform platform;
+
+    // Create two systems: one with a GBSAOBCForce, and one using a CustomGBForce to implement the same interaction.
+
+    System standardSystem;
+    System customSystem;
+    for (int i = 0; i < numParticles; i++) {
+        standardSystem.addParticle(1.0);
+        customSystem.addParticle(1.0);
+    }
+    standardSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
+    customSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
+    GBSAOBCForce* obc = new GBSAOBCForce();
+    CustomGBForce* custom = new CustomGBForce();
+    obc->setCutoffDistance(cutoff);
+    custom->setCutoffDistance(cutoff);
+    custom->addPerParticleParameter("q");
+    custom->addPerParticleParameter("radius");
+    custom->addPerParticleParameter("scale");
+    custom->addGlobalParameter("solventDielectric", obc->getSolventDielectric());
+    custom->addGlobalParameter("soluteDielectric", obc->getSoluteDielectric());
+    custom->addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r-sr2*sr2/r)+0.5*log(L/U)/r+C);"
+                                  "U=r+sr2;"
+                                  "C=2*(1/or1-1/L)*step(sr2-r-or1);"
+                                  "L=max(or1, D);"
+                                  "D=abs(r-sr2);"
+                                  "sr2 = scale2*or2;"
+                                  "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePairNoExclusions);
+    custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
+                                  "psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle);
+    custom->addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce::SingleParticle);
+    string invCutoffString = "";
+    if (obcMethod != GBSAOBCForce::NoCutoff) {
+        stringstream s;
+        s<<(1.0/cutoff);
+        invCutoffString = s.str();
+    }
+    custom->addEnergyTerm("138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2*("+invCutoffString+"-1/f);"
+                          "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePairNoExclusions);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    vector<double> params(3);
+    for (int i = 0; i < numMolecules; i++) {
+        if (i < numMolecules/2) {
+            obc->addParticle(1.0, 0.2, 0.5);
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.5;
+            custom->addParticle(params);
+            obc->addParticle(-1.0, 0.1, 0.5);
+            params[0] = -1.0;
+            params[1] = 0.1;
+            custom->addParticle(params);
+        }
+        else {
+            obc->addParticle(1.0, 0.2, 0.8);
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.8;
+            custom->addParticle(params);
+            obc->addParticle(-1.0, 0.1, 0.8);
+            params[0] = -1.0;
+            params[1] = 0.1;
+            custom->addParticle(params);
+        }
+        positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
+        velocities[2*i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
+        velocities[2*i+1] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
+    }
+    obc->setNonbondedMethod(obcMethod);
+    custom->setNonbondedMethod(customMethod);
+    standardSystem.addForce(obc);
+    customSystem.addForce(custom);
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    Context context1(standardSystem, integrator1, platform);
+    context1.setPositions(positions);
+    context1.setVelocities(velocities);
+    State state1 = context1.getState(State::Forces | State::Energy);
+    Context context2(customSystem, integrator2, platform);
+    context2.setPositions(positions);
+    context2.setVelocities(velocities);
+    State state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
+    }
+    
+    // Try changing the particle parameters and make sure it's still correct.
+    
+    for (int i = 0; i < numMolecules/2; i++) {
+        obc->setParticleParameters(2*i, 1.1, 0.3, 0.6);
+        params[0] = 1.1;
+        params[1] = 0.3;
+        params[2] = 0.6;
+        custom->setParticleParameters(2*i, params);
+        obc->setParticleParameters(2*i+1, -1.1, 0.2, 0.4);
+        params[0] = -1.1;
+        params[1] = 0.2;
+        params[2] = 0.4;
+        custom->setParticleParameters(2*i+1, params);
+    }
+    obc->updateParametersInContext(context1);
+    custom->updateParametersInContext(context2);
+    state1 = context1.getState(State::Forces | State::Energy);
+    state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
+    }
+}
+
+void testMembrane() {
+    const int numMolecules = 70;
+    const int numParticles = numMolecules*2;
+    const double boxSize = 10.0;
+    CpuPlatform platform;
+
+    // Create a system with an implicit membrane.
+
+    System system;
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+    }
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
+    CustomGBForce* custom = new CustomGBForce();
+    custom->setCutoffDistance(2.0);
+    custom->addPerParticleParameter("q");
+    custom->addPerParticleParameter("radius");
+    custom->addPerParticleParameter("scale");
+    custom->addGlobalParameter("thickness", 3);
+    custom->addGlobalParameter("solventDielectric", 78.3);
+    custom->addGlobalParameter("soluteDielectric", 1);
+    custom->addComputedValue("Imol", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r-sr2*sr2/r)+0.5*log(L/U)/r+C);"
+                             "U=r+sr2;"
+                             "C=2*(1/or1-1/L)*step(sr2-r-or1);"
+                             "L=max(or1, D);"
+                             "D=abs(r-sr2);"
+                             "sr2 = scale2*or2;"
+                             "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePairNoExclusions);
+    custom->addComputedValue("Imem", "(1/radius+2*log(2)/thickness)/(1+exp(7.2*(abs(z)+radius-0.5*thickness)))", CustomGBForce::SingleParticle);
+    custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
+                             "psi=max(Imol,Imem)*or; or=radius-0.009", CustomGBForce::SingleParticle);
+    custom->addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce::SingleParticle);
+    custom->addEnergyTerm("-138.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
+                          "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePairNoExclusions);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<double> params(3);
+    for (int i = 0; i < numMolecules; i++) {
+        if (i < numMolecules/2) {
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.5;
+            custom->addParticle(params);
+            params[0] = -1.0;
+            params[1] = 0.1;
+            custom->addParticle(params);
+        }
+        else {
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.8;
+            custom->addParticle(params);
+            params[0] = -1.0;
+            params[1] = 0.1;
+            custom->addParticle(params);
+        }
+        positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
+        velocities[2*i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
+        velocities[2*i+1] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
+    }
+    system.addForce(custom);
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+
+    // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
+
+    double norm = 0.0;
+    for (int i = 0; i < (int) forces.size(); ++i)
+        norm += forces[i].dot(forces[i]);
+    norm = std::sqrt(norm);
+    const double stepSize = 1e-2;
+    double step = 0.5*stepSize/norm;
+    vector<Vec3> positions2(numParticles), positions3(numParticles);
+    for (int i = 0; i < (int) positions.size(); ++i) {
+        Vec3 p = positions[i];
+        Vec3 f = forces[i];
+        positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+        positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
+    }
+    context.setPositions(positions2);
+    State state2 = context.getState(State::Energy);
+    context.setPositions(positions3);
+    State state3 = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-3);
+}
+
+void testTabulatedFunction() {
+    CpuPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomGBForce* force = new CustomGBForce();
+    force->addComputedValue("a", "0", CustomGBForce::ParticlePair);
+    force->addEnergyTerm("fn(r)+1", CustomGBForce::ParticlePair);
+    force->addParticle(vector<double>());
+    force->addParticle(vector<double>());
+    vector<double> table;
+    for (int i = 0; i < 21; i++)
+        table.push_back(std::sin(0.25*i));
+    force->addTabulatedFunction("fn", new Continuous1DFunction(table, 1.0, 6.0));
+    system.addForce(force);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    for (int i = 1; i < 30; i++) {
+        double x = (7.0/30.0)*i;
+        positions[1] = Vec3(x, 0, 0);
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces = state.getForces();
+        double force = (x < 1.0 || x > 6.0 ? 0.0 : -std::cos(x-1.0));
+        double energy = (x < 1.0 || x > 6.0 ? 0.0 : std::sin(x-1.0))+1.0;
+        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], 0.1);
+        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], 0.1);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.02);
+    }
+}
+
+void testMultipleChainRules() {
+    CpuPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomGBForce* force = new CustomGBForce();
+    force->addComputedValue("a", "2*r", CustomGBForce::ParticlePair);
+    force->addComputedValue("b", "a+1", CustomGBForce::SingleParticle);
+    force->addComputedValue("c", "2*b+a", CustomGBForce::SingleParticle);
+    force->addEnergyTerm("0.1*a+1*b+10*c", CustomGBForce::SingleParticle); // 0.1*(2*r) + 2*r+1 + 10*(3*a+2) = 0.2*r + 2*r+1 + 40*r+20+20*r = 62.2*r+21
+    force->addParticle(vector<double>());
+    force->addParticle(vector<double>());
+    system.addForce(force);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    for (int i = 1; i < 5; i++) {
+        positions[1] = Vec3(i, 0, 0);
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces = state.getForces();
+        ASSERT_EQUAL_VEC(Vec3(124.4, 0, 0), forces[0], 1e-4);
+        ASSERT_EQUAL_VEC(Vec3(-124.4, 0, 0), forces[1], 1e-4);
+        ASSERT_EQUAL_TOL(2*(62.2*i+21), state.getPotentialEnergy(), 0.02);
+    }
+}
+
+void testPositionDependence() {
+    CpuPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomGBForce* force = new CustomGBForce();
+    force->addComputedValue("a", "r", CustomGBForce::ParticlePair);
+    force->addComputedValue("b", "a+x*y", CustomGBForce::SingleParticle);
+    force->addEnergyTerm("b*z", CustomGBForce::SingleParticle);
+    force->addEnergyTerm("b1+b2", CustomGBForce::ParticlePair); // = 2*r+x1*y1+x2*y2
+    force->addParticle(vector<double>());
+    force->addParticle(vector<double>());
+    system.addForce(force);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    vector<Vec3> forces(2);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i < 5; i++) {
+        positions[0] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
+        positions[1] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces = state.getForces();
+        Vec3 delta = positions[0]-positions[1];
+        double r = sqrt(delta.dot(delta));
+        double energy = 2*r+positions[0][0]*positions[0][1]+positions[1][0]*positions[1][1];
+        for (int j = 0; j < 2; j++)
+            energy += positions[j][2]*(r+positions[j][0]*positions[j][1]);
+        Vec3 force1(-(1+positions[0][2])*delta[0]/r-(1+positions[0][2])*positions[0][1]-(1+positions[1][2])*delta[0]/r,
+                    -(1+positions[0][2])*delta[1]/r-(1+positions[0][2])*positions[0][0]-(1+positions[1][2])*delta[1]/r,
+                    -(1+positions[0][2])*delta[2]/r-(r+positions[0][0]*positions[0][1])-(1+positions[1][2])*delta[2]/r);
+        Vec3 force2((1+positions[0][2])*delta[0]/r+(1+positions[1][2])*delta[0]/r-(1+positions[1][2])*positions[1][1],
+                    (1+positions[0][2])*delta[1]/r+(1+positions[1][2])*delta[1]/r-(1+positions[1][2])*positions[1][0],
+                    (1+positions[0][2])*delta[2]/r+(1+positions[1][2])*delta[2]/r-(r+positions[1][0]*positions[1][1]));
+        ASSERT_EQUAL_VEC(force1, forces[0], 1e-4);
+        ASSERT_EQUAL_VEC(force2, forces[1], 1e-4);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.02);
+
+        // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
+
+        double norm = 0.0;
+        for (int i = 0; i < (int) forces.size(); ++i)
+            norm += forces[i].dot(forces[i]);
+        norm = std::sqrt(norm);
+        const double stepSize = 1e-3;
+        double step = 0.5*stepSize/norm;
+        vector<Vec3> positions2(2), positions3(2);
+        for (int i = 0; i < (int) positions.size(); ++i) {
+            Vec3 p = positions[i];
+            Vec3 f = forces[i];
+            positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+            positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
+        }
+        context.setPositions(positions2);
+        State state2 = context.getState(State::Energy);
+        context.setPositions(positions3);
+        State state3 = context.getState(State::Energy);
+        ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-3);
+    }
+}
+
+void testExclusions() {
+    CpuPlatform platform;
+    for (int i = 3; i < 4; i++) {
+        System system;
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        VerletIntegrator integrator(0.01);
+        CustomGBForce* force = new CustomGBForce();
+        force->addComputedValue("a", "r", i < 2 ? CustomGBForce::ParticlePair : CustomGBForce::ParticlePairNoExclusions);
+        force->addEnergyTerm("a", CustomGBForce::SingleParticle);
+        force->addEnergyTerm("(1+a1+a2)*r", i%2 == 0 ? CustomGBForce::ParticlePair : CustomGBForce::ParticlePairNoExclusions);
+        force->addParticle(vector<double>());
+        force->addParticle(vector<double>());
+        force->addExclusion(0, 1);
+        system.addForce(force);
+        Context context(system, integrator, platform);
+        vector<Vec3> positions(2);
+        positions[0] = Vec3(0, 0, 0);
+        positions[1] = Vec3(1, 0, 0);
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces = state.getForces();
+        double f, energy;
+        switch (i)
+        {
+            case 0: // e = 0
+                f = 0;
+                energy = 0;
+                break;
+            case 1: // e = r
+                f = 1;
+                energy = 1;
+                break;
+            case 2: // e = 2r
+                f = 2;
+                energy = 2;
+                break;
+            case 3: // e = 3r + 2r^2
+                f = 7;
+                energy = 5;
+                break;
+            default:
+                ASSERT(false);
+        }
+        ASSERT_EQUAL_VEC(Vec3(f, 0, 0), forces[0], 1e-4);
+        ASSERT_EQUAL_VEC(Vec3(-f, 0, 0), forces[1], 1e-4);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 1e-4);
+
+        // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
+
+        double norm = 0.0;
+        for (int i = 0; i < (int) forces.size(); ++i)
+            norm += forces[i].dot(forces[i]);
+        norm = std::sqrt(norm);
+        const double stepSize = 1e-3;
+        double step = stepSize/norm;
+        for (int i = 0; i < (int) positions.size(); ++i) {
+            Vec3 p = positions[i];
+            Vec3 f = forces[i];
+            positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+        }
+        context.setPositions(positions);
+        State state2 = context.getState(State::Energy);
+        ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/stepSize, 1e-3*abs(state.getPotentialEnergy()));
+    }
+}
+
+int main() {
+    try {
+        testOBC(GBSAOBCForce::NoCutoff, CustomGBForce::NoCutoff);
+        testOBC(GBSAOBCForce::CutoffNonPeriodic, CustomGBForce::CutoffNonPeriodic);
+        testOBC(GBSAOBCForce::CutoffPeriodic, CustomGBForce::CutoffPeriodic);
+        testMembrane();
+        testTabulatedFunction();
+        testMultipleChainRules();
+        testPositionDependence();
+        testExclusions();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}