Converted more forces to use OpenCLBondedUtilities

platforms/opencl/src/OpenCLKernels.cpp

@@ -290,8 +290,6 @@ private:
 OpenCLCalcCustomBondForceKernel::~OpenCLCalcCustomBondForceKernel() {
     if (params != NULL)
         delete params;
-    if (indices != NULL)
-        delete indices;
     if (globals != NULL)
         delete globals;
 }
@@ -303,31 +301,18 @@ void OpenCLCalcCustomBondForceKernel::initialize(const System& system, const Cus
     numBonds = endIndex-startIndex;
     if (numBonds == 0)
         return;
+    vector<vector<int> > atoms(numBonds, vector<int>(2));
     params = new OpenCLParameterSet(cl, force.getNumPerBondParameters(), 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<vector<cl_float> > paramVector(numBonds);
-    vector<mm_int4> indicesVector(numBonds);
     for (int i = 0; i < numBonds; i++) {
-        int particle1, particle2;
         vector<double> parameters;
-        force.getBondParameters(startIndex+i, particle1, particle2, parameters);
+        force.getBondParameters(startIndex+i, atoms[i][0], atoms[i][1], 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_int4(particle1, particle2, forceBufferCounter[particle1]++, forceBufferCounter[particle2]++);
     }
     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 OpenCLCustomBondForceInfo(maxBuffers, force));
+    cl.addForce(new OpenCLCustomBondForceInfo(0, force));
 
     // Record information for the expressions.
 
@@ -337,8 +322,6 @@ void OpenCLCalcCustomBondForceKernel::initialize(const System& system, const Cus
         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;
@@ -353,29 +336,30 @@ void OpenCLCalcCustomBondForceKernel::initialize(const System& system, const Cus
         const string& name = force.getPerBondParameterName(i);
         variables[name] = "bondParams"+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;
+    if (force.getNumGlobalParameters() > 0) {
+        globals = new OpenCLArray<cl_float>(cl, force.getNumGlobalParameters(), "customBondGlobals", false, CL_MEM_READ_ONLY);
+        globals->upload(globalParamValues);
+        string argName = cl.getBondedUtilities().addArgument(globals->getDeviceBuffer(), "float");
+        for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+            const string& name = force.getGlobalParameterName(i);
+            string value = argName+"["+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()<<" bondParams"<<(i+1)<<" = "<<buffer.getName()<<"[index];\n";
+        string argName = cl.getBondedUtilities().addArgument(buffer.getMemory(), buffer.getType());
+        compute<<buffer.getType()<<" bondParams"<<(i+1)<<" = "<<argName<<"[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;
-    cl::Program program = cl.createProgram(cl.replaceStrings(OpenCLKernelSources::customBondForce, replacements));
-    kernel = cl::Kernel(program, "computeCustomBondForces");
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::customBondForce, replacements));
 }
 
 double OpenCLCalcCustomBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    if (numBonds == 0)
-        return 0.0;
     if (globals != NULL) {
         bool changed = false;
         for (int i = 0; i < (int) globalParamNames.size(); i++) {
@@ -387,23 +371,6 @@ double OpenCLCalcCustomBondForceKernel::execute(ContextImpl& context, bool inclu
         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, 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::Memory>(nextIndex++, buffer.getMemory());
-        }
-    }
-    cl.executeKernel(kernel, numBonds);
     return 0.0;
 }
 
@@ -499,8 +466,6 @@ private:
 OpenCLCalcCustomAngleForceKernel::~OpenCLCalcCustomAngleForceKernel() {
     if (params != NULL)
         delete params;
-    if (indices != NULL)
-        delete indices;
     if (globals != NULL)
         delete globals;
 }
@@ -512,32 +477,18 @@ void OpenCLCalcCustomAngleForceKernel::initialize(const System& system, const Cu
     numAngles = endIndex-startIndex;
     if (numAngles == 0)
         return;
+    vector<vector<int> > atoms(numAngles, vector<int>(3));
     params = new OpenCLParameterSet(cl, force.getNumPerAngleParameters(), numAngles, "customAngleParams");
-    indices = new OpenCLArray<mm_int8>(cl, numAngles, "customAngleIndices");
-    string extraArguments;
-    if (force.getNumGlobalParameters() > 0) {
-        globals = new OpenCLArray<cl_float>(cl, force.getNumGlobalParameters(), "customAngleGlobals", false, CL_MEM_READ_ONLY);
-        extraArguments += ", __constant float* globals";
-    }
-    vector<int> forceBufferCounter(system.getNumParticles(), 0);
     vector<vector<cl_float> > paramVector(numAngles);
-    vector<mm_int8> indicesVector(numAngles);
     for (int i = 0; i < numAngles; i++) {
-        int particle1, particle2, particle3;
         vector<double> parameters;
-        force.getAngleParameters(startIndex+i, particle1, particle2, particle3, parameters);
+        force.getAngleParameters(startIndex+i, atoms[i][0], atoms[i][1], atoms[i][2], 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, forceBufferCounter[particle1]++,
-                forceBufferCounter[particle2]++, forceBufferCounter[particle3]++, 0, 0);
     }
     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 OpenCLCustomAngleForceInfo(maxBuffers, force));
+    cl.addForce(new OpenCLCustomAngleForceInfo(0, force));
 
     // Record information for the expressions.
 
@@ -547,8 +498,6 @@ void OpenCLCalcCustomAngleForceKernel::initialize(const System& system, const Cu
         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;
@@ -563,29 +512,30 @@ void OpenCLCalcCustomAngleForceKernel::initialize(const System& system, const Cu
         const string& name = force.getPerAngleParameterName(i);
         variables[name] = "angleParams"+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;
+    if (force.getNumGlobalParameters() > 0) {
+        globals = new OpenCLArray<cl_float>(cl, force.getNumGlobalParameters(), "customAngleGlobals", false, CL_MEM_READ_ONLY);
+        globals->upload(globalParamValues);
+        string argName = cl.getBondedUtilities().addArgument(globals->getDeviceBuffer(), "float");
+        for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+            const string& name = force.getGlobalParameterName(i);
+            string value = argName+"["+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()<<" angleParams"<<(i+1)<<" = "<<buffer.getName()<<"[index];\n";
+        string argName = cl.getBondedUtilities().addArgument(buffer.getMemory(), buffer.getType());
+        compute<<buffer.getType()<<" angleParams"<<(i+1)<<" = "<<argName<<"[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;
-    cl::Program program = cl.createProgram(cl.replaceStrings(OpenCLKernelSources::customAngleForce, replacements));
-    kernel = cl::Kernel(program, "computeCustomAngleForces");
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::customAngleForce, replacements));
 }
 
 double OpenCLCalcCustomAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    if (numAngles == 0)
-        return 0.0;
     if (globals != NULL) {
         bool changed = false;
         for (int i = 0; i < (int) globalParamNames.size(); i++) {
@@ -597,23 +547,6 @@ double OpenCLCalcCustomAngleForceKernel::execute(ContextImpl& context, bool incl
         if (changed)
             globals->upload(globalParamValues);
     }
-    if (!hasInitializedKernel) {
-        hasInitializedKernel = true;
-        kernel.setArg<cl_int>(0, cl.getPaddedNumAtoms());
-        kernel.setArg<cl_int>(1, numAngles);
-        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::Memory>(nextIndex++, buffer.getMemory());
-        }
-    }
-    cl.executeKernel(kernel, numAngles);
     return 0.0;
 }
 
@@ -775,8 +708,6 @@ OpenCLCalcCMAPTorsionForceKernel::~OpenCLCalcCMAPTorsionForceKernel() {
         delete mapPositions;
     if (torsionMaps != NULL)
         delete torsionMaps;
-    if (torsionIndices != NULL)
-        delete torsionIndices;
 }
 
 void OpenCLCalcCMAPTorsionForceKernel::initialize(const System& system, const CMAPTorsionForce& force) {
@@ -805,53 +736,25 @@ void OpenCLCalcCMAPTorsionForceKernel::initialize(const System& system, const CM
             coeffVec.push_back(mm_float4((float) c[j][12], (float) c[j][13], (float) c[j][14], (float) c[j][15]));
         }
     }
-    vector<int> forceBufferCounter(system.getNumParticles(), 0);
+    vector<vector<int> > atoms(numTorsions, vector<int>(8));
     vector<cl_int> torsionMapsVec(numTorsions);
-    vector<mm_int16> torsionIndicesVec(numTorsions);
-    for (int i = 0; i < numTorsions; i++) {
-        mm_int16& ind = torsionIndicesVec[i];
-        force.getTorsionParameters(startIndex+i, torsionMapsVec[i], ind.s0, ind.s1, ind.s2, ind.s3, ind.s4, ind.s5, ind.s6, ind.s7);
-        ind.s8 = forceBufferCounter[ind.s0]++;
-        ind.s9 = forceBufferCounter[ind.s1]++;
-        ind.s10 = forceBufferCounter[ind.s2]++;
-        ind.s11 = forceBufferCounter[ind.s3]++;
-        ind.s12 = forceBufferCounter[ind.s4]++;
-        ind.s13 = forceBufferCounter[ind.s5]++;
-        ind.s14 = forceBufferCounter[ind.s6]++;
-        ind.s15 = forceBufferCounter[ind.s7]++;
-    }
+    for (int i = 0; i < numTorsions; i++)
+        force.getTorsionParameters(startIndex+i, torsionMapsVec[i], atoms[i][0], atoms[i][1], atoms[i][2], atoms[i][3], atoms[i][4], atoms[i][5], atoms[i][6], atoms[i][7]);
     coefficients = new OpenCLArray<mm_float4>(cl, coeffVec.size(), "cmapTorsionCoefficients");
     mapPositions = new OpenCLArray<mm_int2>(cl, numMaps, "cmapTorsionMapPositions");
     torsionMaps = new OpenCLArray<cl_int>(cl, numTorsions, "cmapTorsionMaps");
-    torsionIndices = new OpenCLArray<mm_int16>(cl, numTorsions, "cmapTorsionIndices");
     coefficients->upload(coeffVec);
     mapPositions->upload(mapPositionsVec);
     torsionMaps->upload(torsionMapsVec);
-    torsionIndices->upload(torsionIndicesVec);
-    int maxBuffers = 1;
-    for (int i = 0; i < (int) forceBufferCounter.size(); i++)
-        maxBuffers = max(maxBuffers, forceBufferCounter[i]);
-    cl.addForce(new OpenCLCMAPTorsionForceInfo(maxBuffers, force));
-    cl::Program program = cl.createProgram(OpenCLKernelSources::cmapTorsionForce);
-    kernel = cl::Kernel(program, "computeCMAPTorsionForces");
+    map<string, string> replacements;
+    replacements["COEFF"] = cl.getBondedUtilities().addArgument(coefficients->getDeviceBuffer(), "float4");
+    replacements["MAP_POS"] = cl.getBondedUtilities().addArgument(mapPositions->getDeviceBuffer(), "int2");
+    replacements["MAPS"] = cl.getBondedUtilities().addArgument(torsionMaps->getDeviceBuffer(), "int");
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::cmapTorsionForce, replacements));
+    cl.addForce(new OpenCLCMAPTorsionForceInfo(0, force));
 }
 
 double OpenCLCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    if (numTorsions == 0)
