Began OpenCL implementation of CustomNonbondedForce

platforms/opencl/src/OpenCLArray.h

@@ -51,11 +51,12 @@ public:
      * @param name              the name of the array
      * @param createHostBuffer  specifies whether to create a buffer in host memory for copying data to and from
      *                          the OpenCL Buffer
+     * @param flags             the set of flags to specify when creating the OpenCL Buffer
      */
-    OpenCLArray(OpenCLContext& context, int size, const std::string& name, bool createHostBuffer = false) :
+    OpenCLArray(OpenCLContext& context, int size, const std::string& name, bool createHostBuffer = false, cl_int flags = CL_MEM_READ_WRITE) :
             context(context), size(size), name(name), local(createHostBuffer ? size : 0), ownsBuffer(true) {
         try {
-            buffer = new cl::Buffer(context.getContext(), CL_MEM_READ_WRITE, size*sizeof(T));
+            buffer = new cl::Buffer(context.getContext(), flags, size*sizeof(T));
         }
         catch (cl::Error err) {
             std::stringstream str;

platforms/opencl/src/OpenCLExpressionUtilities.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) 2009 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "OpenCLExpressionUtilities.h"
+#include "openmm/OpenMMException.h"
+#include "lepton/Operation.h"
+#include <sstream>
+
+using namespace OpenMM;
+using namespace Lepton;
+using namespace std;
+
+static string doubleToString(double value) {
+    stringstream s;
+    s.precision(8);
+    s << scientific << value << "f";
+    return s.str();
+}
+
+static string intToString(int value) {
+    stringstream s;
+    s << value;
+    return s.str();
+}
+
+string OpenCLExpressionUtilities::createExpression(const ParsedExpression& expression, const map<string, string>& variables) {
+    return processExpression(expression.getRootNode(), variables);
+}
+
+string OpenCLExpressionUtilities::processExpression(const ExpressionTreeNode& node, const map<string, string>& variables) {
+    switch (node.getOperation().getId()) {
+        case Operation::CONSTANT:
+            return doubleToString(dynamic_cast<const Operation::Constant*>(&node.getOperation())->getValue());
+        case Operation::VARIABLE:
+        {
+            map<string, string>::const_iterator iter = variables.find(node.getOperation().getName());
+            if (iter == variables.end())
+                throw OpenMMException("Unknown variable in expression: "+node.getOperation().getName());
+            return iter->second;
+        }
+        case Operation::ADD:
+             return "("+processExpression(node.getChildren()[0], variables)+")+("+processExpression(node.getChildren()[1], variables)+")";
+        case Operation::SUBTRACT:
+            return "("+processExpression(node.getChildren()[0], variables)+")-("+processExpression(node.getChildren()[1], variables)+")";
+        case Operation::MULTIPLY:
+            return "("+processExpression(node.getChildren()[0], variables)+")*("+processExpression(node.getChildren()[1], variables)+")";
+        case Operation::DIVIDE:
+            return "("+processExpression(node.getChildren()[0], variables)+")/("+processExpression(node.getChildren()[1], variables)+")";
+        case Operation::POWER:
+            return "pow(("+processExpression(node.getChildren()[0], variables)+"), ("+processExpression(node.getChildren()[1], variables)+"))";
+        case Operation::NEGATE:
+            return "-("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::SQRT:
+            return "sqrt("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::EXP:
+            return "exp("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::LOG:
+            return "log("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::SIN:
+            return "sin("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::COS:
+            return "cos("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::SEC:
+            return "1.0f/cos("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::CSC:
+            return "1.0f/sin("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::TAN:
+            return "tan("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::COT:
+            return "1.0f/tan("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::ASIN:
+            return "asin("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::ACOS:
+            return "acos("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::ATAN:
+            return "atan("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::SQUARE:
+            return "pow(("+processExpression(node.getChildren()[0], variables)+"), 2.0f)";
+        case Operation::CUBE:
+            return "pow(("+processExpression(node.getChildren()[0], variables)+"), 3.0f)";
+        case Operation::RECIPROCAL:
+            return "1.0f/("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::ADD_CONSTANT:
+            return doubleToString(dynamic_cast<const Operation::AddConstant*>(&node.getOperation())->getValue())+"+("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::MULTIPLY_CONSTANT:
+            return doubleToString(dynamic_cast<const Operation::MultiplyConstant*>(&node.getOperation())->getValue())+"*("+processExpression(node.getChildren()[0], variables)+")";
+        case Operation::POWER_CONSTANT:
+            return "pow(("+processExpression(node.getChildren()[0], variables)+"), "+doubleToString(dynamic_cast<const Operation::PowerConstant*>(&node.getOperation())->getValue())+")";
+    }
+    throw OpenMMException("Internal error: Unknown operation in user-defined expression: "+node.getOperation().getName());
+}

platforms/opencl/src/OpenCLExpressionUtilities.h

+#ifndef OPENMM_OPENCLEXPRESSIONUTILITIES_H_
+#define OPENMM_OPENCLEXPRESSIONUTILITIES_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   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) 2009 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "lepton/ExpressionTreeNode.h"
+#include "lepton/ParsedExpression.h"
+#include <map>
+#include <string>
+
+namespace OpenMM {
+
+/**
+ * This class is used by various classes to generate OpenCL source code implementing
+ * user defined mathematical expressions.
+ */
+
+class OpenCLExpressionUtilities {
+public:
+    static std::string createExpression(const Lepton::ParsedExpression& expression, const std::map<std::string, std::string>& variables);
+private:
+    static std::string processExpression(const Lepton::ExpressionTreeNode& node, const std::map<std::string, std::string>& variables);
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_OPENCLEXPRESSIONUTILITIES_H_*/

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -48,8 +48,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
         return new OpenCLCalcRBTorsionForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcNonbondedForceKernel::Name())
         return new OpenCLCalcNonbondedForceKernel(name, platform, cl, context.getSystem());
-//    if (name == CalcCustomNonbondedForceKernel::Name())
-//        return new OpenCLCalcCustomNonbondedForceKernel(name, platform, cl, context.getSystem());
+    if (name == CalcCustomNonbondedForceKernel::Name())
+        return new OpenCLCalcCustomNonbondedForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcGBSAOBCForceKernel::Name())
         return new OpenCLCalcGBSAOBCForceKernel(name, platform, cl);
     if (name == IntegrateVerletStepKernel::Name())

platforms/opencl/src/OpenCLKernels.cpp

@@ -30,8 +30,11 @@
 #include "openmm/Context.h"
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
+#include "OpenCLExpressionUtilities.h"
 #include "OpenCLIntegrationUtilities.h"
 #include "OpenCLNonbondedUtilities.h"
+#include "lepton/Parser.h"
+#include "lepton/ParsedExpression.h"
 #include <cmath>
 
 using namespace OpenMM;
@@ -666,122 +669,254 @@ double OpenCLCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
     return ewaldSelfEnergy;
 }
 
