OpenCL version of CustomHbondForce is now complete (except for exclusions)

openmmapi/include/openmm/internal/CustomHbondForceImpl.h

@@ -35,6 +35,7 @@
 #include "ForceImpl.h"
 #include "openmm/CustomHbondForce.h"
 #include "openmm/Kernel.h"
+#include "lepton/CustomFunction.h"
 #include "lepton/ExpressionTreeNode.h"
 #include "lepton/ParsedExpression.h"
 #include <utility>
@@ -68,15 +69,16 @@ public:
      * as follows: 0=a1, 1=a2, 2=a3, 3=d1, 4=d2, 5=d3.
      *
      * @param force     the CustomHbondForce to process
-     * @param distances on exist, this will contain an entry for each distance used in the expression.  The key is the name
+     * @param functions definitions of custom function that may appear in the expression
+     * @param distances on exit, this will contain an entry for each distance used in the expression.  The key is the name
      *                  of the corresponding variable, and the value is the list of particle indices.
-     * @param angles    on exist, this will contain an entry for each angle used in the expression.  The key is the name
+     * @param angles    on exit, this will contain an entry for each angle used in the expression.  The key is the name
      *                  of the corresponding variable, and the value is the list of particle indices.
-     * @param dihedrals on exist, this will contain an entry for each dihedral used in the expression.  The key is the name
+     * @param dihedrals on exit, this will contain an entry for each dihedral used in the expression.  The key is the name
      *                  of the corresponding variable, and the value is the list of particle indices.
      * @return a Parsed expression for the energy
      */
