Merge pull request #644 from peastman/cpugb

Created CPU implementation 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 "CompiledExpressionSet.h"
+#include "CpuNeighborList.h"
+#include "lepton/CompiledExpression.h"
+#include "openmm/CustomGBForce.h"
+#include "openmm/internal/ThreadPool.h"
+#include "openmm/internal/vectorize.h"
+#include <map>
+#include <set>
+#include <vector>
+
+namespace OpenMM {
+
+class CpuCustomGBForce {
+private:
+    class ComputeForceTask;
+    class ThreadData;
+
+    bool cutoff;
+    bool periodic;
+    const CpuNeighborList* neighborList;
+    float periodicBoxSize[3];
+    float cutoffDistance, cutoffDistance2;
+    const std::vector<std::set<int> > exclusions;
+    std::vector<std::string> valueNames;
+    std::vector<CustomGBForce::ComputationType> valueTypes;
+    std::vector<std::string> paramNames;
+    std::vector<CustomGBForce::ComputationType> energyTypes;
+    ThreadPool& threads;
+    std::vector<ThreadData*> threadData;
+    std::vector<double> threadEnergy;
+    // Workspace vectors
+    std::vector<std::vector<float> > values, dEdV;
+    // The following variables are used to make information accessible to the individual threads.
+    int numberOfAtoms;
+    float* posq;
+    RealOpenMM** atomParameters;
+    const std::map<std::string, double>* globalParameters;
+    std::vector<AlignedArray<float> >* threadForce;
+    bool includeForce, includeEnergy;
+    void* atomicCounter;
+    
+    /**
+     * This routine contains the code executed by each thread.
+     */
+    void threadComputeForce(ThreadPool& threads, int threadIndex);
+
+    /**
+     * Calculate a computed value that is based on particle pairs
+     * 
+     * @param index            the index of the value to compute
+     * @param data             workspace for the current thread
+     * @param numAtoms         number of atoms
+     * @param posq             atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param useExclusions    specifies whether to use exclusions
+     */
+
+    void calculateParticlePairValue(int index, ThreadData& data, int numAtoms, float* posq, RealOpenMM** atomParameters,
+                                    bool useExclusions, const fvec4& boxSize, const fvec4& invBoxSize);
+
+    /**
+     * 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 data             workspace for the current thread
+     * @param posq             atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     */
+
+    void calculateOnePairValue(int index, int atom1, int atom2, ThreadData& data, float* posq, RealOpenMM** atomParameters,
+                               std::vector<float>& valueArray, const fvec4& boxSize, const fvec4& invBoxSize);
+
+    /**
+     * Calculate an energy term of type SingleParticle
+     * 
+     * @param index            the index of the value to compute
+     * @param data             workspace for the current thread
+     * @param numAtoms         number of atoms
+     * @param posq             atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param forces           forces on atoms are added to this
+     * @param totalEnergy      the energy contribution is added to this
+     */
+
+    void calculateSingleParticleEnergyTerm(int index, ThreadData& data, int numAtoms, float* posq, RealOpenMM** atomParameters, float* forces, double& totalEnergy);
+
+    /**
+     * Calculate an energy term that is based on particle pairs
+     * 
+     * @param index            the index of the term to compute
+     * @param data             workspace for the current thread
+     * @param numAtoms         number of atoms
+     * @param posq             atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @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
+     */
+
+    void calculateParticlePairEnergyTerm(int index, ThreadData& data, int numAtoms, float* posq, RealOpenMM** atomParameters,
+                                    bool useExclusions, float* forces, double& totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize);
+
+    /**
+     * 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 data             workspace for the current thread
+     * @param posq             atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param forces           forces on atoms are added to this
+     * @param totalEnergy      the energy contribution is added to this
+     */
+
+    void calculateOnePairEnergyTerm(int index, int atom1, int atom2, ThreadData& data, float* posq, RealOpenMM** atomParameters,
+                               float* forces, double& totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize);
+
+    /**
+     * Apply the chain rule to compute forces on atoms
+     * 
+     * @param data             workspace for the current thread
+     * @param numAtoms         number of atoms
+     * @param posq             atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param forces           forces on atoms are added to this
+     */
+
+    void calculateChainRuleForces(ThreadData& data, int numAtoms, float* posq, RealOpenMM** atomParameters,
+                                    float* forces, const fvec4& boxSize, const fvec4& invBoxSize);
+
+    /**
+     * 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 data             workspace for the current thread
+     * @param posq             atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param forces           forces on atoms are added to this
+     * @param isExcluded       specifies whether this is an excluded pair
+     */
+
+    void calculateOnePairChainRule(int atom1, int atom2, ThreadData& data, float* posq, RealOpenMM** atomParameters,
+                               float* forces, bool isExcluded, const fvec4& boxSize, const fvec4& invBoxSize);
+
+    /**
+     * Compute the displacement and squared distance between two points, optionally using
+     * periodic boundary conditions.
+     */
+    void getDeltaR(const fvec4& posI, const fvec4& posJ, fvec4& deltaR, float& r2, bool periodic, const fvec4& boxSize, const fvec4& invBoxSize) const;
+
+public:
+
+    /**
+     * Construct a new CpuCustomGBForce.
+     */
+
+     CpuCustomGBForce(int numAtoms, const std::vector<std::set<int> >& exclusions,
+                        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<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<CustomGBForce::ComputationType>& energyTypes,
+                        const std::vector<std::string>& parameterNames, ThreadPool& threads);
+
+     ~CpuCustomGBForce();
+
+    /**
+     * Set the force to use a cutoff.
+     * 
+     * @param distance            the cutoff distance
+     * @param neighbors           the neighbor list to use
+     */
+
+    void setUseCutoff(float distance, const CpuNeighborList& 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(RealVec& boxSize);
+
+    /**
+     * Calculate custom GB ixn
+     * 
+     * @param numberOfAtoms    number of atoms
+     * @param posq             atom coordinates
+     * @param atomParameters   atomParameters[atomIndex][paramterIndex]
+     * @param globalParameters the values of global parameters
+     * @param forces           force array (forces added)
+     * @param totalEnergy      total energy
+     */
+
+    void calculateIxn(int numberOfAtoms, float* posq, RealOpenMM** atomParameters,
+                     std::map<std::string, double>& globalParameters, std::vector<AlignedArray<float> >& threadForce, bool includeForce, bool includeEnergy, double& totalEnergy);
+};
+
+class CpuCustomGBForce::ThreadData {
+public:
+    ThreadData(int numAtoms, int numThreads, int threadIndex,
+               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<Lepton::CompiledExpression>& energyExpressions,
+               const std::vector<std::vector<Lepton::CompiledExpression> >& energyDerivExpressions,
+               const std::vector<std::vector<Lepton::CompiledExpression> >& energyGradientExpressions,
+               const std::vector<std::string>& parameterNames);
+    CompiledExpressionSet expressionSet;
+    std::vector<Lepton::CompiledExpression> valueExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > valueDerivExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > valueGradientExpressions;
+    std::vector<int> valueIndex;
+    std::vector<Lepton::CompiledExpression> energyExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > energyDerivExpressions;
+    std::vector<std::vector<Lepton::CompiledExpression> > energyGradientExpressions;
+    std::vector<int> paramIndex;
+    std::vector<int> particleParamIndex;
+    std::vector<int> particleValueIndex;
+    int xindex, yindex, zindex, rindex;
+    int firstAtom, lastAtom;
+    // Workspace vectors
+    std::vector<float> value0, dVdR1, dVdR2, dVdX, dVdY, dVdZ;
+    std::vector<std::vector<float> > dEdV;
+};
+
+} // namespace OpenMM
+
+#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,53 @@ 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;
+    CpuCustomGBForce* ixn;
+    std::vector<std::set<int> > exclusions;
+    std::vector<std::string> particleParameterNames, globalParameterNames, valueNames;
+    std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
+    std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
+    NonbondedMethod nonbondedMethod;
+    CpuNeighborList* 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"
+#include "gmx_atomic.h"
+
+using namespace OpenMM;
+using namespace std;
+
+class CpuCustomGBForce::ComputeForceTask : public ThreadPool::Task {
+public:
+    ComputeForceTask(CpuCustomGBForce& owner) : owner(owner) {
+    }
+    void execute(ThreadPool& threads, int threadIndex) {
+        owner.threadComputeForce(threads, threadIndex);
+    }
+    CpuCustomGBForce& owner;
+};
+
+CpuCustomGBForce::ThreadData::ThreadData(int numAtoms, int numThreads, int threadIndex,
+                      const vector<Lepton::CompiledExpression>& valueExpressions,
+                      const vector<vector<Lepton::CompiledExpression> >& valueDerivExpressions,
+                      const vector<vector<Lepton::CompiledExpression> >& valueGradientExpressions,
+                      const vector<string>& valueNames,
+                      const vector<Lepton::CompiledExpression>& energyExpressions,
+                      const vector<vector<Lepton::CompiledExpression> >& energyDerivExpressions,
+                      const vector<vector<Lepton::CompiledExpression> >& energyGradientExpressions,
+                      const vector<string>& parameterNames) :
+            valueExpressions(valueExpressions), valueDerivExpressions(valueDerivExpressions), valueGradientExpressions(valueGradientExpressions),
+            energyExpressions(energyExpressions), energyDerivExpressions(energyDerivExpressions), energyGradientExpressions(energyGradientExpressions) {
+    firstAtom = (threadIndex*(long long) numAtoms)/numThreads;
+    lastAtom = ((threadIndex+1)*(long long) numAtoms)/numThreads;
+    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) parameterNames.size(); i++) {
+        paramIndex.push_back(expressionSet.getVariableIndex(parameterNames[i]));
+        for (int j = 1; j < 3; j++) {
+            stringstream name;
+            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;
+            particleValueIndex.push_back(expressionSet.getVariableIndex(name.str()));
+        }
+    }
+    value0.resize(numAtoms);
+    dEdV.resize(valueNames.size());
+    for (int i = 0; i < (int) dEdV.size(); i++)
+        dEdV[i].resize(numAtoms);
+    dVdX.resize(valueDerivExpressions.size());
+    dVdY.resize(valueDerivExpressions.size());
+    dVdZ.resize(valueDerivExpressions.size());
+    dVdR1.resize(valueDerivExpressions.size());
+    dVdR2.resize(valueDerivExpressions.size());
+}
+
+CpuCustomGBForce::CpuCustomGBForce(int numAtoms, const std::vector<std::set<int> >& exclusions,
+                     const vector<Lepton::CompiledExpression>& valueExpressions,
+                     const vector<vector<Lepton::CompiledExpression> >& valueDerivExpressions,
+                     const vector<vector<Lepton::CompiledExpression> >& valueGradientExpressions,
+                     const vector<string>& valueNames,
+                     const vector<CustomGBForce::ComputationType>& valueTypes,
+                     const vector<Lepton::CompiledExpression>& energyExpressions,
+                     const vector<vector<Lepton::CompiledExpression> >& energyDerivExpressions,
+                     const vector<vector<Lepton::CompiledExpression> >& energyGradientExpressions,
+                     const vector<CustomGBForce::ComputationType>& energyTypes,
+                     const vector<string>& parameterNames, ThreadPool& threads) :
+            exclusions(exclusions), cutoff(false), periodic(false), valueNames(valueNames), valueTypes(valueTypes),
+            energyTypes(energyTypes), paramNames(parameterNames), threads(threads) {
+    for (int i = 0; i < threads.getNumThreads(); i++)
+        threadData.push_back(new ThreadData(numAtoms, threads.getNumThreads(), i, valueExpressions, valueDerivExpressions, valueGradientExpressions, valueNames,
+                      energyExpressions, energyDerivExpressions, energyGradientExpressions, parameterNames));
+    values.resize(valueNames.size());
+    dEdV.resize(valueNames.size());
+    for (int i = 0; i < (int) values.size(); i++) {
+        values[i].resize(numAtoms);
+        dEdV[i].resize(numAtoms);
+    }
+}
+
+CpuCustomGBForce::~CpuCustomGBForce() {
+    for (int i = 0; i < (int) threadData.size(); i++)
+        delete threadData[i];
+}
+
+void CpuCustomGBForce::setUseCutoff(float distance, const CpuNeighborList& neighbors) {
+    cutoff = true;
+    cutoffDistance = distance;
+    cutoffDistance2 = distance*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, float* posq, RealOpenMM** atomParameters,
+                                           map<string, double>& globalParameters, vector<AlignedArray<float> >& threadForce,
+                                           bool includeForce, bool includeEnergy, double& totalEnergy) {
+    // Record the parameters for the threads.
+    
+    this->numberOfAtoms = numberOfAtoms;
+    this->posq = posq;
+    this->atomParameters = atomParameters;
+    this->globalParameters = &globalParameters;
+    this->threadForce = &threadForce;
+    this->includeForce = includeForce;
+    this->includeEnergy = includeEnergy;
+    threadEnergy.resize(threads.getNumThreads());
+    gmx_atomic_t counter;
+    this->atomicCounter = &counter;
+
+    // Calculate the first computed value.
+
+    ComputeForceTask task(*this);
+    gmx_atomic_set(&counter, 0);
+    threads.execute(task);
+    threads.waitForThreads();
+    
+    // Calculate the remaining computed values.
+    
+    threads.resumeThreads();
+    threads.waitForThreads();
+
+    // Calculate the energy terms.
+
+    for (int i = 0; i < (int) threadData[0]->energyExpressions.size(); i++) {
+        gmx_atomic_set(&counter, 0);
+        threads.execute(task);
+        threads.waitForThreads();
+    }
+
+    // Sum the energy derivatives.
+
+    threads.resumeThreads();
+    threads.waitForThreads();
+    
+    // Apply the chain rule to evaluate forces.
+
+    gmx_atomic_set(&counter, 0);
+    threads.resumeThreads();
+    threads.waitForThreads();
+
+    // Combine the energies from all the threads.
+    
+    if (includeEnergy) {
+        int numThreads = threads.getNumThreads();
+        for (int i = 0; i < numThreads; i++)
+            totalEnergy += threadEnergy[i];
+    }
+}
+
+void CpuCustomGBForce::threadComputeForce(ThreadPool& threads, int threadIndex) {
+    // Compute this thread's subset of interactions.
+
+    int numThreads = threads.getNumThreads();
+    threadEnergy[threadIndex] = 0;
+    double& energy = threadEnergy[threadIndex];
+    float* forces = &(*threadForce)[threadIndex][0];
+    ThreadData& data = *threadData[threadIndex];
+    fvec4 boxSize(periodicBoxSize[0], periodicBoxSize[1], periodicBoxSize[2], 0);
+    fvec4 invBoxSize((1/periodicBoxSize[0]), (1/periodicBoxSize[1]), (1/periodicBoxSize[2]), 0);
+    for (map<string, double>::const_iterator iter = globalParameters->begin(); iter != globalParameters->end(); ++iter)
+        data.expressionSet.setVariable(data.expressionSet.getVariableIndex(iter->first), iter->second);
+
+    // Calculate the first computed value.
+
+    for (int i = 0; i < (int) data.value0.size(); i++)
+        data.value0[i] = 0.0f;
+    if (valueTypes[0] == CustomGBForce::ParticlePair)
+        calculateParticlePairValue(0, data, numberOfAtoms, posq, atomParameters, true, boxSize, invBoxSize);
+    else
+        calculateParticlePairValue(0, data, numberOfAtoms, posq, atomParameters, false, boxSize, invBoxSize);
+    threads.syncThreads();
+
+    // Sum the first computed value and calculate the remaining ones.
+
+    int numValues = valueTypes.size();
+    for (int atom = data.firstAtom; atom < data.lastAtom; atom++) {
+        float sum = 0.0f;
+        for (int j = 0; j < (int) threadData.size(); j++)
+            sum += threadData[j]->value0[atom];
+        values[0][atom] = sum;
+        data.expressionSet.setVariable(data.xindex, posq[4*atom]);
+        data.expressionSet.setVariable(data.yindex, posq[4*atom+1]);
+        data.expressionSet.setVariable(data.zindex, posq[4*atom+2]);
+        for (int j = 0; j < (int) paramNames.size(); j++)
+            data.expressionSet.setVariable(data.paramIndex[j], atomParameters[atom][j]);
+        for (int i = 1; i < numValues; i++) {
+            data.expressionSet.setVariable(data.valueIndex[i-1], values[i-1][atom]);
+            values[i][atom] = (float) data.valueExpressions[i].evaluate();
+        }
+    }
+    threads.syncThreads();
+
+    // Now calculate the energy and its derivatives.
+
+    for (int i = 0; i < (int) data.dEdV.size(); i++)
+        for (int j = 0; j < (int) data.dEdV[i].size(); j++)
+            data.dEdV[i][j] = 0.0;
+    for (int termIndex = 0; termIndex < (int) data.energyExpressions.size(); termIndex++) {
+        if (energyTypes[termIndex] == CustomGBForce::SingleParticle)
+            calculateSingleParticleEnergyTerm(termIndex, data, numberOfAtoms, posq, atomParameters, forces, energy);
+        else if (energyTypes[termIndex] == CustomGBForce::ParticlePair)
+            calculateParticlePairEnergyTerm(termIndex, data, numberOfAtoms, posq, atomParameters, true, forces, energy, boxSize, invBoxSize);
+        else
+            calculateParticlePairEnergyTerm(termIndex, data, numberOfAtoms, posq, atomParameters, false, forces, energy, boxSize, invBoxSize);
+        threads.syncThreads();
+    }
+
+    // Sum the energy derivatives.
+
+    for (int atom = data.firstAtom; atom < data.lastAtom; atom++) {
+        for (int i = 0; i < (int) dEdV.size(); i++) {
+            float sum = 0.0f;
+            for (int j = 0; j < (int) threadData.size(); j++)
+                sum += threadData[j]->dEdV[i][atom];
+            dEdV[i][atom] = sum;
+        }
+    }
+    threads.syncThreads();
+
+    // Apply the chain rule to evaluate forces.
+
+    calculateChainRuleForces(data, numberOfAtoms, posq, atomParameters, forces, boxSize, invBoxSize);
+}
+
+void CpuCustomGBForce::calculateParticlePairValue(int index, ThreadData& data, int numAtoms, float* posq, RealOpenMM** atomParameters,
+        bool useExclusions, const fvec4& boxSize, const fvec4& invBoxSize) {
+    for (int i = 0; i < numAtoms; i++)
+        values[index][i] = 0.0f;
+    vector<float>& valueArray = (index == 0 ? data.value0 : values[index]);
+    if (cutoff) {
+        // Loop over all pairs in the neighbor list.
+
+        while (true) {
+            int blockIndex = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
+            if (blockIndex >= neighborList->getNumBlocks())
+                break;
+            const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
+            const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
+            const vector<char>& blockExclusions = neighborList->getBlockExclusions(blockIndex);
+            for (int i = 0; i < (int) neighbors.size(); i++) {
+                int first = neighbors[i];
+                for (int k = 0; k < 4; k++) {
+                    if ((blockExclusions[i] & (1<<k)) == 0) {
+                        int second = blockAtom[k];
+                        if (useExclusions && exclusions[first].find(second) != exclusions[first].end())
+                            continue;
+                        calculateOnePairValue(index, first, second, data, posq, atomParameters, valueArray, boxSize, invBoxSize);
+                        calculateOnePairValue(index, second, first, data, posq, atomParameters, valueArray, boxSize, invBoxSize);
+                    }
+                }
+            }
+        }
+    }
+    else {
+        // Perform an O(N^2) loop over all atom pairs.
+
+        while (true) {
+            int i = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
+            if (i >= numAtoms)
+                break;
+            for (int j = i+1; j < numAtoms; j++) {
+                if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
+                    continue;
+                calculateOnePairValue(index, i, j, data, posq, atomParameters, valueArray, boxSize, invBoxSize);
+                calculateOnePairValue(index, j, i, data, posq, atomParameters, valueArray, boxSize, invBoxSize);
+           }
+        }
+    }
+}
+
+void CpuCustomGBForce::calculateOnePairValue(int index, int atom1, int atom2, ThreadData& data, float* posq, RealOpenMM** atomParameters,
+        vector<float>& valueArray, const fvec4& boxSize, const fvec4& invBoxSize) {
+    fvec4 deltaR;
+    fvec4 pos1(posq+4*atom1);
+    fvec4 pos2(posq+4*atom2);
+    float r2;
+    getDeltaR(pos2, pos1, deltaR, r2, periodic, boxSize, invBoxSize);
+    if (cutoff && r2 >= cutoffDistance2)
+        return;
+    float r = sqrtf(r2);
+    for (int i = 0; i < (int) paramNames.size(); i++) {
+        data.expressionSet.setVariable(data.particleParamIndex[i*2], atomParameters[atom1][i]);
+        data.expressionSet.setVariable(data.particleParamIndex[i*2+1], atomParameters[atom2][i]);
+    }
+    data.expressionSet.setVariable(data.rindex, r);
+    for (int i = 0; i < index; i++) {
+        data.expressionSet.setVariable(data.particleValueIndex[i*2], values[i][atom1]);
+        data.expressionSet.setVariable(data.particleValueIndex[i*2+1], values[i][atom2]);
+    }
+    valueArray[atom1] += (float) data.valueExpressions[index].evaluate();
+}
+
+void CpuCustomGBForce::calculateSingleParticleEnergyTerm(int index, ThreadData& data, int numAtoms, float* posq,
+        RealOpenMM** atomParameters, float* forces, double& totalEnergy) {
+    for (int i = data.firstAtom; i < data.lastAtom; i++) {
+        data.expressionSet.setVariable(data.xindex, posq[4*i]);
+        data.expressionSet.setVariable(data.yindex, posq[4*i+1]);
+        data.expressionSet.setVariable(data.zindex, posq[4*i+2]);
+        for (int j = 0; j < (int) paramNames.size(); j++)
+            data.expressionSet.setVariable(data.paramIndex[j], atomParameters[i][j]);
+        for (int j = 0; j < (int) valueNames.size(); j++)
+            data.expressionSet.setVariable(data.valueIndex[j], values[j][i]);
+        if (includeEnergy)
+            totalEnergy += (float) data.energyExpressions[index].evaluate();
+        for (int j = 0; j < (int) valueNames.size(); j++)
+            data.dEdV[j][i] += (float) data.energyDerivExpressions[index][j].evaluate();
+        forces[4*i+0] -= (float) data.energyGradientExpressions[index][0].evaluate();
+        forces[4*i+1] -= (float) data.energyGradientExpressions[index][1].evaluate();
+        forces[4*i+2] -= (float) data.energyGradientExpressions[index][2].evaluate();
+    }
+}
+
+void CpuCustomGBForce::calculateParticlePairEnergyTerm(int index, ThreadData& data, int numAtoms, float* posq, RealOpenMM** atomParameters,
+        bool useExclusions, float* forces, double& totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize) {
+    if (cutoff) {
+        // Loop over all pairs in the neighbor list.
+
+        while (true) {
+            int blockIndex = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
+            if (blockIndex >= neighborList->getNumBlocks())
+                break;
+            const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
+            const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
+            const vector<char>& blockExclusions = neighborList->getBlockExclusions(blockIndex);
+            for (int i = 0; i < (int) neighbors.size(); i++) {
+                int first = neighbors[i];
+                for (int k = 0; k < 4; k++) {
+                    if ((blockExclusions[i] & (1<<k)) == 0) {
+                        int second = blockAtom[k];
+                        if (useExclusions && exclusions[first].find(second) != exclusions[first].end())
+                            continue;
+                        calculateOnePairEnergyTerm(index, first, second, data, posq, atomParameters, forces, totalEnergy, boxSize, invBoxSize);
+                    }
+                }
+            }
+        }
+    }
+    else {
+        // Perform an O(N^2) loop over all atom pairs.
+
+        while (true) {
+            int i = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
+            if (i >= numAtoms)
+                break;
+            for (int j = i+1; j < numAtoms; j++) {
+                if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
+                    continue;
+                calculateOnePairEnergyTerm(index, i, j, data, posq, atomParameters, forces, totalEnergy, boxSize, invBoxSize);
+           }
+        }
+    }
+}
+
+void CpuCustomGBForce::calculateOnePairEnergyTerm(int index, int atom1, int atom2, ThreadData& data, float* posq, RealOpenMM** atomParameters,
+        float* forces, double& totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize) {
+    // Compute the displacement.
+
+    fvec4 deltaR;
+    fvec4 pos1(posq+4*atom1);
+    fvec4 pos2(posq+4*atom2);
+    float r2;
+    getDeltaR(pos2, pos1, deltaR, r2, periodic, boxSize, invBoxSize);
+    if (cutoff && r2 >= cutoffDistance2)
+        return;
+    float r = sqrtf(r2);
+
+    // Record variables for evaluating expressions.
+
+    for (int i = 0; i < (int) paramNames.size(); i++) {
+        data.expressionSet.setVariable(data.particleParamIndex[i*2], atomParameters[atom1][i]);
+        data.expressionSet.setVariable(data.particleParamIndex[i*2+1], atomParameters[atom2][i]);
+    }
+    data.expressionSet.setVariable(data.rindex, r);
+    for (int i = 0; i < (int) valueNames.size(); i++) {
+        data.expressionSet.setVariable(data.particleValueIndex[i*2], values[i][atom1]);
+        data.expressionSet.setVariable(data.particleValueIndex[i*2+1], values[i][atom2]);
+    }
+
+    // Evaluate the energy and its derivatives.
+
+    if (includeEnergy)
+        totalEnergy += (float) data.energyExpressions[index].evaluate();
+    float dEdR = (float) data.energyDerivExpressions[index][0].evaluate();
+    dEdR *= 1/r;
+    fvec4 result = deltaR*dEdR;
+    (fvec4(forces+4*atom1)-result).store(forces+4*atom1);
+    (fvec4(forces+4*atom2)+result).store(forces+4*atom2);
+    for (int i = 0; i < (int) valueNames.size(); i++) {
+        data.dEdV[i][atom1] += (float) data.energyDerivExpressions[index][2*i+1].evaluate();
+        data.dEdV[i][atom2] += (float) data.energyDerivExpressions[index][2*i+2].evaluate();
+    }
+}
+
+void CpuCustomGBForce::calculateChainRuleForces(ThreadData& data, int numAtoms, float* posq, RealOpenMM** atomParameters,
+        float* forces, const fvec4& boxSize, const fvec4& invBoxSize) {
+    if (cutoff) {
+        // Loop over all pairs in the neighbor list.
+
+        while (true) {
+            int blockIndex = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
+            if (blockIndex >= neighborList->getNumBlocks())
+                break;
+            const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
+            const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
+            const vector<char>& blockExclusions = neighborList->getBlockExclusions(blockIndex);
+            for (int i = 0; i < (int) neighbors.size(); i++) {
+                int first = neighbors[i];
+                for (int k = 0; k < 4; k++) {
+                    if ((blockExclusions[i] & (1<<k)) == 0) {
+                        int second = blockAtom[k];
+                        bool isExcluded = (exclusions[first].find(second) != exclusions[first].end());
+                        calculateOnePairChainRule(first, second, data, posq, atomParameters, forces, isExcluded, boxSize, invBoxSize);
+                        calculateOnePairChainRule(second, first, data, posq, atomParameters, forces, isExcluded, boxSize, invBoxSize);
+                    }
+                }
+            }
+        }
+    }
+    else {
+        // Perform an O(N^2) loop over all atom pairs.
+
+        while (true) {
+            int i = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
+            if (i >= numAtoms)
+                break;
+            for (int j = i+1; j < numAtoms; j++) {
+                bool isExcluded = (exclusions[i].find(j) != exclusions[i].end());
+                calculateOnePairChainRule(i, j, data, posq, atomParameters, forces, isExcluded, boxSize, invBoxSize);
+                calculateOnePairChainRule(j, i, data, posq, atomParameters, forces, isExcluded, boxSize, invBoxSize);
+           }
+        }
+    }
+
+    // Compute chain rule terms for computed values that depend explicitly on particle coordinates.
+
+    for (int i = data.firstAtom; i < data.lastAtom; i++) {
+        data.expressionSet.setVariable(data.xindex, posq[4*i]);
+        data.expressionSet.setVariable(data.yindex, posq[4*i+1]);
+        data.expressionSet.setVariable(data.zindex, posq[4*i+2]);
+        for (int j = 0; j < (int) paramNames.size(); j++)
+            data.expressionSet.setVariable(data.paramIndex[j], atomParameters[i][j]);
+        for (int j = 1; j < (int) valueNames.size(); j++) {
+            data.expressionSet.setVariable(data.valueIndex[j-1], values[j-1][i]);
+            data.dVdX[j] = 0.0;
+            data.dVdY[j] = 0.0;
+            data.dVdZ[j] = 0.0;
+            for (int k = 1; k < j; k++) {
+                float dVdV = (float) data.valueDerivExpressions[j][k].evaluate();
+                data.dVdX[j] += dVdV*data.dVdX[k];
+                data.dVdY[j] += dVdV*data.dVdY[k];
+                data.dVdZ[j] += dVdV*data.dVdZ[k];
+            }
+            data.dVdX[j] += (float) data.valueGradientExpressions[j][0].evaluate();
+            data.dVdY[j] += (float) data.valueGradientExpressions[j][1].evaluate();
+            data.dVdZ[j] += (float) data.valueGradientExpressions[j][2].evaluate();
+            forces[4*i+0] -= dEdV[j][i]*data.dVdX[j];
+            forces[4*i+1] -= dEdV[j][i]*data.dVdY[j];
+            forces[4*i+2] -= dEdV[j][i]*data.dVdZ[j];
+        }
+    }
+}
+
+void CpuCustomGBForce::calculateOnePairChainRule(int atom1, int atom2, ThreadData& data, float* posq, RealOpenMM** atomParameters,
+        float* forces, bool isExcluded, const fvec4& boxSize, const fvec4& invBoxSize) {
+    // Compute the displacement.
+
+    fvec4 deltaR;
+    fvec4 pos1(posq+4*atom1);
+    fvec4 pos2(posq+4*atom2);
+    float r2;
+    getDeltaR(pos2, pos1, deltaR, r2, periodic, boxSize, invBoxSize);
+    if (cutoff && r2 >= cutoffDistance2)
+        return;
+    float r = sqrtf(r2);
+
+    // Record variables for evaluating expressions.
+
+    for (int i = 0; i < (int) paramNames.size(); i++) {
+        data.expressionSet.setVariable(data.particleParamIndex[i*2], atomParameters[atom1][i]);
+        data.expressionSet.setVariable(data.particleParamIndex[i*2+1], atomParameters[atom2][i]);
+        data.expressionSet.setVariable(data.paramIndex[i], atomParameters[atom1][i]);
+    }
+    data.expressionSet.setVariable(data.valueIndex[0], values[0][atom1]);
+    data.expressionSet.setVariable(data.xindex, posq[4*atom1]);
+    data.expressionSet.setVariable(data.yindex, posq[4*atom1+1]);
+    data.expressionSet.setVariable(data.zindex, posq[4*atom1+2]);
+    data.expressionSet.setVariable(data.rindex, r);
+    data.expressionSet.setVariable(data.particleValueIndex[0], values[0][atom1]);
+    data.expressionSet.setVariable(data.particleValueIndex[1], values[0][atom2]);
+
+    // Evaluate the derivative of each parameter with respect to position and apply forces.
+
+    float rinv = 1/r;
+    deltaR *= rinv;
+    fvec4 f1(0.0f), f2(0.0f);
+    if (!isExcluded || valueTypes[0] != CustomGBForce::ParticlePair) {
+        data.dVdR1[0] = (float) data.valueDerivExpressions[0][0].evaluate();
+        data.dVdR2[0] = -data.dVdR1[0];
+        f1 -= deltaR*(dEdV[0][atom1]*data.dVdR1[0]);
+        f2 -= deltaR*(dEdV[0][atom1]*data.dVdR2[0]);
+    }
+    for (int i = 1; i < (int) valueNames.size(); i++) {
+        data.expressionSet.setVariable(data.valueIndex[i], values[i][atom1]);
+        data.dVdR1[i] = 0.0;
+        data.dVdR2[i] = 0.0;
+        for (int j = 0; j < i; j++) {
+            float dVdV = (float) data.valueDerivExpressions[i][j].evaluate();
+            data.dVdR1[i] += dVdV*data.dVdR1[j];
+            data.dVdR2[i] += dVdV*data.dVdR2[j];
+        }
+        f1 -= deltaR*(dEdV[i][atom1]*data.dVdR1[i]);
+        f2 -= deltaR*(dEdV[i][atom1]*data.dVdR2[i]);
+    }
+    (fvec4(forces+4*atom1)+f1).store(forces+4*atom1);
+    (fvec4(forces+4*atom2)+f2).store(forces+4*atom2);
+}
+
+void CpuCustomGBForce::getDeltaR(const fvec4& posI, const fvec4& posJ, fvec4& deltaR, float& r2, bool periodic, const fvec4& boxSize, const fvec4& invBoxSize) const {
+    deltaR = posJ-posI;
+    if (periodic) {
+        fvec4 base = round(deltaR*invBoxSize)*boxSize;
+        deltaR = deltaR-base;
+    }
+    r2 = dot3(deltaR, deltaR);
+}

