Created OpenCL implementation of CustomBondForce

platforms/cuda/tests/TestCudaCustomNonbondedForce.cpp

@@ -7,7 +7,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,9 +35,11 @@
  */
 
 #include "../../../tests/AssertionUtilities.h"
+#include "../src/sfmt/SFMT.h"
 #include "openmm/Context.h"
 #include "CudaPlatform.h"
 #include "openmm/CustomNonbondedForce.h"
+#include "openmm/NonbondedForce.h"
 #include "openmm/System.h"
 #include "openmm/VerletIntegrator.h"
 #include <iostream>
@@ -234,6 +236,82 @@ void testTabulatedFunction(bool interpolating) {
     }
 }
 
+void testCoulombLennardJones() {
+    const int numMolecules = 300;
+    const int numParticles = numMolecules*2;
+    const double boxSize = 20.0;
+    CudaPlatform platform;
+
+    // Create two systems: one with a NonbondedForce, and one using a CustomNonbondedForce to implement the same interaction.
+
+    System standardSystem;
+    System customSystem;
+    for (int i = 0; i < numParticles; i++) {
+        standardSystem.addParticle(1.0);
+        customSystem.addParticle(1.0);
+    }
+    NonbondedForce* standardNonbonded = new NonbondedForce();
+    CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6) + 138.935485*q/r");
+    customNonbonded->addParameter("q", "q1*q2");
+    customNonbonded->addParameter("sigma", "0.5*(sigma1+sigma2)");
+    customNonbonded->addParameter("eps", "sqrt(eps1*eps2)");
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    init_gen_rand(0);
+    vector<double> params(3);
+    for (int i = 0; i < numMolecules; i++) {
+        if (i < numMolecules/2) {
+            standardNonbonded->addParticle(1.0, 0.2, 0.1);
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.1;
+            customNonbonded->addParticle(params);
+            standardNonbonded->addParticle(1.0, 0.1, 0.1);
+            params[1] = 0.1;
+            customNonbonded->addParticle(params);
+        }
+        else {
+            standardNonbonded->addParticle(1.0, 0.2, 0.2);
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.2;
+            customNonbonded->addParticle(params);
+            standardNonbonded->addParticle(1.0, 0.1, 0.2);
+            params[1] = 0.1;
+            customNonbonded->addParticle(params);
+        }
+        positions[2*i] = Vec3(boxSize*genrand_real2(), boxSize*genrand_real2(), boxSize*genrand_real2());
+        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
+        velocities[2*i] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
+        velocities[2*i+1] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
+        standardNonbonded->addException(2*i, 2*i+1, 0.0, 1.0, 0.0);
+        customNonbonded->addException(2*i, 2*i+1, vector<double>());
+    }
+    standardNonbonded->setNonbondedMethod(NonbondedForce::NoCutoff);
+    customNonbonded->setNonbondedMethod(CustomNonbondedForce::NoCutoff);
+    standardSystem.addForce(standardNonbonded);
+    customSystem.addForce(customNonbonded);
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    Context context1(standardSystem, integrator1, platform);
+//    Context context2(customSystem, integrator2, platform);
+    context1.setPositions(positions);
+//    context2.setPositions(positions);
+    context1.setVelocities(velocities);
+//    context2.setVelocities(velocities);
+    State state1 = context1.getState(State::Forces | State::Energy);
+//    State state2 = context2.getState(State::Forces | State::Energy);
+    for (int i = 0; i < numParticles; i++)
+        std::cout << i<<": "<<state1.getForces()[i]<< std::endl;
+//    for (int i = 0; i < numParticles; i++)
+//        std::cout << i<<": "<<state1.getForces()[i]<<" "<<state2.getForces()[i]<< std::endl;
+//    std::cout <<state1.getPotentialEnergy()<<" "<<state2.getPotentialEnergy() << std::endl;
+//    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
+//    for (int i = 0; i < numParticles; i++) {
+//        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
+//    }
+}
+
 int main() {
     try {
         testSimpleExpression();
@@ -243,6 +321,7 @@ int main() {
         testPeriodic();
         testTabulatedFunction(true);
         testTabulatedFunction(false);
+//        testCoulombLennardJones();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -40,6 +40,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
         return new OpenCLUpdateStateDataKernel(name, platform, cl);
     if (name == CalcHarmonicBondForceKernel::Name())
         return new OpenCLCalcHarmonicBondForceKernel(name, platform, cl, context.getSystem());
+    if (name == CalcCustomBondForceKernel::Name())
+        return new OpenCLCalcCustomBondForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcHarmonicAngleForceKernel::Name())
         return new OpenCLCalcHarmonicAngleForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcPeriodicTorsionForceKernel::Name())

platforms/opencl/src/OpenCLKernels.cpp

@@ -249,6 +249,154 @@ double OpenCLCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
     return 0.0;
 }
 
