Created OpenCL implementation of CustomTorsionForce

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -50,6 +50,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
         return new OpenCLCalcPeriodicTorsionForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcRBTorsionForceKernel::Name())
         return new OpenCLCalcRBTorsionForceKernel(name, platform, cl, context.getSystem());
+    if (name == CalcCustomTorsionForceKernel::Name())
+        return new OpenCLCalcCustomTorsionForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcNonbondedForceKernel::Name())
         return new OpenCLCalcNonbondedForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcCustomNonbondedForceKernel::Name())

platforms/opencl/src/OpenCLKernels.cpp

@@ -571,7 +571,7 @@ void OpenCLCalcCustomAngleForceKernel::initialize(const System& system, const Cu
     Lepton::ParsedExpression forceExpression = energyExpression.differentiate("theta").optimize();
     map<string, Lepton::ParsedExpression> expressions;
     expressions["energy += "] = energyExpression;
-    expressions["float dEdR = "] = forceExpression;
+    expressions["float dEdAngle = "] = forceExpression;
 
     // Create the kernels.
 
@@ -807,6 +807,163 @@ double OpenCLCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
     return 0.0;
 }
 
+class OpenCLCustomTorsionForceInfo : public OpenCLForceInfo {
+public:
+    OpenCLCustomTorsionForceInfo(int requiredBuffers, const CustomTorsionForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
+    }
+    int getNumParticleGroups() {
+        return force.getNumTorsions();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2, particle3, particle4;
+        vector<double> parameters;
+        force.getTorsionParameters(index, particle1, particle2, particle3, particle4, parameters);
+        particles.resize(3);
+        particles[0] = particle1;
+        particles[1] = particle2;
+        particles[2] = particle3;
+        particles[3] = particle4;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2, particle3, particle4;
+        vector<double> parameters1, parameters2;
+        force.getTorsionParameters(group1, particle1, particle2, particle3, particle4, parameters1);
+        force.getTorsionParameters(group2, particle1, particle2, particle3, particle4, parameters2);
+        for (int i = 0; i < (int) parameters1.size(); i++)
+            if (parameters1[i] != parameters2[i])
+                return false;
+        return true;
+    }
+private:
+    const CustomTorsionForce& force;
+};
+
+OpenCLCalcCustomTorsionForceKernel::~OpenCLCalcCustomTorsionForceKernel() {
+    if (params != NULL)
+        delete params;
+    if (indices != NULL)
+        delete indices;
+    if (globals != NULL)
+        delete globals;
+}
+
+void OpenCLCalcCustomTorsionForceKernel::initialize(const System& system, const CustomTorsionForce& force) {
+    numTorsions = force.getNumTorsions();
+    if (numTorsions == 0)
+        return;
+    params = new OpenCLParameterSet(cl, force.getNumPerTorsionParameters(), numTorsions, "customTorsionParams");
+    indices = new OpenCLArray<mm_int8>(cl, numTorsions, "customTorsionIndices");
+    string extraArguments;
+    if (force.getNumGlobalParameters() > 0) {
+        globals = new OpenCLArray<cl_float>(cl, force.getNumGlobalParameters(), "customTorsionGlobals", false, CL_MEM_READ_ONLY);
+        extraArguments += ", __constant float* globals";
+    }
+    vector<int> forceBufferCounter(system.getNumParticles(), 0);
+    vector<vector<cl_float> > paramVector(numTorsions);
+    vector<mm_int8> indicesVector(numTorsions);
+    for (int i = 0; i < numTorsions; i++) {
+        int particle1, particle2, particle3, particle4;
+        vector<double> parameters;
+        force.getTorsionParameters(i, particle1, particle2, particle3, particle4, parameters);
+        paramVector[i].resize(parameters.size());
+        for (int j = 0; j < (int) parameters.size(); j++)
+            paramVector[i][j] = (cl_float) parameters[j];
+        indicesVector[i] = mm_int8(particle1, particle2, particle3, particle4, forceBufferCounter[particle1]++,
+                forceBufferCounter[particle2]++, forceBufferCounter[particle3]++, forceBufferCounter[particle4]++);
+    }
+    params->setParameterValues(paramVector);
+    indices->upload(indicesVector);
+    int maxBuffers = 1;
+    for (int i = 0; i < (int) forceBufferCounter.size(); i++)
+        maxBuffers = max(maxBuffers, forceBufferCounter[i]);
+    cl.addForce(new OpenCLCustomTorsionForceInfo(maxBuffers, force));
+
+    // Record information for the expressions.
+
+    vector<string> paramNames;
+    for (int i = 0; i < force.getNumPerTorsionParameters(); i++)
+        paramNames.push_back(force.getPerTorsionParameterName(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("theta").optimize();
+    map<string, Lepton::ParsedExpression> expressions;
+    expressions["energy += "] = energyExpression;
+    expressions["float dEdAngle = "] = forceExpression;
+
+    // Create the kernels.
+
+    map<string, string> variables;
+    variables["theta"] = "theta";
+    for (int i = 0; i < force.getNumPerTorsionParameters(); i++) {
+        const string& name = force.getPerTorsionParameterName(i);
+        variables[name] = "torsionParams"+params->getParameterSuffix(i);
+    }
+    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+        const string& name = force.getGlobalParameterName(i);
+        string value = "globals["+intToString(i)+"]";
+        variables[name] = value;
+    }
+    stringstream compute;
+    for (int i = 0; i < (int) params->getBuffers().size(); i++) {
+        const OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
+        extraArguments += ", __global "+buffer.getType()+"* "+buffer.getName();
+        compute<<buffer.getType()<<" torsionParams"<<(i+1)<<" = "<<buffer.getName()<<"[index];\n";
+    }
+    vector<pair<string, string> > functions;
+    compute << OpenCLExpressionUtilities::createExpressions(expressions, variables, functions, "temp", "");
+    map<string, string> replacements;
+    replacements["COMPUTE_FORCE"] = compute.str();
+    replacements["EXTRA_ARGUMENTS"] = extraArguments;
+    replacements["M_PI"] = doubleToString(M_PI);
+    cl::Program program = cl.createProgram(cl.replaceStrings(OpenCLKernelSources::customTorsionForce, replacements));
+    kernel = cl::Kernel(program, "computeCustomTorsionForces");
+}
+
+void OpenCLCalcCustomTorsionForceKernel::executeForces(ContextImpl& context) {
+    if (numTorsions == 0)
+        return;
+    if (globals != NULL) {
+        bool changed = false;
+        for (int i = 0; i < (int) 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, numTorsions);
+        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, indices->getDeviceBuffer());
+        int nextIndex = 6;
+        if (globals != NULL)
+            kernel.setArg<cl::Buffer>(nextIndex++, globals->getDeviceBuffer());
+        for (int i = 0; i < (int) params->getBuffers().size(); i++) {
+            const OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
+            kernel.setArg<cl::Buffer>(nextIndex++, buffer.getBuffer());
+        }
+    }
+    cl.executeKernel(kernel, numTorsions);
+}
+
+double OpenCLCalcCustomTorsionForceKernel::executeEnergy(ContextImpl& context) {
+    executeForces(context);
+    return 0.0;
+}
+
 class OpenCLNonbondedForceInfo : public OpenCLForceInfo {
 public:
     OpenCLNonbondedForceInfo(int requiredBuffers, const NonbondedForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {

platforms/opencl/src/OpenCLKernels.h

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

platforms/opencl/src/OpenCLPlatform.cpp

@@ -52,6 +52,7 @@ OpenCLPlatform::OpenCLPlatform() {
     registerKernelFactory(CalcCustomAngleForceKernel::Name(), factory);
     registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory);
     registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);

platforms/opencl/src/kernels/customAngleForce.cl

@@ -25,8 +25,8 @@ __kernel void computeCustomAngleForces(int numAtoms, int numAngles, __global flo
         float cosine = dot/sqrt(r21*r23);
         float theta = acos(cosine);
         COMPUTE_FORCE
-        float4 c21 = cross(v0, cp)*(dEdR/(r21*rp));
-        float4 c23 = cross(cp, v1)*(dEdR/(r23*rp));
+        float4 c21 = cross(v0, cp)*(dEdAngle/(r21*rp));
+        float4 c23 = cross(cp, v1)*(dEdAngle/(r23*rp));
 
         // Record the force on each of the three atoms.
 

platforms/opencl/src/kernels/customTorsionForce.cl

+/**
+ * Compute custom torsion forces.
+ */
+
+__kernel void computeCustomTorsionForces(int numAtoms, int numTorsions, __global float4* forceBuffers, __global float* energyBuffer,
+        __global float4* posq, __global int8* indices
+        EXTRA_ARGUMENTS) {
+    float energy = 0.0f;
+    for (int index = get_global_id(0); index < numTorsions; index += get_global_size(0)) {
+        int8 atoms = indices[index];
+        float4 a1 = posq[atoms.s0];
+        float4 a2 = posq[atoms.s1];
+        float4 a3 = posq[atoms.s2];
+        float4 a4 = posq[atoms.s3];
+
+        // Compute the force.
+
+        float4 v0 = (float4) (a1.xyz-a2.xyz, 0.0f);
+        float4 v1 = (float4) (a3.xyz-a2.xyz, 0.0f);
+        float4 v2 = (float4) (a3.xyz-a4.xyz, 0.0f);
+        float4 cp0 = cross(v0, v1);
+        float4 cp1 = cross(v1, v2);
+        float cosangle = dot(normalize(cp0), normalize(cp1));
+        float theta;
+        if (cosangle > 0.99f || cosangle < -0.99f) {
+            // We're close to the singularity in acos(), so take the cross product and use asin() instead.
+
+            float4 cross_prod = cross(cp0, cp1);
+            float scale = dot(cp0, cp0)*dot(cp1, cp1);
+            theta = asin(sqrt(dot(cross_prod, cross_prod)/scale));
+            if (cosangle < 0.0f)
+                theta = M_PI-theta;
+        }
+        else
+           theta = acos(cosangle);
+        theta = (dot(v0, cp1) >= 0 ? theta : -theta);
+        COMPUTE_FORCE
+        float normCross1 = dot(cp0, cp0);
+        float normSqrBC = dot(v1, v1);
+        float normBC = sqrt(normSqrBC);
+        float normCross2 = dot(cp1, cp1);
+        float dp = 1.0f/normSqrBC;
+        float4 ff = (float4) ((-dEdAngle*normBC)/normCross1, dot(v0, v1)*dp, dot(v2, v1)*dp, (dEdAngle*normBC)/normCross2);
+        float4 internalF0 = ff.x*cp0;
+        float4 internalF3 = ff.w*cp1;
+        float4 s = ff.y*internalF0 - ff.z*internalF3;
+
+        // Record the force on each of the four atoms.
+
+        unsigned int offsetA = atoms.s0+atoms.s4*numAtoms;
+        unsigned int offsetB = atoms.s1+atoms.s5*numAtoms;
+        unsigned int offsetC = atoms.s2+atoms.s6*numAtoms;
+        unsigned int offsetD = atoms.s3+atoms.s7*numAtoms;
+        float4 forceA = forceBuffers[offsetA];
+        float4 forceB = forceBuffers[offsetB];
+        float4 forceC = forceBuffers[offsetC];
+        float4 forceD = forceBuffers[offsetD];
+        forceA.xyz += internalF0.xyz;
+        forceB.xyz += s.xyz-internalF0.xyz;
+        forceC.xyz += -s.xyz-internalF3.xyz;
+        forceD.xyz += internalF3.xyz;
+        forceBuffers[offsetA] = forceA;
+        forceBuffers[offsetB] = forceB;
+        forceBuffers[offsetC] = forceC;
+        forceBuffers[offsetD] = forceD;
+    }
+    energyBuffer[get_global_id(0)] += energy;
+}

platforms/opencl/tests/TestOpenCLCustomTorsionForce.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 CustomTorsionForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "OpenCLPlatform.h"
+#include "openmm/CustomTorsionForce.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 testTorsions() {
+    OpenCLPlatform platform;
+
+    // Create a system using a CustomTorsionForce.
+
+    System customSystem;
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    CustomTorsionForce* custom = new CustomTorsionForce("k*(1+cos(n*theta-theta0))");
+    custom->addPerTorsionParameter("theta0");
+    custom->addPerTorsionParameter("n");
+    custom->addGlobalParameter("k", 0.5);
+    vector<double> parameters(2);
+    parameters[0] = 1.5;
+    parameters[1] = 1;
+    custom->addTorsion(0, 1, 2, 3, parameters);
+    parameters[0] = 2.0;
+    parameters[1] = 2;
+    custom->addTorsion(1, 2, 3, 4, parameters);
+    customSystem.addForce(custom);
+
+    // Create an identical system using a PeriodicTorsionForce.
+
+    System harmonicSystem;
+    harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    PeriodicTorsionForce* periodic = new PeriodicTorsionForce();
+    periodic->addTorsion(0, 1, 2, 3, 1, 1.5, 0.5);
+    periodic->addTorsion(1, 2, 3, 4, 2, 2.0, 0.5);
+    harmonicSystem.addForce(periodic);
+
+    // 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(harmonicSystem, integrator2, platform);
+        for (int j = 0; j < (int) positions.size(); j++)
+            positions[j] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.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(s1.getForces()[i], s2.getForces()[i], TOL);
+        ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
+    }
+}
+
+int main() {
+    try {
+        testTorsions();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+