-    static Lepton::ParsedExpression prepareExpression(const CustomHbondForce& force, std::map<std::string, std::vector<int> >& distances,
+    static Lepton::ParsedExpression prepareExpression(const CustomHbondForce& force, const std::map<std::string, Lepton::CustomFunction*>& functions, std::map<std::string, std::vector<int> >& distances,
             std::map<std::string, std::vector<int> >& angles, std::map<std::string, std::vector<int> >& dihedrals);
 private:
     class FunctionPlaceholder;

openmmapi/src/CustomHbondForceImpl.cpp

@@ -199,24 +199,16 @@ map<string, double> CustomHbondForceImpl::getDefaultParameters() {
     return parameters;
 }
 
-ParsedExpression CustomHbondForceImpl::prepareExpression(const CustomHbondForce& force, map<string, vector<int> >& distances,
+ParsedExpression CustomHbondForceImpl::prepareExpression(const CustomHbondForce& force, const map<string, CustomFunction*>& customFunctions, map<string, vector<int> >& distances,
         map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
     CustomHbondForceImpl::FunctionPlaceholder custom(1);
     CustomHbondForceImpl::FunctionPlaceholder distance(2);
     CustomHbondForceImpl::FunctionPlaceholder angle(3);
     CustomHbondForceImpl::FunctionPlaceholder dihedral(4);
-    map<string, CustomFunction*> functions;
+    map<string, CustomFunction*> functions = customFunctions;
     functions["distance"] = &distance;
     functions["angle"] = &angle;
     functions["dihedral"] = &dihedral;
-    for (int i = 0; i < force.getNumFunctions(); i++) {
-        string name;
-        vector<double> values;
-        double min, max;
-        bool interpolating;
-        force.getFunctionParameters(i, name, values, min, max, interpolating);
-        functions[name] = &custom;
-    }
     ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions);
     map<string, int> atoms;
     atoms["a1"] = 0;

platforms/opencl/src/OpenCLKernels.cpp

@@ -29,6 +29,7 @@
 #include "openmm/LangevinIntegrator.h"
 #include "openmm/Context.h"
 #include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/CustomHbondForceImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
 #include "OpenCLExpressionUtilities.h"
 #include "OpenCLIntegrationUtilities.h"
@@ -39,6 +40,7 @@
 #include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
 #include "lepton/Operation.h"
 #include <cmath>
+#include <set>
 
 using namespace OpenMM;
 using namespace std;
@@ -56,13 +58,6 @@ static string intToString(int value) {
     return s.str();
 }
 
-static bool isZeroExpression(const Lepton::ParsedExpression& expression) {
-    const Lepton::Operation& op = expression.getRootNode().getOperation();
-    if (op.getId() != Lepton::Operation::CONSTANT)
-        return false;
-    return (dynamic_cast<const Lepton::Operation::Constant&>(op).getValue() == 0.0);
-}
-
 void OpenCLCalcForcesAndEnergyKernel::initialize(const System& system) {
 }
 
@@ -2459,33 +2454,41 @@ public:
         vector<double> parameters;
         if (index < force.getNumDonors()) {
             force.getDonorParameters(index, p1, p2, p3, parameters);
-            particles.resize(3);
-            particles[0] = p1;
-            particles[1] = p2;
-            particles[2] = p3;
+            particles.clear();
+            particles.push_back(p1);
+            if (p2 > -1)
+                particles.push_back(p2);
+            if (p3 > -1)
+                particles.push_back(p3);
             return;
         }
         index -= force.getNumDonors();
         if (index < force.getNumAcceptors()) {
             force.getAcceptorParameters(index, p1, p2, p3, parameters);
-            particles.resize(3);
-            particles[0] = p1;
-            particles[1] = p2;
-            particles[2] = p3;
+            particles.clear();
+            particles.push_back(p1);
+            if (p2 > -1)
+                particles.push_back(p2);
+            if (p3 > -1)
+                particles.push_back(p3);
             return;
         }
         index -= force.getNumAcceptors();
         int donor, acceptor;
         force.getExclusionParticles(index, donor, acceptor);
-        particles.resize(6);
+        particles.clear();
         force.getDonorParameters(donor, p1, p2, p3, parameters);
-        particles[0] = p1;
-        particles[1] = p2;
-        particles[2] = p3;
+        particles.push_back(p1);
+        if (p2 > -1)
+            particles.push_back(p2);
+        if (p3 > -1)
+            particles.push_back(p3);
         force.getAcceptorParameters(acceptor, p1, p2, p3, parameters);
-        particles[3] = p1;
-        particles[4] = p2;
-        particles[5] = p3;
+        particles.push_back(p1);
+        if (p2 > -1)
+            particles.push_back(p2);
+        if (p3 > -1)
+            particles.push_back(p3);
     }
     bool areGroupsIdentical(int group1, int group2) {
         int p1, p2, p3;
@@ -2533,12 +2536,20 @@ OpenCLCalcCustomHbondForceKernel::~OpenCLCalcCustomHbondForceKernel() {
         delete tabulatedFunctions[i];
 }
 
-void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const CustomHbondForce& force) {
-    int forceIndex;
-    for (forceIndex = 0; forceIndex < system.getNumForces() && &system.getForce(forceIndex) != &force; ++forceIndex)
-        ;
-    string prefix = "custom"+intToString(forceIndex)+"_";
+static void addDonorAndAcceptorCode(stringstream& computeDonor, stringstream& computeAcceptor, const string& value) {
+    computeDonor << value;
+    computeAcceptor << value;
+}
+
+static void applyDonorAndAcceptorForces(stringstream& applyToDonor, stringstream& applyToAcceptor, int atom, const string& value) {
+    string forceNames[] = {"f1", "f2", "f3"};
+    if (atom < 3)
+        applyToAcceptor << forceNames[atom]<<".xyz += "<<value<<";\n";
+    else
+        applyToDonor << forceNames[atom-3]<<".xyz += "<<value<<";\n";
+}
 
+void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const CustomHbondForce& force) {
     // Record the lists of donors and acceptors, and the parameters for each one.
 
     numDonors = force.getNumDonors();
@@ -2573,19 +2584,26 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
     acceptors->upload(acceptorVector);
     acceptorParams->setParameterValues(acceptorParamVector);
 
-    // Select a output buffer indices for each donor and acceptor.
+    // Select an output buffer index for each donor and acceptor.
 
     donorBufferIndices = new OpenCLArray<mm_int4>(cl, numDonors, "customHbondDonorBuffers");
     acceptorBufferIndices = new OpenCLArray<mm_int4>(cl, numAcceptors, "customHbondAcceptorBuffers");
     vector<mm_int4> donorBufferVector(numDonors);
     vector<mm_int4> acceptorBufferVector(numAcceptors);
-    vector<int> particleBuffers(numParticles, 0);
+    vector<int> donorBufferCounter(numParticles, 0);
     for (int i = 0; i < numDonors; i++)
-        donorBufferVector[i] = mm_int4(particleBuffers[donorVector[i].x]++, particleBuffers[donorVector[i].y]++, particleBuffers[donorVector[i].z]++, 0);
+        donorBufferVector[i] = mm_int4(donorBufferCounter[donorVector[i].x]++, donorBufferCounter[donorVector[i].y]++, donorBufferCounter[donorVector[i].z]++, 0);
+    vector<int> acceptorBufferCounter(numParticles, 0);
     for (int i = 0; i < numAcceptors; i++)
-        acceptorBufferVector[i] = mm_int4(particleBuffers[acceptorVector[i].x]++, particleBuffers[acceptorVector[i].y]++, particleBuffers[acceptorVector[i].z]++, 0);
+        acceptorBufferVector[i] = mm_int4(acceptorBufferCounter[acceptorVector[i].x]++, acceptorBufferCounter[acceptorVector[i].y]++, acceptorBufferCounter[acceptorVector[i].z]++, 0);
     donorBufferIndices->upload(donorBufferVector);
     acceptorBufferIndices->upload(acceptorBufferVector);
+    int maxBuffers = 1;
+    for (int i = 0; i < (int) donorBufferCounter.size(); i++)
+        maxBuffers = max(maxBuffers, donorBufferCounter[i]);
+    for (int i = 0; i < (int) acceptorBufferCounter.size(); i++)
+        maxBuffers = max(maxBuffers, acceptorBufferCounter[i]);
+    cl.addForce(new OpenCLCustomHbondForceInfo(maxBuffers, force));
 
     // Record exclusions.
 
@@ -2609,7 +2627,7 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
         double min, max;
         bool interpolating;
         force.getFunctionParameters(i, name, values, min, max, interpolating);
-        string arrayName = prefix+"table"+intToString(i);
+        string arrayName = "table"+intToString(i);
         functionDefinitions.push_back(make_pair(name, arrayName));
         functions[name] = &fp;
         tabulatedFunctionParamsVec[i] = mm_float4((float) min, (float) max, (float) ((values.size()-1)/(max-min)), 0.0f);
@@ -2621,10 +2639,10 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
     if (force.getNumFunctions() > 0) {
         tabulatedFunctionParams = new OpenCLArray<mm_float4>(cl, tabulatedFunctionParamsVec.size(), "tabulatedFunctionParameters", false, CL_MEM_READ_ONLY);
         tabulatedFunctionParams->upload(tabulatedFunctionParamsVec);
-        tableArgs << ", __constant float4* " << prefix << "functionParams";
+        tableArgs << ", __constant float4* functionParams";
     }
 
-    // Record information for the expressions.
+    // Record information about parameters.
 
     globalParamNames.resize(force.getNumGlobalParameters());
     globalParamValues.resize(force.getNumGlobalParameters());
@@ -2634,27 +2652,7 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
     }
     if (globals != NULL)
         globals->upload(globalParamValues);
-    bool useCutoff = (force.getNonbondedMethod() != CustomHbondForce::NoCutoff);
-    bool usePeriodic = (force.getNonbondedMethod() != CustomHbondForce::NoCutoff && force.getNonbondedMethod() != CustomHbondForce::CutoffNonPeriodic);
-    Lepton::ParsedExpression energyExpression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
-    Lepton::ParsedExpression rForceExpression = energyExpression.differentiate("r").optimize();
-    Lepton::ParsedExpression thetaForceExpression = energyExpression.differentiate("theta").optimize();
-    Lepton::ParsedExpression psiForceExpression = energyExpression.differentiate("psi").optimize();
-    Lepton::ParsedExpression chiForceExpression = energyExpression.differentiate("chi").optimize();
-    map<string, Lepton::ParsedExpression> forceExpressions;
-    forceExpressions["energy += "] = energyExpression;
-    forceExpressions["float dEdR = "] = rForceExpression;
-    forceExpressions["float dEdTheta = "] = thetaForceExpression;
-    forceExpressions["float dEdPsi = "] = psiForceExpression;
-    forceExpressions["float dEdChi = "] = chiForceExpression;
-
-    // Create the kernels.
-
     map<string, string> variables;
-    variables["r"] = "r";
-    variables["theta"] = "theta";
-    variables["psi"] = "psi";
-    variables["chi"] = "chi";
     for (int i = 0; i < force.getNumPerDonorParameters(); i++) {
         const string& name = force.getPerDonorParameterName(i);
         variables[name] = "donorParams"+donorParams->getParameterSuffix(i);
@@ -2667,38 +2665,172 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
         const string& name = force.getGlobalParameterName(i);
         variables[name] = "globals["+intToString(i)+"]";
     }
-    stringstream compute, extraArgs;
+
+    // Now to generate the kernel.  First, it needs to calculate all distances, angles,
+    // and dihedrals the expression depends on.
+
+    map<string, vector<int> > distances;
+    map<string, vector<int> > angles;
+    map<string, vector<int> > dihedrals;
+    Lepton::ParsedExpression energyExpression = CustomHbondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
+    map<string, Lepton::ParsedExpression> forceExpressions;
+    set<string> computedDeltas;
+    computedDeltas.insert("D1A1");
+    string atomNames[] = {"A1", "A2", "A3", "D1", "D2", "D3"};
+    string atomNamesLower[] = {"a1", "a2", "a3", "d1", "d2", "d3"};
+    stringstream computeDonor, computeAcceptor, extraArgs;
+    int index = 0;
+    for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName = atomNames[atoms[0]]+atomNames[atoms[1]];
+        if (computedDeltas.count(deltaName) == 0) {
+            addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 delta"+deltaName+" = delta("+atomNamesLower[atoms[0]]+", "+atomNamesLower[atoms[1]]+");\n");
+            computedDeltas.insert(deltaName);
+        }
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float r_"+deltaName+" = sqrt(delta"+deltaName+".w);\n");
+        variables[iter->first] = "r_"+deltaName;
+        forceExpressions["float dEdDistance"+intToString(index)+" = "] = energyExpression.differentiate(iter->first).optimize();
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName1 = atomNames[atoms[1]]+atomNames[atoms[0]];
+        string deltaName2 = atomNames[atoms[1]]+atomNames[atoms[2]];
+        string angleName = "angle_"+atomNames[atoms[0]]+atomNames[atoms[1]]+atomNames[atoms[2]];
+        if (computedDeltas.count(deltaName1) == 0) {
+            addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 delta"+deltaName1+" = delta("+atomNamesLower[atoms[1]]+", "+atomNamesLower[atoms[0]]+");\n");
+            computedDeltas.insert(deltaName1);
+        }
+        if (computedDeltas.count(deltaName2) == 0) {
+            addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 delta"+deltaName2+" = delta("+atomNamesLower[atoms[1]]+", "+atomNamesLower[atoms[2]]+");\n");
+            computedDeltas.insert(deltaName2);
+        }
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float "+angleName+" = computeAngle(delta"+deltaName1+", delta"+deltaName2+");\n");
+        variables[iter->first] = angleName;
+        forceExpressions["float dEdAngle"+intToString(index)+" = "] = energyExpression.differentiate(iter->first).optimize();
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName1 = atomNames[atoms[0]]+atomNames[atoms[1]];
+        string deltaName2 = atomNames[atoms[2]]+atomNames[atoms[1]];
+        string deltaName3 = atomNames[atoms[2]]+atomNames[atoms[3]];
+        string crossName1 = "cross_"+deltaName1+"_"+deltaName2;
+        string crossName2 = "cross_"+deltaName2+"_"+deltaName3;
+        string dihedralName = "dihedral_"+atomNames[atoms[0]]+atomNames[atoms[1]]+atomNames[atoms[2]]+atomNames[atoms[3]];
+        if (computedDeltas.count(deltaName1) == 0) {
+            addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 delta"+deltaName1+" = delta("+atomNamesLower[atoms[0]]+", "+atomNamesLower[atoms[1]]+");\n");
+            computedDeltas.insert(deltaName1);
+        }
+        if (computedDeltas.count(deltaName2) == 0) {
+            addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 delta"+deltaName2+" = delta("+atomNamesLower[atoms[2]]+", "+atomNamesLower[atoms[1]]+");\n");
+            computedDeltas.insert(deltaName2);
+        }
+        if (computedDeltas.count(deltaName3) == 0) {
+            addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 delta"+deltaName3+" = delta("+atomNamesLower[atoms[2]]+", "+atomNamesLower[atoms[3]]+");\n");
+            computedDeltas.insert(deltaName3);
+        }
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 "+crossName1+" = computeCross(delta"+deltaName1+", delta"+deltaName2+");\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 "+crossName2+" = computeCross(delta"+deltaName2+", delta"+deltaName3+");\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float "+dihedralName+" = computeAngle("+crossName1+", "+crossName2+");\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, dihedralName+" *= (delta"+deltaName1+".x*"+crossName2+".x + delta"+deltaName1+".y*"+crossName2+".y + delta"+deltaName1+".z*"+crossName2+".z < 0 ? -1 : 1);\n");
+        variables[iter->first] = dihedralName;
+        forceExpressions["float dEdDihedral"+intToString(index)+" = "] = energyExpression.differentiate(iter->first).optimize();
+    }
+
+    // Next it needs to load parameters from global memory.
+
     if (force.getNumGlobalParameters() > 0)
         extraArgs << ", __constant float* globals";
     for (int i = 0; i < (int) donorParams->getBuffers().size(); i++) {
         const OpenCLNonbondedUtilities::ParameterInfo& buffer = donorParams->getBuffers()[i];
         extraArgs << ", __global "+buffer.getType()+"* donor"+buffer.getName();
-        compute<<buffer.getType()<<" donorParams"<<(i+1)<<" = donor"<<buffer.getName()<<"[index];\n";
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, buffer.getType()+" donorParams"+intToString(i+1)+" = donor"+buffer.getName()+"[index];\n");
     }
     for (int i = 0; i < (int) acceptorParams->getBuffers().size(); i++) {
         const OpenCLNonbondedUtilities::ParameterInfo& buffer = acceptorParams->getBuffers()[i];
         extraArgs << ", __global "+buffer.getType()+"* acceptor"+buffer.getName();
-        compute<<buffer.getType()<<" acceptorParams"<<(i+1)<<" = acceptor"<<buffer.getName()<<"[index];\n";
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, buffer.getType()+" acceptorParams"+intToString(i+1)+" = acceptor"+buffer.getName()+"[index];\n");
     }