+class OpenCLCustomBondForceInfo : public OpenCLForceInfo {
+public:
+    OpenCLCustomBondForceInfo(int requiredBuffers, const CustomBondForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
+    }
+    int getNumParticleGroups() {
+        return force.getNumBonds();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2;
+        vector<double> parameters;
+        force.getBondParameters(index, particle1, particle2, parameters);
+        particles.resize(2);
+        particles[0] = particle1;
+        particles[1] = particle2;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2;
+        vector<double> parameters1, parameters2;
+        force.getBondParameters(group1, particle1, particle2, parameters1);
+        force.getBondParameters(group2, particle1, particle2, parameters2);
+        for (int i = 0; i < (int) parameters1.size(); i++)
+            if (parameters1[i] != parameters2[i])
+                return false;
+        return true;
+    }
+private:
+    const CustomBondForce& force;
+};
+
+OpenCLCalcCustomBondForceKernel::~OpenCLCalcCustomBondForceKernel() {
+    if (params != NULL)
+        delete params;
+    if (indices != NULL)
+        delete indices;
+    if (globals != NULL)
+        delete globals;
+}
+
+void OpenCLCalcCustomBondForceKernel::initialize(const System& system, const CustomBondForce& force) {
+    if (force.getNumPerBondParameters() > 4)
+        throw OpenMMException("OpenCLPlatform only supports four per-bond parameters for custom bonded forces");
+    numBonds = force.getNumBonds();
+    params = new OpenCLArray<mm_float4>(cl, numBonds, "customBondParams");
+    indices = new OpenCLArray<mm_int4>(cl, numBonds, "customBondIndices");
+    string extraArguments;
+    if (force.getNumGlobalParameters() > 0) {
+        globals = new OpenCLArray<cl_float>(cl, force.getNumGlobalParameters(), "customBondGlobals", false, CL_MEM_READ_ONLY);
+        extraArguments += ", __constant float* globals";
+    }
+    vector<int> forceBufferCounter(system.getNumParticles(), 0);
+    vector<mm_float4> paramVector(numBonds);
+    vector<mm_int4> indicesVector(numBonds);
+    for (int i = 0; i < numBonds; i++) {
+        int particle1, particle2;
+        vector<double> parameters;
+        force.getBondParameters(i, particle1, particle2, parameters);
+        if (parameters.size() > 0)
+            paramVector[i].x = (cl_float) parameters[0];
+        if (parameters.size() > 1)
+            paramVector[i].y = (cl_float) parameters[1];
+        if (parameters.size() > 2)
+            paramVector[i].z = (cl_float) parameters[2];
+        if (parameters.size() > 3)
+            paramVector[i].w = (cl_float) parameters[3];
+        indicesVector[i] = (mm_int4) {particle1, particle2, forceBufferCounter[particle1]++, forceBufferCounter[particle2]++};
+    }
+    params->upload(paramVector);
+    indices->upload(indicesVector);
+    int maxBuffers = 1;
+    for (int i = 0; i < forceBufferCounter.size(); i++)
+        maxBuffers = max(maxBuffers, forceBufferCounter[i]);
+    cl.addForce(new OpenCLCustomBondForceInfo(maxBuffers, force));
+
+    // Record information for the expressions.
+
+    vector<string> paramNames;
+    for (int i = 0; i < force.getNumPerBondParameters(); i++)
+        paramNames.push_back(force.getPerBondParameterName(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);
+    Lepton::ParsedExpression energyExpression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
+    Lepton::ParsedExpression forceExpression = energyExpression.differentiate("r").optimize();
+    map<string, Lepton::ParsedExpression> expressions;
+    expressions["energy += "] = energyExpression;
+    expressions["float dEdR = "] = forceExpression;
+
+    // Create the kernels.
+
+    map<string, string> variables;
+    variables["r"] = "r";
+    string suffixes[] = {".x", ".y", ".z", ".w"};
+    for (int i = 0; i < force.getNumPerBondParameters(); i++) {
+        const string& name = force.getPerBondParameterName(i);
+        variables[name] = "exceptionParams"+suffixes[i];
+    }
+    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+        const string& name = force.getGlobalParameterName(i);
+        string value = "globals["+intToString(i)+"]";
+        variables[name] = value;
+    }
+    map<string, string> replacements;
+    stringstream compute;
+    vector<pair<string, string> > functions;
+    compute << OpenCLExpressionUtilities::createExpressions(expressions, variables, functions, "temp", "");
+    replacements["COMPUTE_FORCE"] = compute.str();
+    replacements["EXTRA_ARGUMENTS"] = extraArguments;
+    cl::Program program = cl.createProgram(cl.loadSourceFromFile("customBondForce.cl", replacements));
+    kernel = cl::Kernel(program, "computeCustomBondForces");
+}
+
+void OpenCLCalcCustomBondForceKernel::executeForces(ContextImpl& context) {
+    if (globals != NULL) {
+        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);
+    }
+    if (!hasInitializedKernel) {
+        hasInitializedKernel = true;
+        kernel.setArg<cl_int>(0, cl.getPaddedNumAtoms());
+        kernel.setArg<cl_int>(1, numBonds);
+        kernel.setArg<cl::Buffer>(2, cl.getForceBuffers().getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(3, cl.getEnergyBuffer().getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(4, cl.getPosq().getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(5, params->getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(6, indices->getDeviceBuffer());
+        if (globals != NULL)
+            kernel.setArg<cl::Buffer>(7, globals->getDeviceBuffer());
+    }
+    cl.executeKernel(kernel, numBonds);
+}
+
+double OpenCLCalcCustomBondForceKernel::executeEnergy(ContextImpl& context) {
+    executeForces(context);
+    return 0.0;
+}
+
 class OpenCLAngleForceInfo : public OpenCLForceInfo {
 public:
     OpenCLAngleForceInfo(int requiredBuffers, const HarmonicAngleForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
@@ -987,6 +1135,17 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
 }
 
 void OpenCLCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
+    if (globals != NULL) {
+        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);
+    }
     if (exceptionParams != NULL) {
         if (!hasInitializedKernel) {
             hasInitializedKernel = true;
@@ -1002,17 +1161,6 @@ void OpenCLCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
         }
         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) {

platforms/opencl/src/OpenCLKernels.h

@@ -184,6 +184,48 @@ private:
     cl::Kernel kernel;
 };
 
+/**
+ * This kernel is invoked by CustomBondForce to calculate the forces acting on the system and the energy of the system.
+ */
+class OpenCLCalcCustomBondForceKernel : public CalcCustomBondForceKernel {
+public:
+    OpenCLCalcCustomBondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomBondForceKernel(name, platform),
+            hasInitializedKernel(false), cl(cl), system(system), params(NULL), indices(NULL), globals(NULL) {
+    }
+    ~OpenCLCalcCustomBondForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomBondForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomBondForce& 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 CustomBondForce
+     */
+    double executeEnergy(ContextImpl& context);
+private:
+    int numBonds;
+    bool hasInitializedKernel;
+    OpenCLContext& cl;
+    System& system;
+    OpenCLArray<mm_float4>* params;
+    OpenCLArray<mm_int4>* indices;
+    OpenCLArray<cl_float>* globals;
+    std::vector<std::string> globalParamNames;
+    std::vector<cl_float> globalParamValues;
+    cl::Kernel kernel;
+};
+
 /**
  * This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
  */

platforms/opencl/src/OpenCLPlatform.cpp

@@ -48,6 +48,7 @@ OpenCLPlatform::OpenCLPlatform() {
     registerKernelFactory(CalcForcesAndEnergyKernel::Name(), factory);
     registerKernelFactory(UpdateStateDataKernel::Name(), factory);
     registerKernelFactory(CalcHarmonicBondForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomBondForceKernel::Name(), factory);
     registerKernelFactory(CalcHarmonicAngleForceKernel::Name(), factory);
     registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);

platforms/opencl/src/kernels/customBondForce.cl

+/**
+ * Compute custom bond forces.
+ */
+
+__kernel void computeCustomBondForces(int numAtoms, int numBonds, __global float4* forceBuffers, __global float* energyBuffer,
+        __global float4* posq, __global float4* params, __global int4* indices
+        EXTRA_ARGUMENTS) {
+    int index = get_global_id(0);
+    float energy = 0.0f;
+    while (index < numBonds) {
+        // Look up the data for this exception.
+
+        int4 atoms = indices[index];
+        float4 exceptionParams = params[index];
+        float4 delta = posq[atoms.y]-posq[atoms.x];
+
+        // Compute the force.
+
+        float r = sqrt(delta.x*delta.x + delta.y*delta.y + delta.z*delta.z);
+        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/TestOpenCLCustomBondForce.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 the OpenCL implementation of CustomBondForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "OpenCLPlatform.h"
+#include "openmm/CustomBondForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testBonds() {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomBondForce* forceField = new CustomBondForce("scale*k*(r-r0)^2");
+    forceField->addPerBondParameter("r0");
+    forceField->addPerBondParameter("k");
+    forceField->addGlobalParameter("scale", 0.5);
+    vector<double> parameters(2);
+    parameters[0] = 1.5;
+    parameters[1] = 0.8;
+    forceField->addBond(0, 1, parameters);
+    parameters[0] = 1.2;
+    parameters[1] = 0.7;
+    forceField->addBond(1, 2, parameters);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 2, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(0, -0.8*0.5, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0.7*0.2, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-forces[0][0]-forces[2][0], -forces[0][1]-forces[2][1], -forces[0][2]-forces[2][2]), forces[1], TOL);
+    ASSERT_EQUAL_TOL(0.5*0.8*0.5*0.5 + 0.5*0.7*0.2*0.2, state.getPotentialEnergy(), TOL);
+}
+
+int main() {
+    try {
+        testBonds();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+
+

platforms/reference/tests/TestReferenceCustomBondForce.cpp

@@ -30,7 +30,7 @@
  * -------------------------------------------------------------------------- */
 
 /**
- * This tests all the different force terms in the reference implementation of HarmonicBondForce.
+ * This tests the reference implementation of CustomBondForce.
  */
 
 #include "../../../tests/AssertionUtilities.h"