-//OpenCLCalcCustomNonbondedForceKernel::~OpenCLCalcCustomNonbondedForceKernel() {
-//}
-//
-//void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
-//    data.primaryKernel = this; // This must always be the primary kernel so it can update the global parameters
-//    data.hasCustomNonbonded = true;
-//    numParticles = force.getNumParticles();
-//    _gpuContext* gpu = data.gpu;
-//
-//    // Identify which exceptions are actual interactions.
-//
-//    vector<pair<int, int> > exclusions;
-//    vector<int> exceptions;
-//    {
-//        vector<double> parameters;
-//        for (int i = 0; i < force.getNumExceptions(); i++) {
-//            int particle1, particle2;
-//            force.getExceptionParameters(i, particle1, particle2, parameters);
-//            exclusions.push_back(pair<int, int>(particle1, particle2));
-//            if (parameters.size() > 0)
-//                exceptions.push_back(i);
-//        }
-//    }
-//
-//    // Initialize nonbonded interactions.
-//
-//    vector<int> particle(numParticles);
-//    vector<vector<double> > parameters(numParticles);
-//    vector<vector<int> > exclusionList(numParticles);
-//    for (int i = 0; i < numParticles; i++) {
-//        force.getParticleParameters(i, parameters[i]);
-//        particle[i] = i;
-//        exclusionList[i].push_back(i);
-//    }
-//    for (int i = 0; i < (int)exclusions.size(); i++) {
-//        exclusionList[exclusions[i].first].push_back(exclusions[i].second);
-//        exclusionList[exclusions[i].second].push_back(exclusions[i].first);
-//    }
-//    Vec3 boxVectors[3];
-//    system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
-//    gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
-//    OpenCLNonbondedMethod method = NO_CUTOFF;
-//    if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff)
-//        method = CUTOFF;
-//    if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic) {
-//        method = PERIODIC;
-//    }
-//    data.customNonbondedMethod = method;
-//
-//    // Initialize exceptions.
-//
-//    int numExceptions = exceptions.size();
-//    vector<int> exceptionParticle1(numExceptions);
-//    vector<int> exceptionParticle2(numExceptions);
-//    vector<vector<double> > exceptionParams(numExceptions);
-//    for (int i = 0; i < numExceptions; i++)
-//        force.getExceptionParameters(exceptions[i], exceptionParticle1[i], exceptionParticle2[i], exceptionParams[i]);
-//
-//    // Record the tabulated functions.
-//
-//    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);
+class OpenCLCustomNonbondedForceInfo : public OpenCLForceInfo {
+public:
+    OpenCLCustomNonbondedForceInfo(int requiredBuffers, const CustomNonbondedForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
+    }
+    bool areParticlesIdentical(int particle1, int particle2) {
+        vector<double> params1;
+        vector<double> params2;
+        force.getParticleParameters(particle1, params1);
+        force.getParticleParameters(particle2, params2);
+        for (int i = 0; i < params1.size(); i++)
+            if (params1[i] != params2[i])
+                return false;
+        return true;
+    }
+    int getNumParticleGroups() {
+        return force.getNumExceptions();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2;
+        vector<double> params;
+        force.getExceptionParameters(index, particle1, particle2, params);
+        particles.resize(2);
+        particles[0] = particle1;
+        particles[1] = particle2;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2;
+        vector<double> params1;
+        vector<double> params2;
+        force.getExceptionParameters(group1, particle1, particle2, params1);
+        force.getExceptionParameters(group2, particle1, particle2, params2);
+        for (int i = 0; i < params1.size(); i++)
+            if (params1[i] != params2[i])
+                return false;
+        return true;
+    }
+private:
+    const CustomNonbondedForce& force;
+};
+
+OpenCLCalcCustomNonbondedForceKernel::~OpenCLCalcCustomNonbondedForceKernel() {
+    if (params != NULL)
+        delete params;
+    if (globals != NULL)
+        delete globals;
+    if (exceptionParams != NULL)
+        delete exceptionParams;
+    if (exceptionIndices != NULL)
+        delete exceptionIndices;
+}
+
+void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
+    if (force.getNumParameters() > 4)
+        throw OpenMMException("OpenCLPlatform only supports four per-atom parameters for custom nonbonded forces");
+    int forceIndex;
+    for (forceIndex = 0; forceIndex < system.getNumForces() && &system.getForce(forceIndex) != &force; ++forceIndex)
+        ;
+    string prefix = "custom"+intToString(forceIndex)+"_";
+
+    // Identify which exceptions are actual interactions.
+
+    vector<pair<int, int> > exclusions;
+    vector<int> exceptions;
+    {
+        vector<double> parameters;
+        for (int i = 0; i < force.getNumExceptions(); i++) {
+            int particle1, particle2;
+            force.getExceptionParameters(i, particle1, particle2, parameters);
+            exclusions.push_back(pair<int, int>(particle1, particle2));
+            if (parameters.size() > 0)
+                exceptions.push_back(i);
+        }
+    }
+
+    // Record parameters and exclusions.
+
+    int numParticles = force.getNumParticles();
+    params = new OpenCLArray<mm_float4>(cl, numParticles, "customNonbondedParameters");
+    if (force.getNumGlobalParameters() > 0)
+        globals = new OpenCLArray<cl_float>(cl, force.getNumGlobalParameters(), "customNonbondedGlobals", false, CL_MEM_READ_ONLY);
+    vector<mm_float4> paramVec(numParticles);
+    vector<vector<int> > exclusionList(numParticles);
+    for (int i = 0; i < numParticles; i++) {
+        vector<double> parameters;
+        force.getParticleParameters(i, parameters);
+        if (parameters.size() > 0)
+            paramVec[i].x = (cl_float) parameters[0];
+        if (parameters.size() > 1)
+            paramVec[i].y = (cl_float) parameters[1];
+        if (parameters.size() > 2)
+            paramVec[i].z = (cl_float) parameters[2];
+        if (parameters.size() > 3)
+            paramVec[i].w = (cl_float) parameters[3];
+        exclusionList[i].push_back(i);
+    }
+    for (int i = 0; i < (int)exclusions.size(); i++) {
+        exclusionList[exclusions[i].first].push_back(exclusions[i].second);
+        exclusionList[exclusions[i].second].push_back(exclusions[i].first);
+    }
+    params->upload(paramVec);
+
+    // Record the tabulated functions.
+
+    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);
 //        gpuSetTabulatedFunction(gpu, i, name, values, min, max, interpolating);