-    compute << OpenCLExpressionUtilities::createExpressions(forceExpressions, variables, functionDefinitions, "temp", "functionParams");
+
+    // Now evaluate the expressions.
+
+    computeAcceptor << OpenCLExpressionUtilities::createExpressions(forceExpressions, variables, functionDefinitions, "temp", "functionParams");
+    forceExpressions["energy += "] = energyExpression;
+    computeDonor << OpenCLExpressionUtilities::createExpressions(forceExpressions, variables, functionDefinitions, "temp", "functionParams");
+
+    // Finally, apply forces to atoms.
+
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName = atomNames[atoms[0]]+atomNames[atoms[1]];
+        string value = "(dEdDistance"+intToString(index)+"/r_"+deltaName+")*delta"+deltaName+".xyz";
+        applyDonorAndAcceptorForces(computeDonor, computeAcceptor, atoms[0], "-"+value);
+        applyDonorAndAcceptorForces(computeDonor, computeAcceptor, atoms[1], value);
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName1 = atomNames[atoms[1]]+atomNames[atoms[0]];
+        string deltaName2 = atomNames[atoms[1]]+atomNames[atoms[2]];
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "{\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 crossProd = cross(delta"+deltaName2+", delta"+deltaName1+");\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float lengthCross = max(length(crossProd), 1e-6f);\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 deltaCross0 = -cross(delta"+deltaName1+", crossProd)*dEdAngle"+intToString(index)+"/(delta"+deltaName1+".w*lengthCross);\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 deltaCross2 = cross(delta"+deltaName2+", crossProd)*dEdAngle"+intToString(index)+"/(delta"+deltaName2+".w*lengthCross);\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 deltaCross1 = -(deltaCross0+deltaCross2);\n");
+        applyDonorAndAcceptorForces(computeDonor, computeAcceptor, atoms[0], "deltaCross0.xyz");
+        applyDonorAndAcceptorForces(computeDonor, computeAcceptor, atoms[1], "deltaCross1.xyz");
+        applyDonorAndAcceptorForces(computeDonor, computeAcceptor, atoms[2], "deltaCross2.xyz");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "}\n");
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName1 = atomNames[atoms[0]]+atomNames[atoms[1]];
+        string deltaName2 = atomNames[atoms[2]]+atomNames[atoms[1]];
+        string deltaName3 = atomNames[atoms[2]]+atomNames[atoms[3]];
+        string crossName1 = "cross_"+deltaName1+"_"+deltaName2;
+        string crossName2 = "cross_"+deltaName2+"_"+deltaName3;
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "{\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float r = sqrt(delta"+deltaName2+".w);\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 ff;\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "ff.x = (-dEdDihedral"+intToString(index)+"*r)/"+crossName1+".w;\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "ff.y = (delta"+deltaName1+".x*delta"+deltaName2+".x + delta"+deltaName1+".y*delta"+deltaName2+".y + delta"+deltaName1+".z*delta"+deltaName2+".z)/delta"+deltaName2+".w;\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "ff.z = (delta"+deltaName3+".x*delta"+deltaName2+".x + delta"+deltaName3+".y*delta"+deltaName2+".y + delta"+deltaName3+".z*delta"+deltaName2+".z)/delta"+deltaName2+".w;\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "ff.w = (dEdDihedral"+intToString(index)+"*r)/"+crossName2+".w;\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 internalF0 = ff.x*"+crossName1+";\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 internalF3 = ff.w*"+crossName2+";\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "float4 s = ff.y*internalF0 - ff.z*internalF3;\n");
+        applyDonorAndAcceptorForces(computeDonor, computeAcceptor, atoms[0], "internalF0.xyz");
+        applyDonorAndAcceptorForces(computeDonor, computeAcceptor, atoms[1], "s.xyz-internalF0.xyz");
+        applyDonorAndAcceptorForces(computeDonor, computeAcceptor, atoms[2], "-s.xyz-internalF3.xyz");
+        applyDonorAndAcceptorForces(computeDonor, computeAcceptor, atoms[3], "internalF3.xyz");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, "}\n");
+    }
+
+    // Generate the kernels.
+
     map<string, string> replacements;