-        return 0.0;
-    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, coefficients->getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(6, mapPositions->getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(7, torsionIndices->getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(8, torsionMaps->getDeviceBuffer());
-    }
-    cl.executeKernel(kernel, numTorsions);
     return 0.0;
 }
 
@@ -889,8 +792,6 @@ private:
 OpenCLCalcCustomTorsionForceKernel::~OpenCLCalcCustomTorsionForceKernel() {
     if (params != NULL)
         delete params;
-    if (indices != NULL)
-        delete indices;
     if (globals != NULL)
         delete globals;
 }
@@ -902,32 +803,18 @@ void OpenCLCalcCustomTorsionForceKernel::initialize(const System& system, const
     numTorsions = endIndex-startIndex;
     if (numTorsions == 0)
         return;
+    vector<vector<int> > atoms(numTorsions, vector<int>(4));
     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(startIndex+i, particle1, particle2, particle3, particle4, parameters);
+        force.getTorsionParameters(startIndex+i, atoms[i][0], atoms[i][1], atoms[i][2], atoms[i][3], 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));
+    cl.addForce(new OpenCLCustomTorsionForceInfo(0, force));
 
     // Record information for the expressions.
 
@@ -937,8 +824,6 @@ void OpenCLCalcCustomTorsionForceKernel::initialize(const System& system, const
         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;
@@ -953,30 +838,31 @@ void OpenCLCalcCustomTorsionForceKernel::initialize(const System& system, const
         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;
+    if (force.getNumGlobalParameters() > 0) {
+        globals = new OpenCLArray<cl_float>(cl, force.getNumGlobalParameters(), "customTorsionGlobals", false, CL_MEM_READ_ONLY);
+        globals->upload(globalParamValues);
+        string argName = cl.getBondedUtilities().addArgument(globals->getDeviceBuffer(), "float");
+        for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+            const string& name = force.getGlobalParameterName(i);
+            string value = argName+"["+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";
+        string argName = cl.getBondedUtilities().addArgument(buffer.getMemory(), buffer.getType());
+        compute<<buffer.getType()<<" torsionParams"<<(i+1)<<" = "<<argName<<"[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");
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::customTorsionForce, replacements));
 }
 
 double OpenCLCalcCustomTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    if (numTorsions == 0)