-//    }
-//
-//    // Record information for the expressions.
-//
-//    vector<string> paramNames;
-//    vector<string> combiningRules;
-//    for (int i = 0; i < force.getNumParameters(); i++) {
-//        paramNames.push_back(force.getParameterName(i));
-//        combiningRules.push_back(force.getParameterCombiningRule(i));
-//    }
-//    globalParamNames.resize(force.getNumGlobalParameters());
-//    globalParamValues.resize(force.getNumGlobalParameters());
-//    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
-//        globalParamNames[i] = force.getGlobalParameterName(i);
-//        globalParamValues[i] = force.getGlobalParameterDefaultValue(i);
-//    }
-//    gpuSetCustomNonbondedParameters(gpu, parameters, exclusionList, exceptionParticle1, exceptionParticle2, exceptionParams, method,
-//            (float)force.getCutoffDistance(), force.getEnergyFunction(), combiningRules, paramNames, globalParamNames);
-//    if (globalParamValues.size() > 0)
-//        SetCustomNonbondedGlobalParams(&globalParamValues[0]);
-//}
-//
-//void OpenCLCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
-//    if (data.primaryKernel == this) {
-//        updateGlobalParams(context);
-//        calcForces(context, data);
-//    }
-//}
-//
-//double OpenCLCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
-//    if (data.primaryKernel == this) {
-//        updateGlobalParams(context);
-//        return calcEnergy(context, data, system);
-//    }
-//    return 0.0;
-//}
-//
-//void OpenCLCalcCustomNonbondedForceKernel::updateGlobalParams(ContextImpl& context) {
-//    bool changed = false;
-//    for (int i = 0; i < globalParamNames.size(); i++) {
-//        float value = (float) context.getParameter(globalParamNames[i]);
-//        if (value != globalParamValues[i])
-//            changed = true;
-//        globalParamValues[i] = value;
-//    }
-//    if (changed)
-//        SetCustomNonbondedGlobalParams(&globalParamValues[0]);
-//}
-//
+    }
+
+    // Record information for the expressions.
+
+    vector<string> paramNames;
+    vector<string> combiningRules;
+    for (int i = 0; i < force.getNumParameters(); i++) {
+        paramNames.push_back(force.getParameterName(i));
+        combiningRules.push_back(force.getParameterCombiningRule(i));
+    }
+    globalParamNames.resize(force.getNumGlobalParameters());
+    globalParamValues.resize(force.getNumGlobalParameters());
+    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+        globalParamNames[i] = force.getGlobalParameterName(i);
+        globalParamValues[i] = (cl_float) force.getGlobalParameterDefaultValue(i);
+    }
+    if (globals != NULL)
+        globals->upload(globalParamValues);
+    bool useCutoff = (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff);
+    bool usePeriodic = (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff && force.getNonbondedMethod() != CustomNonbondedForce::CutoffNonPeriodic);
+    Lepton::ParsedExpression energyExpression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
+    Lepton::ParsedExpression forceExpression = energyExpression.differentiate("r").optimize();
+
+    // Create the kernels.
+
+    map<string, string> paramVariables;
+    map<string, string> forceVariables;
+    map<string, string> exceptionVariables;
+    forceVariables["r"] = "r";
+    exceptionVariables["r"] = "r";
+    string suffixes[] = {".x", ".y", ".z", ".w"};
+    for (int i = 0; i < force.getNumParameters(); i++) {
+        const string& name = force.getParameterName(i);
+        paramVariables[name+"1"] = prefix+"params1"+suffixes[i];
+        paramVariables[name+"2"] = prefix+"params2"+suffixes[i];
+        forceVariables[name] = prefix+name;
+        exceptionVariables[name] = "exceptionParams"+suffixes[i];
+    }
+    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+        const string& name = force.getGlobalParameterName(i);
+        string value = "globals["+intToString(i)+"]";
+        paramVariables[name] = prefix+value;
+        forceVariables[name] = prefix+value;
+        exceptionVariables[name] = value;
+    }
+    stringstream compute;
+    for (int i = 0; i < force.getNumParameters(); i++) {
+        Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getParameterCombiningRule(i)).optimize();
+        compute << "float " << prefix << force.getParameterName(i) << " = " << OpenCLExpressionUtilities::createExpression(expression, paramVariables) << ";\n";
+    }
+    compute << "tempEnergy += " << OpenCLExpressionUtilities::createExpression(energyExpression, forceVariables) << ";\n";
+    compute << "tempForce -= " << OpenCLExpressionUtilities::createExpression(forceExpression, forceVariables) << ";\n";
+    map<string, string> replacements;
+    replacements["COMPUTE_FORCE"] = compute.str();
+    string source = cl.loadSourceFromFile("customNonbonded.cl", replacements);
+    cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source);
+    cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo(prefix+"params", "float4", sizeof(cl_float4), params->getDeviceBuffer()));
+    if (globals != NULL) {
+        globals->upload(globalParamValues);
+        cl.getNonbondedUtilities().addArgument(OpenCLNonbondedUtilities::ParameterInfo(prefix+"globals", "float", sizeof(cl_float), globals->getDeviceBuffer()));
+    }
+    stringstream computeExceptions;
+    computeExceptions << "energy += " << OpenCLExpressionUtilities::createExpression(energyExpression, exceptionVariables) << ";\n";
+    computeExceptions << "dEdR = " << OpenCLExpressionUtilities::createExpression(forceExpression, exceptionVariables) << ";\n";
+    replacements["COMPUTE_FORCE"] = computeExceptions.str();
+    map<string, string> defines;
+    if (globals != NULL)
+        defines["HAS_GLOBALS"] = "1";
+    cl::Program program = cl.createProgram(cl.loadSourceFromFile("customNonbondedExceptions.cl", replacements), defines);
+    exceptionsKernel = cl::Kernel(program, "computeCustomNonbondedExceptions");
+
+    // Initialize exception parameters.
+
+    int numExceptions = exceptions.size();
+    int maxBuffers = cl.getNonbondedUtilities().getNumForceBuffers();
+    if (numExceptions > 0) {
+        exceptionParams = new OpenCLArray<mm_float4>(cl, numExceptions, "customExceptionParams");
+        exceptionIndices = new OpenCLArray<mm_int4>(cl, numExceptions, "customExceptionIndices");
+        vector<mm_float4> exceptionParamsVector(numExceptions);
+        vector<mm_int4> exceptionIndicesVector(numExceptions);
+        vector<int> forceBufferCounter(system.getNumParticles(), 0);
+        for (int i = 0; i < numExceptions; i++) {
+            int particle1, particle2;
+            vector<double> parameters;
+            force.getExceptionParameters(exceptions[i], particle1, particle2, parameters);
+            if (parameters.size() > 0)
+                exceptionParamsVector[i].x = (cl_float) parameters[0];
+            if (parameters.size() > 1)
+                exceptionParamsVector[i].y = (cl_float) parameters[1];
+            if (parameters.size() > 2)
+                exceptionParamsVector[i].z = (cl_float) parameters[2];
+            if (parameters.size() > 3)
+                exceptionParamsVector[i].w = (cl_float) parameters[3];
+            exceptionIndicesVector[i] = (mm_int4) {particle1, particle2, forceBufferCounter[particle1]++, forceBufferCounter[particle2]++};
+        }
+        exceptionParams->upload(exceptionParamsVector);
+        exceptionIndices->upload(exceptionIndicesVector);
+        for (int i = 0; i < forceBufferCounter.size(); i++)
+            maxBuffers = max(maxBuffers, forceBufferCounter[i]);
+    }
+    cl.addForce(new OpenCLCustomNonbondedForceInfo(maxBuffers, force));
+}
+
+void OpenCLCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
+    if (exceptionParams != NULL) {
+        if (!hasCreatedKernels) {
+            hasCreatedKernels = true;
+            exceptionsKernel.setArg<cl_int>(0, cl.getPaddedNumAtoms());
+            exceptionsKernel.setArg<cl_int>(1, exceptionParams->getSize());
+            exceptionsKernel.setArg<cl_float>(2, cl.getNonbondedUtilities().getCutoffDistance()*cl.getNonbondedUtilities().getCutoffDistance());
+            exceptionsKernel.setArg<mm_float4>(3, cl.getNonbondedUtilities().getPeriodicBoxSize());
+            exceptionsKernel.setArg<cl::Buffer>(4, cl.getForceBuffers().getDeviceBuffer());
+            exceptionsKernel.setArg<cl::Buffer>(5, cl.getEnergyBuffer().getDeviceBuffer());
+            exceptionsKernel.setArg<cl::Buffer>(6, cl.getPosq().getDeviceBuffer());
+            exceptionsKernel.setArg<cl::Buffer>(7, exceptionParams->getDeviceBuffer());
+            exceptionsKernel.setArg<cl::Buffer>(8, exceptionIndices->getDeviceBuffer());
+            if (globals != NULL)
+                exceptionsKernel.setArg<cl::Buffer>(9, globals->getDeviceBuffer());
+        }
+        cl.executeKernel(exceptionsKernel, exceptionIndices->getSize());
+    }
+    if (globals == NULL)
+        return;
+    bool changed = false;
+    for (int i = 0; i < globalParamNames.size(); i++) {
+        cl_float value = (cl_float) context.getParameter(globalParamNames[i]);
+        if (value != globalParamValues[i])
+            changed = true;
+        globalParamValues[i] = value;
+    }
+    if (changed)
+        globals->upload(globalParamValues);
+}
+
+double OpenCLCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
+    executeForces(context);
+    return 0.0;
+}
 
 class OpenCLGBSAOBCForceInfo : public OpenCLForceInfo {
 public:

platforms/opencl/src/OpenCLKernels.h

@@ -335,42 +335,47 @@ private:
     double cutoffSquared, ewaldSelfEnergy;
 };
 