-    replacements["COMPUTE_FORCE"] = compute.str();
+    replacements["COMPUTE_DONOR_FORCE"] = computeDonor.str();
+    replacements["COMPUTE_ACCEPTOR_FORCE"] = computeAcceptor.str();
     replacements["PARAMETER_ARGUMENTS"] = extraArgs.str()+tableArgs.str();
     map<string, string> defines;
     defines["PADDED_NUM_ATOMS"] = intToString(cl.getPaddedNumAtoms());
     defines["NUM_DONORS"] = intToString(force.getNumDonors());
     defines["NUM_ACCEPTORS"] = intToString(force.getNumAcceptors());
     defines["M_PI"] = doubleToString(M_PI);
-    if (!isZeroExpression(rForceExpression))
-        defines["INCLUDE_R"] = "1";
-    if (!isZeroExpression(thetaForceExpression))
-        defines["INCLUDE_THETA"] = "1";
-    if (!isZeroExpression(psiForceExpression))
-        defines["INCLUDE_PSI"] = "1";
-    if (!isZeroExpression(chiForceExpression))
-        defines["INCLUDE_CHI"] = "1";
+    if (force.getNonbondedMethod() != CustomHbondForce::NoCutoff) {
+        defines["USE_CUTOFF"] = "1";
+        defines["CUTOFF_SQUARED"] = doubleToString(force.getCutoffDistance()*force.getCutoffDistance());
+    }
+    if (force.getNonbondedMethod() != CustomHbondForce::NoCutoff && force.getNonbondedMethod() != CustomHbondForce::CutoffNonPeriodic)
+        defines["USE_PERIODIC"] = "1";
+
     cl::Program program = cl.createProgram(cl.replaceStrings(OpenCLKernelSources::customHbondForce, replacements), defines);