-        return 0.0;
     if (globals != NULL) {
         bool changed = false;
         for (int i = 0; i < (int) globalParamNames.size(); i++) {
@@ -988,23 +874,6 @@ double OpenCLCalcCustomTorsionForceKernel::execute(ContextImpl& context, bool in
         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::Memory>(nextIndex++, buffer.getMemory());
-        }
-    }
-    cl.executeKernel(kernel, numTorsions);
     return 0.0;
 }
 
@@ -2632,8 +2501,6 @@ private:
 OpenCLCalcCustomExternalForceKernel::~OpenCLCalcCustomExternalForceKernel() {
     if (params != NULL)
         delete params;
-    if (indices != NULL)
-        delete indices;
     if (globals != NULL)
         delete globals;
 }
@@ -2645,24 +2512,17 @@ void OpenCLCalcCustomExternalForceKernel::initialize(const System& system, const
     numParticles = endIndex-startIndex;
     if (numParticles == 0)
         return;
+    vector<vector<int> > atoms(numParticles, vector<int>(1));
     params = new OpenCLParameterSet(cl, force.getNumPerParticleParameters(), numParticles, "customExternalParams");
-    indices = new OpenCLArray<cl_int>(cl, numParticles, "customExternalIndices");
-    string extraArguments;
-    if (force.getNumGlobalParameters() > 0) {
-        globals = new OpenCLArray<cl_float>(cl, force.getNumGlobalParameters(), "customExternalGlobals", false, CL_MEM_READ_ONLY);
-        extraArguments += ", __constant float* globals";
-    }
     vector<vector<cl_float> > paramVector(numParticles);
-    vector<cl_int> indicesVector(numParticles);
     for (int i = 0; i < numParticles; i++) {
         vector<double> parameters;
-        force.getParticleParameters(startIndex+i, indicesVector[i], parameters);
+        force.getParticleParameters(startIndex+i, atoms[i][0], parameters);
         paramVector[i].resize(parameters.size());
         for (int j = 0; j < (int) parameters.size(); j++)
             paramVector[i][j] = (cl_float) parameters[j];
     }
     params->setParameterValues(paramVector);
-    indices->upload(indicesVector);
     cl.addForce(new OpenCLCustomExternalForceInfo(force, system.getNumParticles()));
 
     // Record information for the expressions.
@@ -2673,8 +2533,6 @@ void OpenCLCalcCustomExternalForceKernel::initialize(const System& system, const
         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 forceExpressionX = energyExpression.differentiate("x").optimize();
     Lepton::ParsedExpression forceExpressionY = energyExpression.differentiate("y").optimize();
@@ -2688,36 +2546,37 @@ void OpenCLCalcCustomExternalForceKernel::initialize(const System& system, const
     // Create the kernels.
 
     map<string, string> variables;
-    variables["x"] = "pos.x";
-    variables["y"] = "pos.y";
-    variables["z"] = "pos.z";
+    variables["x"] = "pos1.x";
+    variables["y"] = "pos1.y";
+    variables["z"] = "pos1.z";
     for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
         const string& name = force.getPerParticleParameterName(i);
         variables[name] = "particleParams"+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;
+    if (force.getNumGlobalParameters() > 0) {
+        globals = new OpenCLArray<cl_float>(cl, force.getNumGlobalParameters(), "customExternalGlobals", false, CL_MEM_READ_ONLY);
+        globals->upload(globalParamValues);
+        string argName = cl.getBondedUtilities().addArgument(globals->getDeviceBuffer(), "float");
+        for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+            const string& name = force.getGlobalParameterName(i);
+            string value = argName+"["+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()<<" particleParams"<<(i+1)<<" = "<<buffer.getName()<<"[index];\n";
+        string argName = cl.getBondedUtilities().addArgument(buffer.getMemory(), buffer.getType());
+        compute<<buffer.getType()<<" particleParams"<<(i+1)<<" = "<<argName<<"[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;
-    cl::Program program = cl.createProgram(cl.replaceStrings(OpenCLKernelSources::customExternalForce, replacements));
-    kernel = cl::Kernel(program, "computeCustomExternalForces");
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::customExternalForce, replacements));
 }
 
 double OpenCLCalcCustomExternalForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    if (numParticles == 0)
-        return 0.0;
     if (globals != NULL) {
         bool changed = false;
         for (int i = 0; i < (int) globalParamNames.size(); i++) {
@@ -2729,22 +2588,6 @@ double OpenCLCalcCustomExternalForceKernel::execute(ContextImpl& context, bool i
         if (changed)
             globals->upload(globalParamValues);
     }
-    if (!hasInitializedKernel) {
-        hasInitializedKernel = true;
-        kernel.setArg<cl_int>(0, numParticles);
-        kernel.setArg<cl::Buffer>(1, cl.getForceBuffers().getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(2, cl.getEnergyBuffer().getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(3, cl.getPosq().getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(4, indices->getDeviceBuffer());
-        int nextIndex = 5;
-        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::Memory>(nextIndex++, buffer.getMemory());
-        }
-    }
-    cl.executeKernel(kernel, numParticles);
     return 0.0;
 }
 

platforms/opencl/src/OpenCLKernels.h

@@ -217,7 +217,7 @@ private:
 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) {
+            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
     }
     ~OpenCLCalcCustomBondForceKernel();
     /**
@@ -242,11 +242,9 @@ private:
     OpenCLContext& cl;
     System& system;
     OpenCLParameterSet* params;
-    OpenCLArray<mm_int4>* indices;
     OpenCLArray<cl_float>* globals;
     std::vector<std::string> globalParamNames;
     std::vector<cl_float> globalParamValues;
-    cl::Kernel kernel;
 };
 
 /**
@@ -288,7 +286,7 @@ private:
 class OpenCLCalcCustomAngleForceKernel : public CalcCustomAngleForceKernel {
 public:
     OpenCLCalcCustomAngleForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomAngleForceKernel(name, platform),
-            hasInitializedKernel(false), cl(cl), system(system), params(NULL), indices(NULL), globals(NULL) {
+            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
     }
     ~OpenCLCalcCustomAngleForceKernel();
     /**
@@ -313,11 +311,9 @@ private:
     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;
 };
 
 /**
@@ -392,7 +388,7 @@ private:
 class OpenCLCalcCMAPTorsionForceKernel : public CalcCMAPTorsionForceKernel {
 public:
     OpenCLCalcCMAPTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCMAPTorsionForceKernel(name, platform),
-            hasInitializedKernel(false), cl(cl), system(system), coefficients(NULL), mapPositions(NULL), torsionIndices(NULL), torsionMaps(NULL) {
+            hasInitializedKernel(false), cl(cl), system(system), coefficients(NULL), mapPositions(NULL), torsionMaps(NULL) {
     }
     ~OpenCLCalcCMAPTorsionForceKernel();
     /**
@@ -418,9 +414,7 @@ private:
     System& system;
     OpenCLArray<mm_float4>* coefficients;
     OpenCLArray<mm_int2>* mapPositions;
-    OpenCLArray<mm_int16>* torsionIndices;
     OpenCLArray<cl_int>* torsionMaps;
-    cl::Kernel kernel;
 };
 
 /**
@@ -429,7 +423,7 @@ private:
 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) {
+            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
     }
     ~OpenCLCalcCustomTorsionForceKernel();
     /**
@@ -454,11 +448,9 @@ private:
     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;
 };
 
 /**
@@ -651,7 +643,7 @@ private:
 class OpenCLCalcCustomExternalForceKernel : public CalcCustomExternalForceKernel {
 public:
     OpenCLCalcCustomExternalForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomExternalForceKernel(name, platform),
-            hasInitializedKernel(false), cl(cl), system(system), params(NULL), indices(NULL), globals(NULL) {
+            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
     }
     ~OpenCLCalcCustomExternalForceKernel();
     /**
@@ -676,11 +668,9 @@ private:
     OpenCLContext& cl;
     System& system;
     OpenCLParameterSet* params;
-    OpenCLArray<cl_int>* indices;
     OpenCLArray<cl_float>* globals;
     std::vector<std::string> globalParamNames;
     std::vector<cl_float> globalParamValues;
-    cl::Kernel kernel;
 };
 
 /**

platforms/opencl/src/kernels/cmapTorsionForce.cl

-/**
- * Compute CNAP torsion forces.
- */
+const float PI = 3.14159265358979323846f;
 
-__kernel void computeCMAPTorsionForces(int numAtoms, int numTorsions, __global float4* forceBuffers, __global float* energyBuffer,
-        __global float4* posq, __global float4* coeff, __global int2* mapPositions, __global int16* indices, __global int* maps) {
-    const float PI = 3.14159265358979323846f;
-    float energy = 0.0f;
-    for (int index = get_global_id(0); index < numTorsions; index += get_global_size(0)) {
-        int16 atoms = indices[index];
-        float4 a1 = posq[atoms.s0];
-        float4 a2 = posq[atoms.s1];
-        float4 a3 = posq[atoms.s2];
-        float4 a4 = posq[atoms.s3];
-        float4 b1 = posq[atoms.s4];
-        float4 b2 = posq[atoms.s5];
-        float4 b3 = posq[atoms.s6];
-        float4 b4 = posq[atoms.s7];
+// Compute the first angle.
 
-        // Compute the first angle.
+float4 v0a = (float4) (pos1.xyz-pos2.xyz, 0.0f);
+float4 v1a = (float4) (pos3.xyz-pos2.xyz, 0.0f);
+float4 v2a = (float4) (pos3.xyz-pos4.xyz, 0.0f);
+float4 cp0a = cross(v0a, v1a);
+float4 cp1a = cross(v1a, v2a);
+float cosangle = dot(normalize(cp0a), normalize(cp1a));
+float angleA;
+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 v0a = (float4) (a1.xyz-a2.xyz, 0.0f);
-        float4 v1a = (float4) (a3.xyz-a2.xyz, 0.0f);
-        float4 v2a = (float4) (a3.xyz-a4.xyz, 0.0f);
-        float4 cp0a = cross(v0a, v1a);
-        float4 cp1a = cross(v1a, v2a);
-        float cosangle = dot(normalize(cp0a), normalize(cp1a));
-        float angleA;
-        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(cp0a, cp1a);
-            float scale = dot(cp0a, cp0a)*dot(cp1a, cp1a);
-            angleA = asin(SQRT(dot(cross_prod, cross_prod)/scale));
-            if (cosangle < 0.0f)
-                angleA = PI-angleA;
-        }
-        else
-           angleA = acos(cosangle);
-        angleA = (dot(v0a, cp1a) >= 0 ? angleA : -angleA);
-        angleA = fmod(angleA+2.0f*PI, 2.0f*PI);
+    float4 cross_prod = cross(cp0a, cp1a);
+    float scale = dot(cp0a, cp0a)*dot(cp1a, cp1a);
+    angleA = asin(SQRT(dot(cross_prod, cross_prod)/scale));
+    if (cosangle < 0.0f)
+        angleA = PI-angleA;
+}
+else
+   angleA = acos(cosangle);
+angleA = (dot(v0a, cp1a) >= 0 ? angleA : -angleA);
+angleA = fmod(angleA+2.0f*PI, 2.0f*PI);
 
-        // Compute the second angle.
+// Compute the second angle.
 
-        float4 v0b = (float4) (b1.xyz-b2.xyz, 0.0f);
-        float4 v1b = (float4) (b3.xyz-b2.xyz, 0.0f);
-        float4 v2b = (float4) (b3.xyz-b4.xyz, 0.0f);
-        float4 cp0b = cross(v0b, v1b);
-        float4 cp1b = cross(v1b, v2b);
-        cosangle = dot(normalize(cp0b), normalize(cp1b));
-        float angleB;
-        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 v0b = (float4) (pos5.xyz-pos6.xyz, 0.0f);
+float4 v1b = (float4) (pos7.xyz-pos6.xyz, 0.0f);
+float4 v2b = (float4) (pos7.xyz-pos8.xyz, 0.0f);
+float4 cp0b = cross(v0b, v1b);
+float4 cp1b = cross(v1b, v2b);
+cosangle = dot(normalize(cp0b), normalize(cp1b));
+float angleB;
+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(cp0b, cp1b);
-            float scale = dot(cp0b, cp0b)*dot(cp1b, cp1b);
-            angleB = asin(SQRT(dot(cross_prod, cross_prod)/scale));
-            if (cosangle < 0.0f)
-                angleB = PI-angleB;
-        }
-        else
-           angleB = acos(cosangle);
-        angleB = (dot(v0b, cp1b) >= 0 ? angleB : -angleB);
-        angleB = fmod(angleB+2.0f*PI, 2.0f*PI);
+    float4 cross_prod = cross(cp0b, cp1b);
+    float scale = dot(cp0b, cp0b)*dot(cp1b, cp1b);
+    angleB = asin(SQRT(dot(cross_prod, cross_prod)/scale));
+    if (cosangle < 0.0f)
+        angleB = PI-angleB;
+}
+else
+   angleB = acos(cosangle);
+angleB = (dot(v0b, cp1b) >= 0 ? angleB : -angleB);
+angleB = fmod(angleB+2.0f*PI, 2.0f*PI);
 
-        // Identify which patch this is in.
+// Identify which patch this is in.
 
-        int2 pos = mapPositions[maps[index]];
-        int size = pos.y;
-        float delta = 2*PI/size;
-        int s = (int) (angleA/delta);
-        int t = (int) (angleB/delta);
-        float4 c[4];
-        int coeffIndex = pos.x+4*(s+size*t);
-        c[0] = coeff[coeffIndex];
-        c[1] = coeff[coeffIndex+1];
-        c[2] = coeff[coeffIndex+2];
-        c[3] = coeff[coeffIndex+3];
-        float da = angleA/delta-s;
-        float db = angleB/delta-t;
+int2 pos = MAP_POS[MAPS[index]];
+int size = pos.y;
+float delta = 2*PI/size;
+int s = (int) (angleA/delta);
+int t = (int) (angleB/delta);
+float4 c[4];
+int coeffIndex = pos.x+4*(s+size*t);
+c[0] = COEFF[coeffIndex];
+c[1] = COEFF[coeffIndex+1];
+c[2] = COEFF[coeffIndex+2];
+c[3] = COEFF[coeffIndex+3];
+float da = angleA/delta-s;
+float db = angleB/delta-t;
 
-        // Evaluate the spline to determine the energy and gradients.
+// Evaluate the spline to determine the energy and gradients.
 
-        float torsionEnergy = 0.0f;
-        float dEdA = 0.0f;
-        float dEdB = 0.0f;
-        torsionEnergy = da*torsionEnergy + ((c[3].w*db + c[3].z)*db + c[3].y)*db + c[3].x;
-        dEdA = db*dEdA + (3.0f*c[3].w*da + 2.0f*c[2].w)*da + c[1].w;
-        dEdB = da*dEdB + (3.0f*c[3].w*db + 2.0f*c[3].z)*db + c[3].y;
-        torsionEnergy = da*torsionEnergy + ((c[2].w*db + c[2].z)*db + c[2].y)*db + c[2].x;
-        dEdA = db*dEdA + (3.0f*c[3].z*da + 2.0f*c[2].z)*da + c[1].z;
-        dEdB = da*dEdB + (3.0f*c[2].w*db + 2.0f*c[2].z)*db + c[2].y;
-        torsionEnergy = da*torsionEnergy + ((c[1].w*db + c[1].z)*db + c[1].y)*db + c[1].x;
-        dEdA = db*dEdA + (3.0f*c[3].y*da + 2.0f*c[2].y)*da + c[1].y;
-        dEdB = da*dEdB + (3.0f*c[1].w*db + 2.0f*c[1].z)*db + c[1].y;
-        torsionEnergy = da*torsionEnergy + ((c[0].w*db + c[0].z)*db + c[0].y)*db + c[0].x;
-        dEdA = db*dEdA + (3.0f*c[3].x*da + 2.0f*c[2].x)*da + c[1].x;
-        dEdB = da*dEdB + (3.0f*c[0].w*db + 2.0f*c[0].z)*db + c[0].y;
-        dEdA /= delta;
-        dEdB /= delta;
-        energy += torsionEnergy;
+float torsionEnergy = 0.0f;
+float dEdA = 0.0f;
+float dEdB = 0.0f;
+torsionEnergy = da*torsionEnergy + ((c[3].w*db + c[3].z)*db + c[3].y)*db + c[3].x;
+dEdA = db*dEdA + (3.0f*c[3].w*da + 2.0f*c[2].w)*da + c[1].w;
+dEdB = da*dEdB + (3.0f*c[3].w*db + 2.0f*c[3].z)*db + c[3].y;
+torsionEnergy = da*torsionEnergy + ((c[2].w*db + c[2].z)*db + c[2].y)*db + c[2].x;
+dEdA = db*dEdA + (3.0f*c[3].z*da + 2.0f*c[2].z)*da + c[1].z;
+dEdB = da*dEdB + (3.0f*c[2].w*db + 2.0f*c[2].z)*db + c[2].y;
+torsionEnergy = da*torsionEnergy + ((c[1].w*db + c[1].z)*db + c[1].y)*db + c[1].x;
+dEdA = db*dEdA + (3.0f*c[3].y*da + 2.0f*c[2].y)*da + c[1].y;
+dEdB = da*dEdB + (3.0f*c[1].w*db + 2.0f*c[1].z)*db + c[1].y;
+torsionEnergy = da*torsionEnergy + ((c[0].w*db + c[0].z)*db + c[0].y)*db + c[0].x;
+dEdA = db*dEdA + (3.0f*c[3].x*da + 2.0f*c[2].x)*da + c[1].x;
+dEdB = da*dEdB + (3.0f*c[0].w*db + 2.0f*c[0].z)*db + c[0].y;
+dEdA /= delta;
+dEdB /= delta;
+energy += torsionEnergy;
 
-        // Apply the force to the first torsion.
+// Apply the force to the first torsion.
 
-        float normCross1 = dot(cp0a, cp0a);
-        float normSqrBC = dot(v1a, v1a);
-        float normBC = SQRT(normSqrBC);
-        float normCross2 = dot(cp1a, cp1a);
-        float dp = 1.0f/normSqrBC;
-        float4 ff = (float4) ((-dEdA*normBC)/normCross1, dot(v0a, v1a)*dp, dot(v2a, v1a)*dp, (dEdA*normBC)/normCross2);
-        float4 internalF0 = ff.x*cp0a;
-        float4 internalF3 = ff.w*cp1a;
-        float4 d = ff.y*internalF0 - ff.z*internalF3;
-        int4 offset = atoms.lo.lo + numAtoms*atoms.hi.lo;
-        float4 forceA = forceBuffers[offset.x];
-        float4 forceB = forceBuffers[offset.y];
-        float4 forceC = forceBuffers[offset.z];
-        float4 forceD = forceBuffers[offset.w];
-        forceA.xyz += internalF0.xyz;
-        forceB.xyz += d.xyz-internalF0.xyz;
-        forceC.xyz += -d.xyz-internalF3.xyz;
-        forceD.xyz += internalF3.xyz;
-        forceBuffers[offset.x] = forceA;
-        forceBuffers[offset.y] = forceB;
-        forceBuffers[offset.z] = forceC;
-        forceBuffers[offset.w] = forceD;
+float normCross1 = dot(cp0a, cp0a);
+float normSqrBC = dot(v1a, v1a);
+float normBC = SQRT(normSqrBC);
+float normCross2 = dot(cp1a, cp1a);
+float dp = 1.0f/normSqrBC;
+float4 ff = (float4) ((-dEdA*normBC)/normCross1, dot(v0a, v1a)*dp, dot(v2a, v1a)*dp, (dEdA*normBC)/normCross2);
+float4 force1 = ff.x*cp0a;
+float4 force4 = ff.w*cp1a;
+float4 d = ff.y*force1 - ff.z*force4;
+float4 force2 = d-force1;
+float4 force3 = -d-force4;
 
-        // Apply the force to the second torsion.
+// Apply the force to the second torsion.
 
-        normCross1 = dot(cp0b, cp0b);
-        normSqrBC = dot(v1b, v1b);
-        normBC = SQRT(normSqrBC);
-        normCross2 = dot(cp1b, cp1b);
-        dp = 1.0f/normSqrBC;
-        ff = (float4) ((-dEdB*normBC)/normCross1, dot(v0b, v1b)*dp, dot(v2b, v1b)*dp, (dEdB*normBC)/normCross2);
-        internalF0 = ff.x*cp0b;
-        internalF3 = ff.w*cp1b;
-        d = ff.y*internalF0 - ff.z*internalF3;
-        offset = atoms.lo.hi + numAtoms*atoms.hi.hi;
-        forceA = forceBuffers[offset.x];
-        forceB = forceBuffers[offset.y];
-        forceC = forceBuffers[offset.z];
-        forceD = forceBuffers[offset.w];
-        forceA.xyz += internalF0.xyz;
-        forceB.xyz += d.xyz-internalF0.xyz;
-        forceC.xyz += -d.xyz-internalF3.xyz;
-        forceD.xyz += internalF3.xyz;
-        forceBuffers[offset.x] = forceA;
-        forceBuffers[offset.y] = forceB;
-        forceBuffers[offset.z] = forceC;
-        forceBuffers[offset.w] = forceD;
-    }
-    energyBuffer[get_global_id(0)] += energy;
-}
+normCross1 = dot(cp0b, cp0b);
+normSqrBC = dot(v1b, v1b);
+normBC = SQRT(normSqrBC);
+normCross2 = dot(cp1b, cp1b);
+dp = 1.0f/normSqrBC;
+ff = (float4) ((-dEdB*normBC)/normCross1, dot(v0b, v1b)*dp, dot(v2b, v1b)*dp, (dEdB*normBC)/normCross2);
+float4 force5 = ff.x*cp0b;
+float4 force8 = ff.w*cp1b;
+d = ff.y*force5 - ff.z*force8;
+float4 force6 = d-force5;
+float4 force7 = -d-force8;