-///**
-// * This kernel is invoked by CustomNonbondedForce to calculate the forces acting on the system.
-// */
-//class OpenCLCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
-//public:
-//    OpenCLCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomNonbondedForceKernel(name, platform), cl(cl), system(system) {
-//    }
-//    ~OpenCLCalcCustomNonbondedForceKernel();
-//    /**
-//     * Initialize the kernel.
-//     *
-//     * @param system     the System this kernel will be applied to
-//     * @param force      the CustomNonbondedForce this kernel will be used for
-//     */
-//    void initialize(const System& system, const CustomNonbondedForce& force);
-//    /**
-//     * Execute the kernel to calculate the forces.
-//     *
-//     * @param context    the context in which to execute this kernel
-//     */
-//    void executeForces(ContextImpl& context);
-//    /**
-//     * Execute the kernel to calculate the energy.
-//     *
-//     * @param context    the context in which to execute this kernel
-//     * @return the potential energy due to the CustomNonbondedForce
-//     */
-//    double executeEnergy(ContextImpl& context);
-//private:
-//    void updateGlobalParams(ContextImpl& context);
-//    OpenCLContext& cl;
-//    int numParticles;
-//    std::vector<std::string> globalParamNames;
-//    std::vector<float> globalParamValues;
-//    System& system;
-//};
+/**
+ * This kernel is invoked by CustomNonbondedForce to calculate the forces acting on the system.
+ */
+class OpenCLCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
+public:
+    OpenCLCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomNonbondedForceKernel(name, platform),
+            hasCreatedKernels(false), cl(cl), params(NULL), globals(NULL), exceptionParams(NULL), exceptionIndices(NULL), system(system) {
+    }
+    ~OpenCLCalcCustomNonbondedForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomNonbondedForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomNonbondedForce& force);
+    /**
+     * Execute the kernel to calculate the forces.
+     *
+     * @param context    the context in which to execute this kernel
+     */
+    void executeForces(ContextImpl& context);
+    /**
+     * Execute the kernel to calculate the energy.
+     *
+     * @param context    the context in which to execute this kernel
+     * @return the potential energy due to the CustomNonbondedForce
+     */
+    double executeEnergy(ContextImpl& context);
+private:
+    bool hasCreatedKernels;
+    OpenCLContext& cl;
+    OpenCLArray<mm_float4>* params;
+    OpenCLArray<cl_float>* globals;
+    OpenCLArray<mm_float4>* exceptionParams;
+    OpenCLArray<mm_int4>* exceptionIndices;
+    cl::Kernel exceptionsKernel;
+    std::vector<std::string> globalParamNames;
+    std::vector<cl_float> globalParamValues;
+    System& system;
+};
 
 /**
  * This kernel is invoked by GBSAOBCForce to calculate the forces acting on the system.

platforms/opencl/src/OpenCLNonbondedUtilities.cpp

@@ -103,6 +103,10 @@ void OpenCLNonbondedUtilities::addParameter(const ParameterInfo& parameter) {
     parameters.push_back(parameter);
 }
 
+void OpenCLNonbondedUtilities::addArgument(const ParameterInfo& parameter) {
+    arguments.push_back(parameter);
+}
+
 void OpenCLNonbondedUtilities::initialize(const System& system) {
     if (cutoff == -1.0)
         return; // There are no nonbonded interactions in the System.
@@ -222,7 +226,7 @@ void OpenCLNonbondedUtilities::initialize(const System& system) {
 
     // Create kernels.
 
-    forceKernel = createInteractionKernel(kernelSource, parameters, true);
+    forceKernel = createInteractionKernel(kernelSource, parameters, arguments, true);
     if (useCutoff) {
         map<string, string> defines;
         if (forceBufferPerAtomBlock)
@@ -274,7 +278,7 @@ void OpenCLNonbondedUtilities::computeInteractions() {
         context.executeKernel(forceKernel, tiles->getSize()*OpenCLContext::TileSize);
 }
 
-cl::Kernel OpenCLNonbondedUtilities::createInteractionKernel(const string& source, const vector<ParameterInfo> params, bool useExclusions) const {
+cl::Kernel OpenCLNonbondedUtilities::createInteractionKernel(const string& source, const vector<ParameterInfo>& params, const vector<ParameterInfo>& arguments, bool useExclusions) const {
     map<string, string> replacements;
     replacements["COMPUTE_INTERACTION"] = source;
     stringstream args;
@@ -288,6 +292,15 @@ cl::Kernel OpenCLNonbondedUtilities::createInteractionKernel(const string& sourc
         args << "* local_";
         args << params[i].getName();
     }
+    for (int i = 0; i < arguments.size(); i++) {
+        if ((arguments[i].getBuffer().getInfo<CL_MEM_FLAGS>() & CL_MEM_READ_ONLY) == 0)
+            args << ", __global ";
+        else
+            args << ", __constant ";
+        args << arguments[i].getType();
+        args << "* ";
+        args << arguments[i].getName();
+    }
     replacements["PARAMETER_ARGUMENTS"] = args.str();
     stringstream loadLocal1;
     for (int i = 0; i < params.size(); i++) {
@@ -383,5 +396,9 @@ cl::Kernel OpenCLNonbondedUtilities::createInteractionKernel(const string& sourc
         kernel.setArg<cl::Buffer>(i*2+paramBase, params[i].getBuffer());
         kernel.setArg(i*2+paramBase+1, OpenCLContext::ThreadBlockSize*params[i].getSize(), NULL);
     }
+    paramBase += 2*params.size();
+    for (int i = 0; i < (int) arguments.size(); i++) {
+        kernel.setArg<cl::Buffer>(i+paramBase, arguments[i].getBuffer());
+    }
     return kernel;
 }

platforms/opencl/src/OpenCLNonbondedUtilities.h

@@ -81,6 +81,10 @@ public:
      * Add a per-atom parameter that the default interaction kernel may depend on.
      */
     void addParameter(const ParameterInfo& parameter);