-    kernel = cl::Kernel(program, "computeHbonds");
+    donorKernel = cl::Kernel(program, "computeDonorForces");
+    acceptorKernel = cl::Kernel(program, "computeAcceptorForces");
 }
 
 void OpenCLCalcCustomHbondForceKernel::executeForces(ContextImpl& context) {
@@ -2716,27 +2848,54 @@ void OpenCLCalcCustomHbondForceKernel::executeForces(ContextImpl& context) {
     if (!hasInitializedKernel) {
         hasInitializedKernel = true;
         int index = 0;
-        kernel.setArg<cl::Buffer>(index++, cl.getForceBuffers().getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(index++, cl.getEnergyBuffer().getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(index++, cl.getPosq().getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(index++, donors->getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(index++, acceptors->getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(index++, donorBufferIndices->getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(index++, acceptorBufferIndices->getDeviceBuffer());
-        kernel.setArg(index++, OpenCLContext::ThreadBlockSize*sizeof(mm_float4), NULL);
-        kernel.setArg(index++, OpenCLContext::ThreadBlockSize*sizeof(mm_float4), NULL);
+        donorKernel.setArg<cl::Buffer>(index++, cl.getForceBuffers().getDeviceBuffer());
+        donorKernel.setArg<cl::Buffer>(index++, cl.getEnergyBuffer().getDeviceBuffer());
+        donorKernel.setArg<cl::Buffer>(index++, cl.getPosq().getDeviceBuffer());
+        donorKernel.setArg<cl::Buffer>(index++, donors->getDeviceBuffer());
+        donorKernel.setArg<cl::Buffer>(index++, acceptors->getDeviceBuffer());
+        donorKernel.setArg<cl::Buffer>(index++, donorBufferIndices->getDeviceBuffer());
+        donorKernel.setArg(index++, 3*OpenCLContext::ThreadBlockSize*sizeof(mm_float4), NULL);
         if (globals != NULL)
-            kernel.setArg<cl::Buffer>(index++, globals->getDeviceBuffer());
+            donorKernel.setArg<cl::Buffer>(index++, globals->getDeviceBuffer());
         for (int i = 0; i < (int) donorParams->getBuffers().size(); i++) {
             const OpenCLNonbondedUtilities::ParameterInfo& buffer = donorParams->getBuffers()[i];
-            kernel.setArg<cl::Buffer>(index++, buffer.getBuffer());
+            donorKernel.setArg<cl::Buffer>(index++, buffer.getBuffer());
         }
         for (int i = 0; i < (int) acceptorParams->getBuffers().size(); i++) {
             const OpenCLNonbondedUtilities::ParameterInfo& buffer = acceptorParams->getBuffers()[i];
-            kernel.setArg<cl::Buffer>(index++, buffer.getBuffer());
+            donorKernel.setArg<cl::Buffer>(index++, buffer.getBuffer());
+        }
+        if (tabulatedFunctionParams != NULL) {
+            for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
+                donorKernel.setArg<cl::Buffer>(index++, tabulatedFunctions[i]->getDeviceBuffer());
+            donorKernel.setArg<cl::Buffer>(index++, tabulatedFunctionParams->getDeviceBuffer());
+        }
+        index = 0;
+        acceptorKernel.setArg<cl::Buffer>(index++, cl.getForceBuffers().getDeviceBuffer());
+        acceptorKernel.setArg<cl::Buffer>(index++, cl.getEnergyBuffer().getDeviceBuffer());
+        acceptorKernel.setArg<cl::Buffer>(index++, cl.getPosq().getDeviceBuffer());
+        acceptorKernel.setArg<cl::Buffer>(index++, donors->getDeviceBuffer());
+        acceptorKernel.setArg<cl::Buffer>(index++, acceptors->getDeviceBuffer());
+        acceptorKernel.setArg<cl::Buffer>(index++, acceptorBufferIndices->getDeviceBuffer());
+        acceptorKernel.setArg(index++, 3*OpenCLContext::ThreadBlockSize*sizeof(mm_float4), NULL);
+        if (globals != NULL)
+            acceptorKernel.setArg<cl::Buffer>(index++, globals->getDeviceBuffer());
+        for (int i = 0; i < (int) donorParams->getBuffers().size(); i++) {
+            const OpenCLNonbondedUtilities::ParameterInfo& buffer = donorParams->getBuffers()[i];
+            acceptorKernel.setArg<cl::Buffer>(index++, buffer.getBuffer());
+        }
+        for (int i = 0; i < (int) acceptorParams->getBuffers().size(); i++) {
+            const OpenCLNonbondedUtilities::ParameterInfo& buffer = acceptorParams->getBuffers()[i];
+            acceptorKernel.setArg<cl::Buffer>(index++, buffer.getBuffer());
+        }
+        if (tabulatedFunctionParams != NULL) {
+            for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
+                acceptorKernel.setArg<cl::Buffer>(index++, tabulatedFunctions[i]->getDeviceBuffer());
+            acceptorKernel.setArg<cl::Buffer>(index++, tabulatedFunctionParams->getDeviceBuffer());
         }
     }
-    cl.executeKernel(kernel, std::max(numDonors, numAcceptors));
+    cl.executeKernel(donorKernel, std::max(numDonors, numAcceptors));
+    cl.executeKernel(acceptorKernel, std::max(numDonors, numAcceptors));
 }
 
 double OpenCLCalcCustomHbondForceKernel::executeEnergy(ContextImpl& context) {

platforms/opencl/src/OpenCLKernels.h

@@ -710,7 +710,7 @@ private:
     std::vector<cl_float> globalParamValues;
     std::vector<OpenCLArray<mm_float4>*> tabulatedFunctions;
     System& system;
-    cl::Kernel kernel;
+    cl::Kernel donorKernel, acceptorKernel;
 };
 
 /**

platforms/opencl/src/kernels/customHbondForce.cl

@@ -25,7 +25,6 @@ float4 deltaPeriodic(float4 vec1, float4 vec2) {
 /**
  * Compute the angle between two vectors.  The w component of each vector should contain the squared magnitude.
  */
-
 float computeAngle(float4 vec1, float4 vec2) {
     float dot = vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
     float cosine = dot/sqrt(vec1.w*vec2.w);
@@ -45,120 +44,165 @@ float computeAngle(float4 vec1, float4 vec2) {
 }
 
 /**
- * Compute hbond interactions.
+ * Compute the cross product of two vectors, setting the fourth component to the squared magnitude.
  */
+float4 computeCross(float4 vec1, float4 vec2) {
+    float4 result = cross(vec1, vec2);
+    result.w = result.x*result.x + result.y*result.y + result.z*result.z;
+    return result;
+}
 
-__kernel void computeHbonds(__global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, /*__global unsigned int* exclusions,
-        __global unsigned int* exclusionIndices, */__global int4* donorAtoms, __global int4* acceptorAtoms, __global int4* donorBufferIndices, __global int4* acceptorBufferIndices, __local float4* posBuffer, __local float4* deltaBuffer
+/**
+ * Compute forces on donors.
+ */
+__kernel void computeDonorForces(__global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, /*__global unsigned int* exclusions,
+        __global unsigned int* exclusionIndices, */__global int4* donorAtoms, __global int4* acceptorAtoms, __global int4* donorBufferIndices, __local float4* posBuffer
         PARAMETER_ARGUMENTS) {
     float energy = 0.0f;
-    unsigned int tgx = get_local_id(0) & (get_local_size(0)-1);
-    unsigned int tbx = get_local_id(0) - tgx;
     float4 f1 = 0;
     float4 f2 = 0;
-    for (int donorIndex = get_global_id(0); donorIndex < NUM_DONORS; donorIndex += get_global_size(0)) {
+    float4 f3 = 0;
+    for (int donorStart = 0; donorStart < NUM_DONORS; donorStart += get_global_size(0)) {
         // Load information about the donor this thread will compute forces on.
 
-        int4 atoms = donorAtoms[donorIndex];
-        float4 d1 = posq[atoms.x];
-        float4 d2 = posq[atoms.y];
-        float4 d3 = posq[atoms.z];
-        float4 deltaD1D2 = delta(d1, d2);
+        int donorIndex = donorStart+get_global_id(0);
+        int4 atoms;
+        float4 d1, d2, d3;
+        if (donorIndex < NUM_DONORS) {
+            atoms = donorAtoms[donorIndex];
+            d1 = posq[atoms.x];
+            d2 = posq[atoms.y];
+            d3 = posq[atoms.z];
+        }
+        else
+            atoms = (int4) (-1, -1, -1, -1);
         for (int acceptorStart = 0; acceptorStart < NUM_ACCEPTORS; acceptorStart += get_local_size(0)) {
             // Load the next block of acceptors into local memory.
 
             int blockSize = min((int) get_local_size(0), NUM_ACCEPTORS-acceptorStart);
-            if (tgx < blockSize) {
-                int4 atoms2 = acceptorAtoms[acceptorStart+tgx];
-                float4 pos1 = posq[atoms2.x];
-                float4 pos2 = posq[atoms2.y];
-                float4 pos3 = posq[atoms2.z];
-                posBuffer[get_local_id(0)] = pos1;
-                deltaBuffer[get_local_id(0)] = delta(pos2, pos1);
+            if (get_local_id(0) < blockSize) {
+                int4 atoms2 = acceptorAtoms[acceptorStart+get_local_id(0)];
+                posBuffer[3*get_local_id(0)] = posq[atoms2.x];
+                posBuffer[3*get_local_id(0)+1] = posq[atoms2.y];
+                posBuffer[3*get_local_id(0)+2] = posq[atoms2.z];
             }
             barrier(CLK_LOCAL_MEM_FENCE);
+            if (donorIndex < NUM_DONORS) {
                 for (int index = 0; index < blockSize; index++) {
                     // Compute the interaction between a donor and an acceptor.
 
-                float4 a1 = posBuffer[index];
+                    float4 a1 = posBuffer[3*index];
+                    float4 a2 = posBuffer[3*index+1];
+                    float4 a3 = posBuffer[3*index+2];
                     float4 deltaD1A1 = deltaPeriodic(d1, a1);
 #ifdef USE_CUTOFF
                     if (deltaD1A1.w < CUTOFF_SQUARED) {
 #endif
-                    // Compute variables the force can depend on.
-
-                    float r = sqrt(deltaD1A1.w);
-                    float4 deltaA2A1 = deltaBuffer[index];
-                    float theta = computeAngle(deltaD1A1, deltaD1D2);
-                    float psi = computeAngle(deltaD1A1, deltaA2A1);
-                    float4 cross1 = cross(deltaA2A1, deltaD1A1);
-                    float4 cross2 = cross(deltaD1A1, deltaD1D2);
-                    cross1.w = cross1.x*cross1.x + cross1.y*cross1.y + cross1.z*cross1.z;
-                    cross2.w = cross2.x*cross2.x + cross2.y*cross2.y + cross2.z*cross2.z;
-                    float chi = computeAngle(cross1, cross2);
-                    chi = (dot(deltaA2A1, cross2) < 0 ? -chi : chi);
-                    COMPUTE_FORCE
-
-#ifdef INCLUDE_R
-                    // Apply forces based on r.
-
-                    f1.xyz -= (dEdR/r)*deltaD1A1.xyz;
-#endif
-
-#ifdef INCLUDE_THETA
-                    // Apply forces based on theta.
-
-                    float4 thetaCross = cross(deltaD1D2, deltaD1A1);
-                    float lengthThetaCross = max(length(thetaCross), 1e-6f);
-                    float4 deltaCross0 = cross(deltaD1D2, thetaCross)*dEdTheta/(deltaD1D2.w*lengthThetaCross);
-                    float4 deltaCross2 = -cross(deltaD1A1, thetaCross)*dEdTheta/(deltaD1A1.w*lengthThetaCross);
-                    float4 deltaCross1 = -(deltaCross0+deltaCross2);
-                    f1.xyz += deltaCross1.xyz;
-                    f2.xyz += deltaCross0.xyz;
+                        COMPUTE_DONOR_FORCE
+#ifdef USE_CUTOFF
+                    }
 #endif
+                }
+            }
+        }
 
-#ifdef INCLUDE_PSI
-                    // Apply forces based on psi.
-
-                    float4 psiCross = cross(deltaA2A1, deltaD1A1);
-                    float lengthPsiCross = max(length(psiCross), 1e-6f);
-                    deltaCross0 = cross(deltaD1A1, psiCross)*dEdPsi/(deltaD1A1.w*lengthPsiCross);
-//                    float4 deltaCross2 = -cross(deltaA2A1, psiCross)*dEdPsi/(deltaA2A1.w*lengthPsiCross);
-//                    float4 deltaCross1 = -(deltaCross0+deltaCross2);
-                    f1.xyz += deltaCross0.xyz;
-#endif
+        // Write results
 
-#ifdef INCLUDE_CHI
-                    // Apply forces based on chi.
+        int4 bufferIndices = donorBufferIndices[donorIndex];
+        if (atoms.x > -1) {
+            unsigned int offset = atoms.x+bufferIndices.x*PADDED_NUM_ATOMS;
+            float4 force = forceBuffers[offset];
+            force.xyz += f1.xyz;
+            forceBuffers[offset] = force;
+        }
+        if (atoms.y > -1) {
+            unsigned int offset = atoms.y+bufferIndices.y*PADDED_NUM_ATOMS;
+            float4 force = forceBuffers[offset];
+            force.xyz += f2.xyz;
+            forceBuffers[offset] = force;
+        }
+        if (atoms.z > -1) {
+            unsigned int offset = atoms.z+bufferIndices.z*PADDED_NUM_ATOMS;
+            float4 force = forceBuffers[offset];
+            force.xyz += f3.xyz;
+            forceBuffers[offset] = force;
+        }
+    }
+    energyBuffer[get_global_id(0)] += energy;
+}
+/**
+ * Compute forces on acceptors.
+ */
+__kernel void computeAcceptorForces(__global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, /*__global unsigned int* exclusions,
+        __global unsigned int* exclusionIndices, */__global int4* donorAtoms, __global int4* acceptorAtoms, __global int4* acceptorBufferIndices, __local float4* posBuffer
+        PARAMETER_ARGUMENTS) {
+    float4 f1 = 0;
+    float4 f2 = 0;
+    float4 f3 = 0;
+    for (int acceptorStart = 0; acceptorStart < NUM_ACCEPTORS; acceptorStart += get_global_size(0)) {
+        // Load information about the acceptor this thread will compute forces on.
+
+        int acceptorIndex = acceptorStart+get_global_id(0);
+        int4 atoms;
+        float4 a1, a2, a3;
+        if (acceptorIndex < NUM_ACCEPTORS) {
+            atoms = acceptorAtoms[acceptorIndex];
+            a1 = posq[atoms.x];
+            a2 = posq[atoms.y];
+            a3 = posq[atoms.z];
+        }
+        else
+            atoms = (int4) (-1, -1, -1, -1);
+        for (int donorStart = 0; donorStart < NUM_DONORS; donorStart += get_local_size(0)) {
+            // Load the next block of donors into local memory.
+
+            int blockSize = min((int) get_local_size(0), NUM_DONORS-donorStart);
+            if (get_local_id(0) < blockSize) {
+                int4 atoms2 = donorAtoms[donorStart+get_local_id(0)];
+                posBuffer[3*get_local_id(0)] = posq[atoms2.x];
+                posBuffer[3*get_local_id(0)+1] = posq[atoms2.y];
+                posBuffer[3*get_local_id(0)+2] = posq[atoms2.z];
+            }
+            barrier(CLK_LOCAL_MEM_FENCE);
+            if (acceptorIndex < NUM_ACCEPTORS) {
+                for (int index = 0; index < blockSize; index++) {
+                    // Compute the interaction between a donor and an acceptor.
 
-                    float4 ff;
-                    ff.x = (-dEdChi*r)/cross1.w;
-                    ff.y = (deltaA2A1.x*deltaD1A1.x + deltaA2A1.y*deltaD1A1.y + deltaA2A1.z*deltaD1A1.z)/deltaD1A1.w;
-                    ff.z = (deltaD1D2.x*deltaD1A1.x + deltaD1D2.y*deltaD1A1.y + deltaD1D2.z*deltaD1A1.z)/deltaD1A1.w;
-                    ff.w = (dEdChi*r)/cross2.w;
-                    float4 internalF0 = ff.x*cross1;
-                    float4 internalF3 = ff.w*cross2;
-                    float4 s = ff.y*internalF0 - ff.z*internalF3;
-                    f1.xyz -= s.xyz+internalF3.xyz;
-                    f2.xyz += internalF3.xyz;
+                    float4 d1 = posBuffer[3*index];
+                    float4 d2 = posBuffer[3*index+1];
+                    float4 d3 = posBuffer[3*index+2];
+                    float4 deltaD1A1 = deltaPeriodic(d1, a1);
+#ifdef USE_CUTOFF
+                    if (deltaD1A1.w < CUTOFF_SQUARED) {
 #endif
+                        COMPUTE_ACCEPTOR_FORCE
 #ifdef USE_CUTOFF
                     }
 #endif
                 }
             }
+        }
 
         // Write results
 
-        int4 bufferIndices = donorBufferIndices[donorIndex];
-        unsigned int offset1 = atoms.x+bufferIndices.x*PADDED_NUM_ATOMS;
-        unsigned int offset2 = atoms.y+bufferIndices.y*PADDED_NUM_ATOMS;
-        float4 force1 = forceBuffers[offset1];
-        float4 force2 = forceBuffers[offset2];
-        force1.xyz += f1.xyz;
-        force2.xyz += f2.xyz;
-        forceBuffers[offset1] = force1;
-        forceBuffers[offset2] = force2;
+        int4 bufferIndices = acceptorBufferIndices[acceptorIndex];
+        if (atoms.x > -1) {
+            unsigned int offset = atoms.x+bufferIndices.x*PADDED_NUM_ATOMS;
+            float4 force = forceBuffers[offset];
+            force.xyz += f1.xyz;
+            forceBuffers[offset] = force;
+        }
+        if (atoms.y > -1) {
+            unsigned int offset = atoms.y+bufferIndices.y*PADDED_NUM_ATOMS;
+            float4 force = forceBuffers[offset];
+            force.xyz += f2.xyz;
+            forceBuffers[offset] = force;
+        }
+        if (atoms.z > -1) {
+            unsigned int offset = atoms.z+bufferIndices.z*PADDED_NUM_ATOMS;
+            float4 force = forceBuffers[offset];
+            force.xyz += f3.xyz;
+            forceBuffers[offset] = force;
+        }
     }
-    energyBuffer[get_global_id(0)] += energy;
 }

platforms/opencl/tests/TestOpenCLCustomHbondForce.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-2010 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 the OpenCL implementation of CustomHbondForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "OpenCLPlatform.h"
+#include "openmm/CustomHbondForce.h"
+#include "openmm/HarmonicAngleForce.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/PeriodicTorsionForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testHbond() {
+    OpenCLPlatform platform;
+
+    // Create a system using a CustomHbondForce.
+
+    System customSystem;
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    CustomHbondForce* custom = new CustomHbondForce("0.5*kr*(distance(d1,a1)-r0)^2 + 0.5*ktheta*(angle(a1,d1,d2)-theta0)^2 + 0.5*kpsi*(angle(d1,a1,a2)-psi0)^2 + kchi*(1+cos(n*dihedral(a3,a2,a1,d1)-chi0))");
+    custom->addPerDonorParameter("r0");
+    custom->addPerDonorParameter("theta0");
+    custom->addPerDonorParameter("psi0");
+    custom->addPerAcceptorParameter("chi0");
+    custom->addPerAcceptorParameter("n");
+    custom->addGlobalParameter("kr", 0.4);
+    custom->addGlobalParameter("ktheta", 0.5);
+    custom->addGlobalParameter("kpsi", 0.6);
+    custom->addGlobalParameter("kchi", 0.7);
+    vector<double> parameters(3);
+    parameters[0] = 1.5;
+    parameters[1] = 1.7;
+    parameters[2] = 1.9;
+    custom->addDonor(1, 0, -1, parameters);
+    parameters.resize(2);
+    parameters[0] = 2.1;
+    parameters[1] = 2;
+    custom->addAcceptor(2, 3, 4, parameters);
+    custom->setCutoffDistance(10.0);
+    customSystem.addForce(custom);
+
+    // Create an identical system using HarmonicBondForce, HarmonicAngleForce, and PeriodicTorsionForce.
+
+    System standardSystem;
+    standardSystem.addParticle(1.0);
+    standardSystem.addParticle(1.0);
+    standardSystem.addParticle(1.0);
+    standardSystem.addParticle(1.0);
+    standardSystem.addParticle(1.0);
+    HarmonicBondForce* bond = new HarmonicBondForce();
+    bond->addBond(1, 2, 1.5, 0.4);
+    standardSystem.addForce(bond);
+    HarmonicAngleForce* angle = new HarmonicAngleForce();
+    angle->addAngle(0, 1, 2, 1.7, 0.5);
+    angle->addAngle(1, 2, 3, 1.9, 0.6);
+    standardSystem.addForce(angle);
+    PeriodicTorsionForce* torsion = new PeriodicTorsionForce();
+    torsion->addTorsion(1, 2, 3, 4, 2, 2.1, 0.7);;
+    standardSystem.addForce(torsion);
+
+    // Set the atoms in various positions, and verify that both systems give identical forces and energy.
+
+    init_gen_rand(0);
+    vector<Vec3> positions(5);
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    for (int i = 0; i < 10; i++) {
+        Context c1(customSystem, integrator1, platform);
+        Context c2(standardSystem, integrator2, platform);
+        for (int j = 0; j < (int) positions.size(); j++)
+            positions[j] = Vec3(2.0*genrand_real2(), 2.0*genrand_real2(), 2.0*genrand_real2());
+        c1.setPositions(positions);
+        c2.setPositions(positions);
+        State s1 = c1.getState(State::Forces | State::Energy);
+        State s2 = c2.getState(State::Forces | State::Energy);
+        for (int i = 0; i < customSystem.getNumParticles(); i++)
+            ASSERT_EQUAL_VEC(s2.getForces()[i], s1.getForces()[i], TOL);
+        ASSERT_EQUAL_TOL(s2.getPotentialEnergy(), s1.getPotentialEnergy(), TOL);
+    }
+}
+
+void testExclusions() {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomHbondForce* custom = new CustomHbondForce("(distance(d1,a1)-1)^2");
+    custom->addDonor(0, 1, -1, vector<double>());
+    custom->addDonor(1, 0, -1, vector<double>());
+    custom->addAcceptor(2, 0, -1, vector<double>());
+    custom->addExclusion(1, 0);
+    system.addForce(custom);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0, 2, 0);
+    positions[2] = Vec3(2, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(2, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-2, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(1.0, state.getPotentialEnergy(), TOL);
+}
+
+void testCutoff() {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomHbondForce* custom = new CustomHbondForce("(distance(d1,a1)-1)^2");
+    custom->addDonor(0, 1, -1, vector<double>());
+    custom->addDonor(1, 0, -1, vector<double>());
+    custom->addAcceptor(2, 0, -1, vector<double>());
+    custom->setNonbondedMethod(CustomHbondForce::CutoffNonPeriodic);
+    custom->setCutoffDistance(2.5);
+    system.addForce(custom);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0, 3, 0);
+    positions[2] = Vec3(2, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(2, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-2, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(1.0, state.getPotentialEnergy(), TOL);
+}
+
+void testCustomFunctions() {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomHbondForce* custom = new CustomHbondForce("foo(distance(d1,a1))");
+    custom->addDonor(1, 0, -1, vector<double>());
+    custom->addDonor(2, 0, -1, vector<double>());
+    custom->addAcceptor(0, 1, -1, vector<double>());
+    vector<double> function(2);
+    function[0] = 0;
+    function[1] = 1;
+    custom->addFunction("foo", function, 0, 10, true);
+    system.addForce(custom);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0, 2, 0);
+    positions[2] = Vec3(2, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(0.1, 0.1, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -0.1, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-0.1, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(0.1*2+0.1*2, state.getPotentialEnergy(), TOL);
+}
+
+int main() {
+    try {
+        testHbond();
+//        testExclusions();
+        testCutoff();
+        testCustomFunctions();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+

platforms/reference/src/ReferenceKernels.cpp

@@ -1314,7 +1314,7 @@ void ReferenceCalcCustomHbondForceKernel::initialize(const System& system, const
     map<string, vector<int> > distances;
     map<string, vector<int> > angles;
     map<string, vector<int> > dihedrals;
-    Lepton::ParsedExpression energyExpression = CustomHbondForceImpl::prepareExpression(force, distances, angles, dihedrals);
+    Lepton::ParsedExpression energyExpression = CustomHbondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
     vector<string> donorParameterNames;
     vector<string> acceptorParameterNames;
     for (int i = 0; i < numDonorParameters; i++)

platforms/reference/tests/TestReferenceCustomHbondForce.cpp

@@ -179,11 +179,42 @@ void testCutoff() {
     ASSERT_EQUAL_TOL(1.0, state.getPotentialEnergy(), TOL);
 }
 
+void testCustomFunctions() {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomHbondForce* custom = new CustomHbondForce("foo(distance(d1,a1))");
+    custom->addDonor(1, 0, -1, vector<double>());
+    custom->addDonor(2, 0, -1, vector<double>());
+    custom->addAcceptor(0, 1, -1, vector<double>());
+    vector<double> function(2);
+    function[0] = 0;
+    function[1] = 1;
+    custom->addFunction("foo", function, 0, 10, true);
+    system.addForce(custom);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0, 2, 0);
+    positions[2] = Vec3(2, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(0.1, 0.1, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -0.1, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-0.1, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(0.1*2+0.1*2, state.getPotentialEnergy(), TOL);
+}
+
 int main() {
     try {
         testHbond();
         testExclusions();
         testCutoff();
+        testCustomFunctions();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;