platforms/opencl/src/kernels/customAngleForce.cl

-/**
- * Compute custom angle forces.
- */
-
-__kernel void computeCustomAngleForces(int numAtoms, int numAngles, __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 < numAngles; index += get_global_size(0)) {
-        int8 atoms = indices[index];
-        float4 a1 = posq[atoms.s0];
-        float4 a2 = posq[atoms.s1];
-        float4 a3 = posq[atoms.s2];
-
-        // Compute the force.
-
-        float4 v0 = a2-a1;
-        float4 v1 = a2-a3;
-        float4 cp = cross(v0, v1);
-        float rp = cp.x*cp.x + cp.y*cp.y + cp.z*cp.z;
-        rp = max(sqrt(rp), 1.0e-06f);
-        float r21 = v0.x*v0.x + v0.y*v0.y + v0.z*v0.z;
-        float r23 = v1.x*v1.x + v1.y*v1.y + v1.z*v1.z;
-        float dot = v0.x*v1.x + v0.y*v1.y + v0.z*v1.z;
-        float cosine = clamp(dot/sqrt(r21*r23), -1.0f, 1.0f);
-        float theta = acos(cosine);
-        COMPUTE_FORCE
-        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.
-
-        unsigned int offsetA = atoms.s0+atoms.s3*numAtoms;
-        unsigned int offsetB = atoms.s1+atoms.s4*numAtoms;
-        unsigned int offsetC = atoms.s2+atoms.s5*numAtoms;
-        float4 forceA = forceBuffers[offsetA];
-        float4 forceB = forceBuffers[offsetB];
-        float4 forceC = forceBuffers[offsetC];
-        forceA.xyz += c21.xyz;
-        forceB.xyz -= c21.xyz+c23.xyz;
-        forceC.xyz += c23.xyz;
-        forceBuffers[offsetA] = forceA;
-        forceBuffers[offsetB] = forceB;
-        forceBuffers[offsetC] = forceC;
-    }
-    energyBuffer[get_global_id(0)] += energy;
-}
+float4 v0 = pos2-pos1;
+float4 v1 = pos2-pos3;
+float4 cp = cross(v0, v1);
+float rp = cp.x*cp.x + cp.y*cp.y + cp.z*cp.z;
+rp = max(sqrt(rp), 1.0e-06f);
+float r21 = v0.x*v0.x + v0.y*v0.y + v0.z*v0.z;
+float r23 = v1.x*v1.x + v1.y*v1.y + v1.z*v1.z;
+float dot = v0.x*v1.x + v0.y*v1.y + v0.z*v1.z;
+float cosine = clamp(dot/sqrt(r21*r23), -1.0f, 1.0f);
+float theta = acos(cosine);
+COMPUTE_FORCE
+float4 force1 = cross(v0, cp)*(dEdAngle/(r21*rp));
+float4 force3 = cross(cp, v1)*(dEdAngle/(r23*rp));
+float4 force2 = -force1-force3;

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 int4* indices
-        EXTRA_ARGUMENTS) {
-    float energy = 0.0f;
-    for (int index = get_global_id(0); index < numBonds; index += get_global_size(0)) {
-        // Look up the data for this bond.
-
-        int4 atoms = indices[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;
-    }
-    energyBuffer[get_global_id(0)] += energy;
-}
+float4 delta = pos2-pos1;
+float r = SQRT(delta.x*delta.x + delta.y*delta.y + delta.z*delta.z);
+COMPUTE_FORCE
+delta.xyz *= -dEdR/r;
+float4 force1 = -delta;
+float4 force2 = delta;

platforms/opencl/src/kernels/customExternalForce.cl

-/**
- * Compute custom external forces.
- */
-
-__kernel void computeCustomExternalForces(int numTerms, __global float4* forceBuffers, __global float* energyBuffer,
-        __global float4* posq, __global int* indices
-        EXTRA_ARGUMENTS) {
-    float energy = 0.0f;
-    for (int index = get_global_id(0); index < numTerms; index += get_global_size(0)) {
-        // Look up the data for this particle.
-
-        int atom = indices[index];
-        float4 pos = posq[atom];
-
-        // Compute the force.
-
-        COMPUTE_FORCE
-
-        // Record the force on the atom.
-
-        float4 force = forceBuffers[atom];
-        force.x -= dEdX;
-        force.y -= dEdY;
-        force.z -= dEdZ;
-        forceBuffers[atom] = force;
-    }
-    energyBuffer[get_global_id(0)] += energy;
-}
+COMPUTE_FORCE
+float4 force1 = (float4) (-dEdX, -dEdY, -dEdZ, 0.0f);

platforms/opencl/src/kernels/customTorsionForce.cl

-/**
- * Compute custom torsion forces.
- */
+float4 v0 = (float4) (pos1.xyz-pos2.xyz, 0.0f);
+float4 v1 = (float4) (pos3.xyz-pos2.xyz, 0.0f);
+float4 v2 = (float4) (pos3.xyz-pos4.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.
 
-__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;
+    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 force1 = ff.x*cp0;
+float4 force4 = ff.w*cp1;
+float4 s = ff.y*force1 - ff.z*force4;
+float4 force2 = s-force1;
+float4 force3 = -s-force4;