+    /**
+     * Add an array (other than a per-atom parameter) that should be passed as an argument to the default interaction kernel.
+     */
+    void addArgument(const ParameterInfo& parameter);
     /**
      * Initialize this object in preparation for a simulation.
      */
@@ -173,9 +177,10 @@ public:
      * 
      * @param source        the source code for evaluating the force and energy
      * @param params        the per-atom parameters this kernel may depend on
+     * @param arguments     arrays (other than per-atom parameters) that should be passed as arguments to the kernel
      * @param useExclusions specifies whether exclusions are applied to this interaction
      */
-    cl::Kernel createInteractionKernel(const std::string& source, const std::vector<ParameterInfo> params, bool useExclusions) const;
+    cl::Kernel createInteractionKernel(const std::string& source, const std::vector<ParameterInfo>& params, const std::vector<ParameterInfo>& arguments, bool useExclusions) const;
 private:
     OpenCLContext& context;
     cl::Kernel forceKernel;
@@ -192,6 +197,7 @@ private:
     OpenCLArray<mm_float4>* blockBoundingBox;
     std::vector<std::vector<int> > atomExclusions;
     std::vector<ParameterInfo> parameters;
+    std::vector<ParameterInfo> arguments;
     OpenCLCompact* compact;
     std::string kernelSource;
     std::map<std::string, std::string> kernelDefines;