platforms/cpu/src/CpuCustomNonbondedForce.cpp

@@ -182,9 +182,7 @@ void CpuCustomNonbondedForce::threadComputeForce(ThreadPool& threads, int thread
             int blockIndex = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
             if (blockIndex >= neighborList->getNumBlocks())
                 break;
-            int blockAtom[4];
-            for (int i = 0; i < 4; i++)
-                blockAtom[i] = neighborList->getSortedAtoms()[4*blockIndex+i];
+            const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
             const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
             const vector<char>& exclusions = neighborList->getBlockExclusions(blockIndex);
             for (int i = 0; i < (int) neighbors.size(); i++) {

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,168 @@ 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;
+    if (ixn != NULL)
+        delete ixn;
+}
+
+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 CpuNeighborList(4);
+
+    // 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.
+
+    vector<vector<Lepton::CompiledExpression> > valueDerivExpressions(force.getNumComputedValues());
+    vector<vector<Lepton::CompiledExpression> > valueGradientExpressions(force.getNumComputedValues());
+    vector<Lepton::CompiledExpression> valueExpressions;
+    vector<Lepton::CompiledExpression> energyExpressions;
+    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.
+
+    vector<vector<Lepton::CompiledExpression> > energyDerivExpressions(force.getNumEnergyTerms());
+    vector<vector<Lepton::CompiledExpression> > energyGradientExpressions(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;
+    ixn = new CpuCustomGBForce(numParticles, exclusions, valueExpressions, valueDerivExpressions, valueGradientExpressions, valueNames, valueTypes, energyExpressions,
+        energyDerivExpressions, energyGradientExpressions, energyTypes, particleParameterNames, data.threads);
+    data.isPeriodic = (force.getNonbondedMethod() == CustomGBForce::CutoffPeriodic);
+}
+
+double CpuCalcCustomGBForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    vector<RealVec>& forceData = extractForces(context);
+    RealOpenMM energy = 0;
+    RealVec& box = extractBoxSize(context);
+    float floatBoxSize[3] = {(float) box[0], (float) box[1], (float) box[2]};
+    if (data.isPeriodic)
+        ixn->setPeriodic(extractBoxSize(context));
+    if (nonbondedMethod != NoCutoff) {
+        vector<set<int> > noExclusions(numParticles);
+        neighborList->computeNeighborList(numParticles, data.posq, exclusions, floatBoxSize, data.isPeriodic, nonbondedCutoff, data.threads);
+        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, &data.posq[0], particleParamArray, globalParameters, data.threadForce, includeForces, includeEnergy, energy);
+    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/CpuNonbondedForceVec4.cpp

@@ -48,13 +48,11 @@ CpuNonbondedForceVec4::CpuNonbondedForceVec4() {
 void CpuNonbondedForceVec4::calculateBlockIxn(int blockIndex, float* forces, double* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize) {
     // Load the positions and parameters of the atoms in the block.
     
-    int blockAtom[4];
+    const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
     fvec4 blockAtomPosq[4];
     fvec4 blockAtomForceX(0.0f), blockAtomForceY(0.0f), blockAtomForceZ(0.0f);
-    for (int i = 0; i < 4; i++) {
-        blockAtom[i] = neighborList->getSortedAtoms()[4*blockIndex+i];
+    for (int i = 0; i < 4; i++)
         blockAtomPosq[i] = fvec4(posq+4*blockAtom[i]);
-    }
     fvec4 blockAtomX = fvec4(blockAtomPosq[0][0], blockAtomPosq[1][0], blockAtomPosq[2][0], blockAtomPosq[3][0]);
     fvec4 blockAtomY = fvec4(blockAtomPosq[0][1], blockAtomPosq[1][1], blockAtomPosq[2][1], blockAtomPosq[3][1]);
     fvec4 blockAtomZ = fvec4(blockAtomPosq[0][2], blockAtomPosq[1][2], blockAtomPosq[2][2], blockAtomPosq[3][2]);
@@ -159,13 +157,11 @@ void CpuNonbondedForceVec4::calculateBlockIxn(int blockIndex, float* forces, dou
 void CpuNonbondedForceVec4::calculateBlockEwaldIxn(int blockIndex, float* forces, double* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize) {
     // Load the positions and parameters of the atoms in the block.
     
-    int blockAtom[4];
+    const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
     fvec4 blockAtomPosq[4];
     fvec4 blockAtomForceX(0.0f), blockAtomForceY(0.0f), blockAtomForceZ(0.0f);
-    for (int i = 0; i < 4; i++) {
-        blockAtom[i] = neighborList->getSortedAtoms()[4*blockIndex+i];
+    for (int i = 0; i < 4; i++)
         blockAtomPosq[i] = fvec4(posq+4*blockAtom[i]);
-    }
     fvec4 blockAtomX = fvec4(blockAtomPosq[0][0], blockAtomPosq[1][0], blockAtomPosq[2][0], blockAtomPosq[3][0]);
     fvec4 blockAtomY = fvec4(blockAtomPosq[0][1], blockAtomPosq[1][1], blockAtomPosq[2][1], blockAtomPosq[3][1]);
     fvec4 blockAtomZ = fvec4(blockAtomPosq[0][2], blockAtomPosq[1][2], blockAtomPosq[2][2], blockAtomPosq[3][2]);

platforms/cpu/src/CpuNonbondedForceVec8.cpp

@@ -74,14 +74,12 @@ CpuNonbondedForceVec8::CpuNonbondedForceVec8() {
 void CpuNonbondedForceVec8::calculateBlockIxn(int blockIndex, float* forces, double* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize) {
     // Load the positions and parameters of the atoms in the block.
     
-    int blockAtom[8];
+    const int* blockAtom = &neighborList->getSortedAtoms()[8*blockIndex];
     fvec4 blockAtomPosq[8];
     fvec8 blockAtomForceX(0.0f), blockAtomForceY(0.0f), blockAtomForceZ(0.0f);
     fvec8 blockAtomX, blockAtomY, blockAtomZ, blockAtomCharge;
-    for (int i = 0; i < 8; i++) {
-        blockAtom[i] = neighborList->getSortedAtoms()[8*blockIndex+i];
+    for (int i = 0; i < 8; i++)
         blockAtomPosq[i] = fvec4(posq+4*blockAtom[i]);
-    }
     transpose(blockAtomPosq[0], blockAtomPosq[1], blockAtomPosq[2], blockAtomPosq[3], blockAtomPosq[4], blockAtomPosq[5], blockAtomPosq[6], blockAtomPosq[7], blockAtomX, blockAtomY, blockAtomZ, blockAtomCharge);
     blockAtomCharge *= ONE_4PI_EPS0;
     fvec8 blockAtomSigma(atomParameters[blockAtom[0]].first, atomParameters[blockAtom[1]].first, atomParameters[blockAtom[2]].first, atomParameters[blockAtom[3]].first, atomParameters[blockAtom[4]].first, atomParameters[blockAtom[5]].first, atomParameters[blockAtom[6]].first, atomParameters[blockAtom[7]].first);
@@ -184,14 +182,12 @@ void CpuNonbondedForceVec8::calculateBlockIxn(int blockIndex, float* forces, dou
 void CpuNonbondedForceVec8::calculateBlockEwaldIxn(int blockIndex, float* forces, double* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize) {
     // Load the positions and parameters of the atoms in the block.
     
-    int blockAtom[8];
+    const int* blockAtom = &neighborList->getSortedAtoms()[8*blockIndex];
     fvec4 blockAtomPosq[8];
     fvec8 blockAtomForceX(0.0f), blockAtomForceY(0.0f), blockAtomForceZ(0.0f);
     fvec8 blockAtomX, blockAtomY, blockAtomZ, blockAtomCharge;
-    for (int i = 0; i < 8; i++) {
-        blockAtom[i] = neighborList->getSortedAtoms()[8*blockIndex+i];
+    for (int i = 0; i < 8; i++)
         blockAtomPosq[i] = fvec4(posq+4*blockAtom[i]);
-    }
     transpose(blockAtomPosq[0], blockAtomPosq[1], blockAtomPosq[2], blockAtomPosq[3], blockAtomPosq[4], blockAtomPosq[5], blockAtomPosq[6], blockAtomPosq[7], blockAtomX, blockAtomY, blockAtomZ, blockAtomCharge);
     blockAtomCharge *= ONE_4PI_EPS0;
     fvec8 blockAtomSigma(atomParameters[blockAtom[0]].first, atomParameters[blockAtom[1]].first, atomParameters[blockAtom[2]].first, atomParameters[blockAtom[3]].first, atomParameters[blockAtom[4]].first, atomParameters[blockAtom[5]].first, atomParameters[blockAtom[6]].first, atomParameters[blockAtom[7]].first);

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;
+}