platforms/opencl/src/OpenCLPlatform.cpp

@@ -52,7 +52,7 @@ OpenCLPlatform::OpenCLPlatform() {
     registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
-//    registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory);
     registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);

platforms/opencl/src/kernels/customNonbonded.cl

+#ifdef USE_CUTOFF
+if (!isExcluded && r2 < CUTOFF_SQUARED) {
+#else
+if (!isExcluded) {
+#endif
+    float tempForce = 0.0f;
+    COMPUTE_FORCE
+    dEdR += tempForce*invR;
+}

platforms/opencl/src/kernels/customNonbondedExceptions.cl

+/**
+ * Compute custom nonbonded exceptions.
+ */
+
+__kernel void computeCustomNonbondedExceptions(int numAtoms, int numExceptions, float cutoffSquared, float4 periodicBoxSize, __global float4* forceBuffers, __global float* energyBuffer,
+        __global float4* posq, __global float4* params, __global int4* indices
+#ifdef HAS_GLOBALS
+        , __constant float* globals) {
+#else
+        ) {
+#endif
+    int index = get_global_id(0);
+    float energy = 0.0f;
+    while (index < numExceptions) {
+        // Look up the data for this exception.
+
+        int4 atoms = indices[index];
+        float4 exceptionParams = params[index];
+        float4 delta = posq[atoms.y]-posq[atoms.x];
+#ifdef USE_PERIODIC
+        delta.x -= floor(delta.x/periodicBoxSize.x+0.5f)*periodicBoxSize.x;
+        delta.y -= floor(delta.y/periodicBoxSize.y+0.5f)*periodicBoxSize.y;
+        delta.z -= floor(delta.z/periodicBoxSize.z+0.5f)*periodicBoxSize.z;
+#endif
+        // Compute the force.
+
+        float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
+#ifdef USE_CUTOFF
+        if (r2 > cutoffSquared) {
+#else
+        {
+#endif
+            float r = sqrt(r2);
+            float dEdR;
+            COMPUTE_FORCE
+            delta.xyz *= -dEdR/r;
+
+            // Record the force on each of the two atoms.
+
+            unsigned int offsetA = atoms.x+atoms.z*numAtoms;
+            unsigned int offsetB = atoms.y+atoms.w*numAtoms;
+            float4 forceA = forceBuffers[offsetA];
+            float4 forceB = forceBuffers[offsetB];
+            forceA.xyz -= delta.xyz;
+            forceB.xyz += delta.xyz;
+            forceBuffers[offsetA] = forceA;
+            forceBuffers[offsetB] = forceB;
+        }
+        index += get_global_size(0);
+    }
+    energyBuffer[get_global_id(0)] += energy;
+}

platforms/opencl/tests/TestOpenCLCustomNonbondedForce.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-2009 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 OpenCL implementation of CustomNonbondedForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "OpenCLPlatform.h"
+#include "openmm/CustomNonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testSimpleExpression() {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomNonbondedForce* forceField = new CustomNonbondedForce("-0.1*r^3");
+    forceField->addParticle(vector<double>());
+    forceField->addParticle(vector<double>());
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double force = 0.1*3*(2*2);
+    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_TOL(-0.1*(2*2*2), state.getPotentialEnergy(), TOL);
+}
+
+void testParameters() {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomNonbondedForce* forceField = new CustomNonbondedForce("scale*a*(r*b)^3");
+    forceField->addParameter("a", "a1*a2");
+    forceField->addParameter("b", "c+b1+b2");
+    forceField->addGlobalParameter("scale", 3.0);
+    forceField->addGlobalParameter("c", -1.0);
+    vector<double> params(2);
+    params[0] = 1.5;
+    params[1] = 2.0;
+    forceField->addParticle(params);
+    params[0] = 2.0;
+    params[1] = 3.0;
+    forceField->addParticle(params);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+    context.setPositions(positions);
+    context.setParameter("scale", 1.0);
+    context.setParameter("c", 0.0);
+    State state = context.getState(State::Forces | State::Energy);
+    vector<Vec3> forces = state.getForces();
+    double force = -3.0*3*5.0*(10*10);
+    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_TOL(3.0*(10*10*10), state.getPotentialEnergy(), TOL);
+    context.setParameter("scale", 1.5);
+    context.setParameter("c", 1.0);
+    state = context.getState(State::Forces | State::Energy);
+    forces = state.getForces();
+    force = -1.5*3.0*3*6.0*(12*12);
+    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_TOL(1.5*3.0*(12*12*12), state.getPotentialEnergy(), TOL);
+}
+
+void testExceptions() {
+    OpenCLPlatform platform;
+    System system;
+    VerletIntegrator integrator(0.01);
+    CustomNonbondedForce* nonbonded = new CustomNonbondedForce("a*r");
+    nonbonded->addParameter("a", "a1+a2");
+    vector<double> params(1);
+    vector<Vec3> positions(4);
+    for (int i = 0; i < 4; i++) {
+        system.addParticle(1.0);
+        params[0] = i+1;
+        nonbonded->addParticle(params);
+        positions[i] = Vec3(i, 0, 0);
+    }
+    nonbonded->addException(0, 1, vector<double>());
+    nonbonded->addException(1, 2, vector<double>());
+    nonbonded->addException(2, 3, vector<double>());
+    params[0] = 0.5;
+    nonbonded->addException(0, 2, params);
+    nonbonded->addException(1, 3, params);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(0.5+1+4, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0.5, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-0.5, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-(0.5+1+4), 0, 0), forces[3], TOL);
+    ASSERT_EQUAL_TOL((1+4)*3+0.5*2+0.5*2, 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);
+    CustomNonbondedForce* forceField = new CustomNonbondedForce("r");
+    forceField->addParticle(vector<double>());
+    forceField->addParticle(vector<double>());
+    forceField->addParticle(vector<double>());
+    forceField->setNonbondedMethod(CustomNonbondedForce::CutoffNonPeriodic);
+    forceField->setCutoffDistance(2.5);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0, 2, 0);
+    positions[2] = Vec3(0, 3, 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), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -1, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(2.0+1.0, state.getPotentialEnergy(), TOL);
+}
+
+void testPeriodic() {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomNonbondedForce* forceField = new CustomNonbondedForce("r");
+    forceField->addParticle(vector<double>());
+    forceField->addParticle(vector<double>());
+    forceField->addParticle(vector<double>());
+    forceField->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
+    forceField->setCutoffDistance(2.0);
+    system.setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0, 2.1, 0);
+    positions[2] = Vec3(0, 3, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(0, -2, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 2, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(1.9+1+0.9, state.getPotentialEnergy(), TOL);
+}
+
+void testTabulatedFunction(bool interpolating) {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomNonbondedForce* forceField = new CustomNonbondedForce("fn(r)+1");
+    forceField->addParticle(vector<double>());
+    forceField->addParticle(vector<double>());
+    vector<double> table;
+    for (int i = 0; i < 21; i++)
+        table.push_back(std::sin(0.25*i));
+    forceField->addFunction("fn", table, 1.0, 6.0, interpolating);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    double tol = 0.01;
+    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);
+    }
+}
+
+int main() {
+    try {
+        testSimpleExpression();
+        testParameters();
+        testExceptions();
+        testCutoff();
+        testPeriodic();
+//        testTabulatedFunction(true);
+//        testTabulatedFunction(false);
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}