Created OpenCLBondedUtilities

platforms/opencl/src/OpenCLBondedUtilities.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) 2011 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "OpenCLBondedUtilities.h"
+#include "OpenCLExpressionUtilities.h"
+#include "openmm/OpenMMException.h"
+#include "OpenCLNonbondedUtilities.h"
+#include <iostream>
+
+using namespace OpenMM;
+using namespace std;
+
+OpenCLBondedUtilities::OpenCLBondedUtilities(OpenCLContext& context) : context(context), numForceBuffers(0), maxBonds(0), hasInitializedKernels(false) {
+}
+
+OpenCLBondedUtilities::~OpenCLBondedUtilities() {
+    for (int i = 0; i < (int) atomIndices.size(); i++)
+        delete atomIndices[i];
+    for (int i = 0; i < (int) bufferIndices.size(); i++)
+        delete bufferIndices[i];
+}
+
+void OpenCLBondedUtilities::addInteraction(const vector<vector<int> >& atoms, const string& source) {
+    if (atoms.size() > 0) {
+        forceAtoms.push_back(atoms);
+        forceSource.push_back(source);
+        int width = 1;
+        while (width < atoms[0].size())
+            width *= 2;
+        indexWidth.push_back(width);
+    }
+}
+
+std::string OpenCLBondedUtilities::addArgument(cl::Memory& data, const string& type) {
+    arguments.push_back(&data);
+    argTypes.push_back(type);
+    return "customArg"+OpenCLExpressionUtilities::intToString(arguments.size());
+}
+
+void OpenCLBondedUtilities::initialize(const System& system) {
+    int numForces = forceAtoms.size();
+    if (numForces == 0)
+        return;
+    
+    // Build the lists of atom indicse and buffer indices.
+    
+    vector<vector<cl_uint> > bufferVec(numForces);
+    vector<vector<int> > bufferCounter(numForces, vector<int>(system.getNumParticles(), 0));
+    vector<int> numBuffers(numForces, 0);
+    for (int i = 0; i < numForces; i++) {
+        int numBonds = forceAtoms[i].size();
+        int numAtoms = forceAtoms[i][0].size();
+        int width = indexWidth[i];
+        vector<cl_uint> indexVec(width*numBonds);
+        for (int bond = 0; bond < numBonds; bond++) {
+            for (int atom = 0; atom < numAtoms; atom++)
+                indexVec[bond*width+atom] = forceAtoms[i][bond][atom];
+        }
+        OpenCLArray<cl_uint>* indices = new OpenCLArray<cl_uint>(context, indexVec.size(), "bondedIndices");
+        indices->upload(indexVec);
+        atomIndices.push_back(indices);
+        bufferVec[i].resize(width*numBonds, 0);
+        for (int bond = 0; bond < numBonds; bond++) {
+            for (int atom = 0; atom < numAtoms; atom++)
+                bufferVec[i][bond*width+atom] = bufferCounter[i][forceAtoms[i][bond][atom]]++;
+        }
+        for (int j = 0; j < (int) bufferCounter[i].size(); j++)
+            numBuffers[i] = max(numBuffers[i], bufferCounter[i][j]);
+    }
+    
+    // Figure out how many force buffers will be required.
+    
+    for (int i = 0; i < numForces; i++)
+        numForceBuffers = max(numForceBuffers, numBuffers[i]);
+    int bufferLimit = max(numForceBuffers, (int) context.getPlatformData().contexts.size());
+    if (context.getNonbondedUtilities().getHasInteractions())
+        bufferLimit = max(bufferLimit, context.getNonbondedUtilities().getNumForceBuffers());
+    
+    // For efficiency, we want to merge multiple forces into a single kernel - but only if that
+    // won't increase the number of force buffers.  Figure out sets of forces that can be merged.
+    
+    vector<int> unmerged(numForces);
+    for (int i = 0; i < numForces; i++)
+        unmerged[i] = i;
+    for (int i = 0; i < numForces; i++)
+        for (int j = i-1; j >= 0; j--) {
+            if (numBuffers[unmerged[j]] <= numBuffers[unmerged[j+1]])
+                break;
+            int temp = unmerged[j+1];
+            unmerged[j+1] = unmerged[j];
+            unmerged[j] = temp;
+        }
+    while (unmerged.size() > 0) {
+        int sum = numBuffers[unmerged.back()];
+        int i;
+        for (i = 0; i < unmerged.size()-1; i++) {
+            if (sum+numBuffers[unmerged[i]] > bufferLimit)
+                break;
+            sum += numBuffers[unmerged[i]];
+        }
+        forceSets.push_back(vector<int>());
+        for (int j = 0; j < i; j++)
+            forceSets.back().push_back(unmerged[j]);
+        forceSets.back().push_back(unmerged.back());
+        for (int j = 0; j < i; j++)
+            unmerged.erase(unmerged.begin());
+        unmerged.pop_back();
+    }
+
+    // Update the buffer indices based on merged sets.
+    
+    bufferIndices.resize(numForces);
+    for (int i = 0; i < (int) forceSets.size(); i++)
+        for (int j = 0; j < forceSets[i].size(); j++) {
+            int force = forceSets[i][j];
+            int numBonds = forceAtoms[force].size();
+            int numAtoms = forceAtoms[force][0].size();
+            int width = indexWidth[force];
+            for (int k = 0; k < j; k++)
+                for (int bond = 0; bond < numBonds; bond++)
+                    for (int atom = 0; atom < numAtoms; atom++)
+                        bufferVec[force][bond*width+atom] += bufferCounter[forceSets[i][k]][forceAtoms[force][bond][atom]];
+            OpenCLArray<cl_uint>* buffers = new OpenCLArray<cl_uint>(context, bufferVec[force].size(), "bondedBufferIndices");
+            buffers->upload(bufferVec[force]);
+            bufferIndices[force] = buffers;
+        }
+
+    // Create the kernels.
+
+    for (vector<vector<int> >::const_iterator iter = forceSets.begin(); iter != forceSets.end(); ++iter) {
+        const vector<int>& set = *iter;
+        int setSize = set.size();
+        stringstream s;
+        s<<"__kernel void computeBondedForces(__global float4* forceBuffers, __global float* energyBuffer, __global float4* posq";
+        for (int i = 0; i < setSize; i++) {
+            int force = set[i];
+            string indexType = "uint"+(indexWidth[force] == 1 ? "" : OpenCLExpressionUtilities::intToString(indexWidth[force]));
+            s<<", __global "<<indexType<<"* atomIndices"<<i;
+            s<<", __global "<<indexType<<"* bufferIndices"<<i;
+        }
+        for (int i = 0; i < (int) arguments.size(); i++)
+            s<<", __global "<<argTypes[i]<<"* customArg"<<(i+1);
+        s<<") {\n";
+        s<<"float energy = 0.0f;\n";
+        for (int i = 0; i < setSize; i++) {
+            int force = set[i];
+            s<<createForceSource(i, forceAtoms[force].size(), forceAtoms[force][0].size(), forceSource[force]);
+        }
+        s<<"energyBuffer[get_global_id(0)] += energy;\n";
+        s<<"}\n";
+        map<string, string> defines;
+        defines["PADDED_NUM_ATOMS"] = OpenCLExpressionUtilities::intToString(context.getPaddedNumAtoms());
+        cl::Program program = context.createProgram(s.str(), defines);
+        kernels.push_back(cl::Kernel(program, "computeBondedForces"));
+    }
+    forceAtoms.clear();
+    forceSource.clear();
+}
+
+string OpenCLBondedUtilities::createForceSource(int forceIndex, int numBonds, int numAtoms, const string& computeForce) {
+    maxBonds = max(maxBonds, numBonds);
+    int width = 1;
+    while (width < numAtoms)
+        width *= 2;
+    string suffix1[] = {""};
+    string suffix4[] = {".x", ".y", ".z", ".w"};
+    string suffix16[] = {".s0", ".s1", ".s2", ".s3", ".s4", ".s5", ".s6", ".s7",
+        ".s8", ".s9", ".s10", ".s11", ".s12", ".s13", ".s14", ".s15"};
+    string* suffix;
+    if (width == 1)
+        suffix = suffix1;
+    else if (width <= 4)
+        suffix = suffix4;
+    else
+        suffix = suffix16;
+    string indexType = "uint"+(width == 1 ? "" : OpenCLExpressionUtilities::intToString(width));
+    stringstream s;
+    s<<"for (unsigned int index = get_global_id(0); index < "<<numBonds<<"; index += get_global_size(0)) {\n";
+    s<<"    "<<indexType<<" atoms = atomIndices"<<forceIndex<<"[index];\n";
+    s<<"    "<<indexType<<" buffers = bufferIndices"<<forceIndex<<"[index];\n";
+    for (int i = 0; i < numAtoms; i++) {
+        s<<"    unsigned int atom"<<(i+1)<<" = atoms"<<suffix[i]<<";\n";
+        s<<"    float4 pos"<<(i+1)<<" = posq[atom"<<(i+1)<<"];\n";
+    }
+    s<<computeForce<<"\n";
+    for (int i = 0; i < numAtoms; i++) {
+        s<<"    {\n";
+        s<<"    unsigned int offset = atom"<<(i+1)<<"+buffers"<<suffix[i]<<"*PADDED_NUM_ATOMS;\n";
+        s<<"    float4 force = forceBuffers[offset];\n";
+        s<<"    force.xyz += force"<<(i+1)<<".xyz;\n";
+        s<<"    forceBuffers[offset] = force;\n";
+        s<<"    }\n";
+    }
+    s<<"}\n";
+    return s.str();
+}
+
+void OpenCLBondedUtilities::computeInteractions() {
+    if (!hasInitializedKernels) {
+        hasInitializedKernels = true;
+        for (int i = 0; i < (int) forceSets.size(); i++) {
+            int index = 0;
+            cl::Kernel& kernel = kernels[i];
+            kernel.setArg<cl::Buffer>(index++, context.getForceBuffers().getDeviceBuffer());
+            kernel.setArg<cl::Buffer>(index++, context.getEnergyBuffer().getDeviceBuffer());
+            kernel.setArg<cl::Buffer>(index++, context.getPosq().getDeviceBuffer());
+            for (int j = 0; j < (int) forceSets[i].size(); j++) {
+                kernel.setArg<cl::Buffer>(index++, atomIndices[forceSets[i][j]]->getDeviceBuffer());
+                kernel.setArg<cl::Buffer>(index++, bufferIndices[forceSets[i][j]]->getDeviceBuffer());
+            }
+            for (int j = 0; j < (int) arguments.size(); j++)
+                kernel.setArg<cl::Memory>(index++, *arguments[j]);
+        }
+    }
+    for (int i = 0; i < (int) kernels.size(); i++)
+        context.executeKernel(kernels[i], maxBonds);
+}

platforms/opencl/src/OpenCLBondedUtilities.h

+#ifndef OPENMM_OPENCLBONDEDUTILITIES_H_
+#define OPENMM_OPENCLBONDEDUTILITIES_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2011 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "OpenCLArray.h"
+#include "OpenCLContext.h"
+#include "openmm/System.h"
+#include <string>
+#include <vector>
+
+namespace OpenMM {
+
+/**
+ * This class provides a generic mechanism for evaluating bonded interactions.  You write only
+ * the source code needed to compute one interaction, and this class takes care of creating
+ * and executing a complete kernel that loops over bonds, evaluates each one, and accumulates
+ * the resulting forces and energies.  This offers two advantages.  First, it simplifies the
+ * task of writing a new Force.  Second, it allows multiple forces to be evaluated by a single
+ * kernel, which reduces overhead and improves performance.
+ * 
+ * A "bonded interaction" means an interaction that affects a small, fixed set of particles.
+ * The interaction energy may depend on the positions of only those particles, and the list of
+ * particles forming a "bond" may not change with time.  Examples of bonded interactions
+ * include HarmonicBondForce, HarmonicAngleForce, and PeriodicTorsionForce.
+ * 
+ * To create a bonded interaction, call addInteraction().  You pass to it a block of source
+ * code for evaluating the interaction.  The inputs and outputs for that source code are as
+ * follows:
+ * 
+ * <ol>
+ * <li>The index of the bond being evaluated will have been stored in the unsigned int variable "index".</li>
+ * <li>The indices of the atoms forming that bond will have been stored in the unsigned int variables "atom1",
+ * "atom2", ....</li>
+ * <li>The positions of those atoms will have been stored in the float4 variables "pos1", "pos2", ....</li>
+ * <li>A float variable called "energy" will exist.  Your code should add the potential energy of the
+ * bond to that variable.</li>
+ * <li>Your code should define float4 variables called "force1", "force2", ... that contain the force to
+ * apply to each atom.</li>
+ * </ol>
+ * 
+ * As a simple example, the following source code would be used to implement a pairwise interaction of
+ * the form E=-r^2:
+ * 
+ * <tt><pre>
+ * float4 delta = pos2-pos1;
+ * energy += delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
+ * float4 force1 = 2.0f*delta;
+ * float4 force2 = -2.0f*delta;
+ * </pre></tt>
+ * 
+ * Interactions will often depend on parameters or other data.  Call addArgument() to provide the data
+ * to this class.  It will be passed to the interaction kernel as an argument, and you can refer to it
+ * from your interaction code.
+ */
+
+class OPENMM_EXPORT OpenCLBondedUtilities {
+public:
+    OpenCLBondedUtilities(OpenCLContext& context);
+    ~OpenCLBondedUtilities();
+    /**
+     * Add a bonded interaction.
+     *
+     * @param atoms    this should have one entry for each bond, and that entry should contain the list
+     *                 of atoms involved in the bond.  Every entry must have the same number of atoms.
+     * @param source   the code to evaluate the interaction
+     */
+    void addInteraction(const std::vector<std::vector<int> >& atoms, const std::string& source);
+    /**
+     * Add an argument that should be passed to the interaction kernel.
+     * 
+     * @param data    the device memory containing the data to pass
+     * @param type    the data type contained in the memory (e.g. "float4")
+     * @return the name that will be used for the argument.  Any code you pass to addInteraction() should
+     * refer to it by this name.
+     */
+    std::string addArgument(cl::Memory& data, const std::string& type);
+    /**
+     * Initialize this object in preparation for a simulation.
+     */
+    void initialize(const System& system);
+    /**
+     * Get the number of force buffers required for bonded forces.
+     */
+    int getNumForceBuffers() {
+        return numForceBuffers;
+    }
+    /**
+     * Compute the bonded interactions.
+     */
+    void computeInteractions();
+private:
+    std::string createForceSource(int forceIndex, int numBonds, int numAtoms, const std::string& computeForce);
+    OpenCLContext& context;
+    std::vector<cl::Kernel> kernels;
+    std::vector<std::vector<std::vector<int> > > forceAtoms;
+    std::vector<int> indexWidth;
+    std::vector<std::string> forceSource;
+    std::vector<std::vector<int> > forceSets;
+    std::vector<cl::Memory*> arguments;
+    std::vector<std::string> argTypes;
+    std::vector<OpenCLArray<cl_uint>*> atomIndices;
+    std::vector<OpenCLArray<cl_uint>*> bufferIndices;
+    int numForceBuffers, maxBonds;
+    bool hasInitializedKernels;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_OPENCLBONDEDUTILITIES_H_*/

platforms/opencl/src/OpenCLContext.cpp

@@ -30,6 +30,7 @@
 #include <cmath>
 #include "OpenCLContext.h"
 #include "OpenCLArray.h"
+#include "OpenCLBondedUtilities.h"
 #include "OpenCLForceInfo.h"
 #include "OpenCLIntegrationUtilities.h"
 #include "OpenCLKernelSources.h"
@@ -62,7 +63,7 @@ static void CL_CALLBACK errorCallback(const char* errinfo, const void* private_i
 OpenCLContext::OpenCLContext(int numParticles, int deviceIndex, OpenCLPlatform::PlatformData& platformData) :
         time(0.0), platformData(platformData), stepCount(0), computeForceCount(0), posq(NULL), velm(NULL),
         forceBuffers(NULL), longForceBuffer(NULL), energyBuffer(NULL), atomIndex(NULL), integration(NULL),
-        nonbonded(NULL), thread(NULL) {
+        bonded(NULL), nonbonded(NULL), thread(NULL) {
     try {
         contextIndex = platformData.contexts.size();
         std::vector<cl::Platform> platforms;
@@ -130,6 +131,7 @@ OpenCLContext::OpenCLContext(int numParticles, int deviceIndex, OpenCLPlatform::
         paddedNumAtoms = TileSize*((numParticles+TileSize-1)/TileSize);
         numAtomBlocks = (paddedNumAtoms+(TileSize-1))/TileSize;
         numThreadBlocks = 6*device.getInfo<CL_DEVICE_MAX_COMPUTE_UNITS>();
+        bonded = new OpenCLBondedUtilities(*this);
         nonbonded = new OpenCLNonbondedUtilities(*this);
         posq = new OpenCLArray<mm_float4>(*this, paddedNumAtoms, "posq", true);
         velm = new OpenCLArray<mm_float4>(*this, paddedNumAtoms, "velm", true);
@@ -207,6 +209,8 @@ OpenCLContext::~OpenCLContext() {
         delete atomIndex;
     if (integration != NULL)
         delete integration;
+    if (bonded != NULL)
+        delete bonded;
     if (nonbonded != NULL)
         delete nonbonded;
     if (thread != NULL)
@@ -217,7 +221,9 @@ void OpenCLContext::initialize(const System& system) {
     for (int i = 0; i < numAtoms; i++)
         (*velm)[i].w = (float) (1.0/system.getParticleMass(i));
     velm->upload();
+    bonded->initialize(system);
     numForceBuffers = platformData.contexts.size();
+    numForceBuffers = std::max(numForceBuffers, bonded->getNumForceBuffers());
     for (int i = 0; i < (int) forces.size(); i++)
         numForceBuffers = std::max(numForceBuffers, forces[i]->getRequiredForceBuffers());
     forceBuffers = new OpenCLArray<mm_float4>(*this, paddedNumAtoms*numForceBuffers, "forceBuffers", false);

platforms/opencl/src/OpenCLContext.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2009 Stanford University and the Authors.           *
+ * Portions copyright (c) 2009-2011 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -46,6 +46,7 @@ template <class T>
 class OpenCLArray;
 class OpenCLForceInfo;
 class OpenCLIntegrationUtilities;
+class OpenCLBondedUtilities;
 class OpenCLNonbondedUtilities;
 class System;
 
@@ -425,6 +426,12 @@ public:
     OpenCLIntegrationUtilities& getIntegrationUtilities() {
         return *integration;
     }
+    /**
+     * Get the OpenCLBondedUtilities for this context.
+     */
+    OpenCLBondedUtilities& getBondedUtilities() {
+        return *bonded;
+    }
     /**
      * Get the OpenCLNonbondedUtilities for this context.
      */
@@ -489,6 +496,7 @@ private:
     std::vector<cl::Memory*> autoclearBuffers;
     std::vector<int> autoclearBufferSizes;
     OpenCLIntegrationUtilities* integration;
+    OpenCLBondedUtilities* bonded;
     OpenCLNonbondedUtilities* nonbonded;
     WorkThread* thread;
 };

platforms/opencl/src/OpenCLKernels.cpp

@@ -33,6 +33,7 @@
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/CustomHbondForceImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
+#include "OpenCLBondedUtilities.h"
 #include "OpenCLExpressionUtilities.h"
 #include "OpenCLIntegrationUtilities.h"
 #include "OpenCLNonbondedUtilities.h"
@@ -81,6 +82,7 @@ void OpenCLCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, boo
 }
 
 double OpenCLCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    cl.getBondedUtilities().computeInteractions();
     cl.getNonbondedUtilities().computeInteractions();
     if (includeForces)
         cl.reduceForces();
@@ -228,8 +230,6 @@ private:
 OpenCLCalcHarmonicBondForceKernel::~OpenCLCalcHarmonicBondForceKernel() {
     if (params != NULL)
         delete params;
-    if (indices != NULL)
-        delete indices;
 }
 
 void OpenCLCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
@@ -239,42 +239,22 @@ void OpenCLCalcHarmonicBondForceKernel::initialize(const System& system, const H
     numBonds = endIndex-startIndex;
     if (numBonds == 0)
         return;
+    vector<vector<int> > atoms(numBonds, vector<int>(2));
     params = new OpenCLArray<mm_float2>(cl, numBonds, "bondParams");
-    indices = new OpenCLArray<mm_int4>(cl, numBonds, "bondIndices");
-    vector<int> forceBufferCounter(system.getNumParticles(), 0);
     vector<mm_float2> paramVector(numBonds);
-    vector<mm_int4> indicesVector(numBonds);
     for (int i = 0; i < numBonds; i++) {
-        int particle1, particle2;
         double length, k;
-        force.getBondParameters(startIndex+i, particle1, particle2, length, k);
+        force.getBondParameters(startIndex+i, atoms[i][0], atoms[i][1], length, k);
         paramVector[i] = mm_float2((cl_float) length, (cl_float) k);
-        indicesVector[i] = mm_int4(particle1, particle2, forceBufferCounter[particle1]++, forceBufferCounter[particle2]++);
     }
     params->upload(paramVector);
-    indices->upload(indicesVector);
-    int maxBuffers = 1;
-    for (int i = 0; i < (int) forceBufferCounter.size(); i++)
-        maxBuffers = max(maxBuffers, forceBufferCounter[i]);
-    cl.addForce(new OpenCLBondForceInfo(maxBuffers, force));
-    cl::Program program = cl.createProgram(OpenCLKernelSources::harmonicBondForce);
-    kernel = cl::Kernel(program, "calcHarmonicBondForce");
+    map<string, string> replacements;
+    replacements["PARAMS"] = cl.getBondedUtilities().addArgument(params->getDeviceBuffer(), "float2");
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::harmonicBondForce, replacements));
+    cl.addForce(new OpenCLBondForceInfo(0, force));
 }
 
 double OpenCLCalcHarmonicBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    if (numBonds == 0)
-        return 0.0;
-    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());
-    }
-    cl.executeKernel(kernel, numBonds);
     return 0.0;
 }
 
@@ -457,8 +437,6 @@ private:
 OpenCLCalcHarmonicAngleForceKernel::~OpenCLCalcHarmonicAngleForceKernel() {
     if (params != NULL)
         delete params;
-    if (indices != NULL)
-        delete indices;
 }
 
 void OpenCLCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
@@ -468,44 +446,23 @@ void OpenCLCalcHarmonicAngleForceKernel::initialize(const System& system, const
     numAngles = endIndex-startIndex;
     if (numAngles == 0)
         return;
+    vector<vector<int> > atoms(numAngles, vector<int>(3));
     params = new OpenCLArray<mm_float2>(cl, numAngles, "angleParams");
-    indices = new OpenCLArray<mm_int8>(cl, numAngles, "angleIndices");
-    vector<int> forceBufferCounter(system.getNumParticles(), 0);
     vector<mm_float2> paramVector(numAngles);
-    vector<mm_int8> indicesVector(numAngles);
     for (int i = 0; i < numAngles; i++) {
-        int particle1, particle2, particle3;
         double angle, k;
-        force.getAngleParameters(startIndex+i, particle1, particle2, particle3, angle, k);
+        force.getAngleParameters(startIndex+i, atoms[i][0], atoms[i][1], atoms[i][2], angle, k);
         paramVector[i] = mm_float2((cl_float) angle, (cl_float) k);
-        indicesVector[i] = mm_int8(particle1, particle2, particle3,
-                forceBufferCounter[particle1]++, forceBufferCounter[particle2]++, forceBufferCounter[particle3]++, 0, 0);
 
     }
     params->upload(paramVector);
-    indices->upload(indicesVector);
-    int maxBuffers = 1;
-    for (int i = 0; i < (int) forceBufferCounter.size(); i++)
-        maxBuffers = max(maxBuffers, forceBufferCounter[i]);
-    cl.addForce(new OpenCLAngleForceInfo(maxBuffers, force));
-    cl::Program program = cl.createProgram(OpenCLKernelSources::harmonicAngleForce);
-    kernel = cl::Kernel(program, "calcHarmonicAngleForce");
+    map<string, string> replacements;
+    replacements["PARAMS"] = cl.getBondedUtilities().addArgument(params->getDeviceBuffer(), "float2");
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::harmonicAngleForce, replacements));
+    cl.addForce(new OpenCLAngleForceInfo(0, force));
 }
 
 double OpenCLCalcHarmonicAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    if (numAngles == 0)
-        return 0.0;
-    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, params->getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(6, indices->getDeviceBuffer());
-    }
-    cl.executeKernel(kernel, numAngles);
     return 0.0;
 }
 
@@ -691,8 +648,6 @@ private:
 OpenCLCalcPeriodicTorsionForceKernel::~OpenCLCalcPeriodicTorsionForceKernel() {
     if (params != NULL)
         delete params;
-    if (indices != NULL)
-        delete indices;
 }
 
 void OpenCLCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
@@ -702,44 +657,24 @@ void OpenCLCalcPeriodicTorsionForceKernel::initialize(const System& system, cons
     numTorsions = endIndex-startIndex;
     if (numTorsions == 0)
         return;
+    vector<vector<int> > atoms(numTorsions, vector<int>(4));
     params = new OpenCLArray<mm_float4>(cl, numTorsions, "periodicTorsionParams");
-    indices = new OpenCLArray<mm_int8>(cl, numTorsions, "periodicTorsionIndices");
-    vector<int> forceBufferCounter(system.getNumParticles(), 0);
     vector<mm_float4> paramVector(numTorsions);
-    vector<mm_int8> indicesVector(numTorsions);
     for (int i = 0; i < numTorsions; i++) {
-        int particle1, particle2, particle3, particle4, periodicity;
+        int periodicity;
         double phase, k;
-        force.getTorsionParameters(startIndex+i, particle1, particle2, particle3, particle4, periodicity, phase, k);
+        force.getTorsionParameters(startIndex+i, atoms[i][0], atoms[i][1], atoms[i][2], atoms[i][3], periodicity, phase, k);
         paramVector[i] = mm_float4((cl_float) k, (cl_float) phase, (cl_float) periodicity, 0.0f);
-        indicesVector[i] = mm_int8(particle1, particle2, particle3, particle4,
-                forceBufferCounter[particle1]++, forceBufferCounter[particle2]++, forceBufferCounter[particle3]++, forceBufferCounter[particle4]++);
 
     }
     params->upload(paramVector);
-    indices->upload(indicesVector);
-    int maxBuffers = 1;
-    for (int i = 0; i < (int) forceBufferCounter.size(); i++)
-        maxBuffers = max(maxBuffers, forceBufferCounter[i]);
-    cl.addForce(new OpenCLPeriodicTorsionForceInfo(maxBuffers, force));
-    cl::Program program = cl.createProgram(OpenCLKernelSources::periodicTorsionForce);
-    kernel = cl::Kernel(program, "calcPeriodicTorsionForce");
+    map<string, string> replacements;
+    replacements["PARAMS"] = cl.getBondedUtilities().addArgument(params->getDeviceBuffer(), "float4");
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::periodicTorsionForce, replacements));
+    cl.addForce(new OpenCLPeriodicTorsionForceInfo(0, force));
 }
 
 double OpenCLCalcPeriodicTorsionForceKernel::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, params->getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(6, indices->getDeviceBuffer());
-    }
-    cl.executeKernel(kernel, numTorsions);
     return 0.0;
 }
 
@@ -774,8 +709,6 @@ private:
 OpenCLCalcRBTorsionForceKernel::~OpenCLCalcRBTorsionForceKernel() {
     if (params != NULL)
         delete params;
-    if (indices != NULL)
-        delete indices;
 }
 
 void OpenCLCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
@@ -785,44 +718,23 @@ void OpenCLCalcRBTorsionForceKernel::initialize(const System& system, const RBTo
     numTorsions = endIndex-startIndex;
     if (numTorsions == 0)
         return;
+    vector<vector<int> > atoms(numTorsions, vector<int>(4));
     params = new OpenCLArray<mm_float8>(cl, numTorsions, "rbTorsionParams");
-    indices = new OpenCLArray<mm_int8>(cl, numTorsions, "rbTorsionIndices");
-    vector<int> forceBufferCounter(system.getNumParticles(), 0);
     vector<mm_float8> paramVector(numTorsions);
-    vector<mm_int8> indicesVector(numTorsions);
     for (int i = 0; i < numTorsions; i++) {
-        int particle1, particle2, particle3, particle4;
         double c0, c1, c2, c3, c4, c5;
-        force.getTorsionParameters(startIndex+i, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
+        force.getTorsionParameters(startIndex+i, atoms[i][0], atoms[i][1], atoms[i][2], atoms[i][3], c0, c1, c2, c3, c4, c5);
         paramVector[i] = mm_float8((cl_float) c0, (cl_float) c1, (cl_float) c2, (cl_float) c3, (cl_float) c4, (cl_float) c5, 0.0f, 0.0f);
-        indicesVector[i] = mm_int8(particle1, particle2, particle3, particle4,
-                forceBufferCounter[particle1]++, forceBufferCounter[particle2]++, forceBufferCounter[particle3]++, forceBufferCounter[particle4]++);
 
     }
     params->upload(paramVector);
-    indices->upload(indicesVector);
-    int maxBuffers = 1;
-    for (int i = 0; i < (int) forceBufferCounter.size(); i++)
-        maxBuffers = max(maxBuffers, forceBufferCounter[i]);
-    cl.addForce(new OpenCLRBTorsionForceInfo(maxBuffers, force));
-    cl::Program program = cl.createProgram(OpenCLKernelSources::rbTorsionForce);
-    kernel = cl::Kernel(program, "calcRBTorsionForce");
+    map<string, string> replacements;
+    replacements["PARAMS"] = cl.getBondedUtilities().addArgument(params->getDeviceBuffer(), "float8");
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::rbTorsionForce, replacements));
+    cl.addForce(new OpenCLRBTorsionForceInfo(0, force));
 }
 
 double OpenCLCalcRBTorsionForceKernel::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, params->getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(6, indices->getDeviceBuffer());
-    }
-    cl.executeKernel(kernel, numTorsions);
     return 0.0;
 }
 
@@ -1133,8 +1045,6 @@ OpenCLCalcNonbondedForceKernel::~OpenCLCalcNonbondedForceKernel() {
         delete sigmaEpsilon;
     if (exceptionParams != NULL)
         delete exceptionParams;
-    if (exceptionIndices != NULL)
-        delete exceptionIndices;
     if (cosSinSums != NULL)
         delete cosSinSums;
     if (pmeGrid != NULL)
@@ -1353,44 +1263,27 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
     int startIndex = cl.getContextIndex()*exceptions.size()/numContexts;
     int endIndex = (cl.getContextIndex()+1)*exceptions.size()/numContexts;
     int numExceptions = endIndex-startIndex;
-    int maxBuffers = cl.getNonbondedUtilities().getNumForceBuffers();
     if (numExceptions > 0) {
+        vector<vector<int> > atoms(numExceptions, vector<int>(2));
         exceptionParams = new OpenCLArray<mm_float4>(cl, numExceptions, "exceptionParams");
-        exceptionIndices = new OpenCLArray<mm_int4>(cl, numExceptions, "exceptionIndices");
         vector<mm_float4> exceptionParamsVector(numExceptions);
-        vector<mm_int4> exceptionIndicesVector(numExceptions);
-        vector<int> forceBufferCounter(system.getNumParticles(), 0);
         for (int i = 0; i < numExceptions; i++) {
-            int particle1, particle2;
             double chargeProd, sigma, epsilon;
-            force.getExceptionParameters(exceptions[startIndex+i], particle1, particle2, chargeProd, sigma, epsilon);
+            force.getExceptionParameters(exceptions[startIndex+i], atoms[i][0], atoms[i][1], chargeProd, sigma, epsilon);
             exceptionParamsVector[i] = mm_float4((float) (ONE_4PI_EPS0*chargeProd), (float) sigma, (float) (4.0*epsilon), 0.0f);
-            exceptionIndicesVector[i] = mm_int4(particle1, particle2, forceBufferCounter[particle1]++, forceBufferCounter[particle2]++);
         }
         exceptionParams->upload(exceptionParamsVector);
-        exceptionIndices->upload(exceptionIndicesVector);
-        for (int i = 0; i < (int) forceBufferCounter.size(); i++)
-            maxBuffers = max(maxBuffers, forceBufferCounter[i]);
+        map<string, string> replacements;
+        replacements["PARAMS"] = cl.getBondedUtilities().addArgument(exceptionParams->getDeviceBuffer(), "float4");
+        cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::nonbondedExceptions, replacements));
     }
-    cl.addForce(new OpenCLNonbondedForceInfo(maxBuffers, force));
-    defines.clear();
-    defines["NUM_ATOMS"] = intToString(cl.getPaddedNumAtoms());
-    defines["NUM_EXCEPTIONS"] = intToString(numExceptions);
-    cl::Program program = cl.createProgram(OpenCLKernelSources::nonbondedExceptions, defines);
-    exceptionsKernel = cl::Kernel(program, "computeNonbondedExceptions");
+    cl.addForce(new OpenCLNonbondedForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force));
 }
 
 double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     bool deviceIsCpu = (cl.getDevice().getInfo<CL_DEVICE_TYPE>() == CL_DEVICE_TYPE_CPU);
     if (!hasInitializedKernel) {
         hasInitializedKernel = true;
-        if (exceptionIndices != NULL) {
-            exceptionsKernel.setArg<cl::Buffer>(0, cl.getForceBuffers().getDeviceBuffer());
-            exceptionsKernel.setArg<cl::Buffer>(1, cl.getEnergyBuffer().getDeviceBuffer());
-            exceptionsKernel.setArg<cl::Buffer>(2, cl.getPosq().getDeviceBuffer());
-            exceptionsKernel.setArg<cl::Buffer>(3, exceptionParams->getDeviceBuffer());
-            exceptionsKernel.setArg<cl::Buffer>(4, exceptionIndices->getDeviceBuffer());
-        }
         if (cosSinSums != NULL) {
             ewaldSumsKernel.setArg<cl::Buffer>(0, cl.getEnergyBuffer().getDeviceBuffer());
             ewaldSumsKernel.setArg<cl::Buffer>(1, cl.getPosq().getDeviceBuffer());
@@ -1435,8 +1328,6 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
             }
        }
     }
-    if (exceptionIndices != NULL)
-        cl.executeKernel(exceptionsKernel, exceptionIndices->getSize());
     if (cosSinSums != NULL && cl.getContextIndex() == 0) {
         mm_float4 boxSize = cl.getPeriodicBoxSize();
         mm_float4 recipBoxSize = mm_float4((float) (2*M_PI/boxSize.x), (float) (2*M_PI/boxSize.y), (float) (2*M_PI/boxSize.z), 0);

platforms/opencl/src/OpenCLKernels.h

@@ -184,7 +184,7 @@ private:
 class OpenCLCalcHarmonicBondForceKernel : public CalcHarmonicBondForceKernel {
 public:
     OpenCLCalcHarmonicBondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcHarmonicBondForceKernel(name, platform),
-            hasInitializedKernel(false), cl(cl), system(system), params(NULL), indices(NULL) {
+            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
     }
     ~OpenCLCalcHarmonicBondForceKernel();
     /**
@@ -209,8 +209,6 @@ private:
     OpenCLContext& cl;
     System& system;
     OpenCLArray<mm_float2>* params;
-    OpenCLArray<mm_int4>* indices;
-    cl::Kernel kernel;
 };
 
 /**
@@ -257,7 +255,7 @@ private:
 class OpenCLCalcHarmonicAngleForceKernel : public CalcHarmonicAngleForceKernel {
 public:
     OpenCLCalcHarmonicAngleForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcHarmonicAngleForceKernel(name, platform),
-            hasInitializedKernel(false), cl(cl), system(system), params(NULL), indices(NULL) {
+            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
     }
     ~OpenCLCalcHarmonicAngleForceKernel();
     /**
@@ -282,8 +280,6 @@ private:
     OpenCLContext& cl;
     System& system;
     OpenCLArray<mm_float2>* params;
-    OpenCLArray<mm_int8>* indices;
-    cl::Kernel kernel;
 };
 
 /**
@@ -330,7 +326,7 @@ private:
 class OpenCLCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
 public:
     OpenCLCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcPeriodicTorsionForceKernel(name, platform),
-            hasInitializedKernel(false), cl(cl), system(system), params(NULL), indices(NULL) {
+            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
     }
     ~OpenCLCalcPeriodicTorsionForceKernel();
     /**
@@ -355,8 +351,6 @@ private:
     OpenCLContext& cl;
     System& system;
     OpenCLArray<mm_float4>* params;
-    OpenCLArray<mm_int8>* indices;
-    cl::Kernel kernel;
 };
 
 /**
@@ -365,7 +359,7 @@ private:
 class OpenCLCalcRBTorsionForceKernel : public CalcRBTorsionForceKernel {
 public:
     OpenCLCalcRBTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcRBTorsionForceKernel(name, platform),
-            hasInitializedKernel(false), cl(cl), system(system), params(NULL), indices(NULL) {
+            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
     }
     ~OpenCLCalcRBTorsionForceKernel();
     /**
@@ -390,8 +384,6 @@ private:
     OpenCLContext& cl;
     System& system;
     OpenCLArray<mm_float8>* params;
-    OpenCLArray<mm_int8>* indices;
-    cl::Kernel kernel;
 };
 
 /**
@@ -475,7 +467,7 @@ private:
 class OpenCLCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
 public:
     OpenCLCalcNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcNonbondedForceKernel(name, platform),
-            hasInitializedKernel(false), cl(cl), sigmaEpsilon(NULL), exceptionParams(NULL), exceptionIndices(NULL), cosSinSums(NULL), pmeGrid(NULL),
+            hasInitializedKernel(false), cl(cl), sigmaEpsilon(NULL), exceptionParams(NULL), cosSinSums(NULL), pmeGrid(NULL),
             pmeGrid2(NULL), pmeBsplineModuliX(NULL), pmeBsplineModuliY(NULL), pmeBsplineModuliZ(NULL), pmeBsplineTheta(NULL), pmeAtomRange(NULL),
             pmeAtomGridIndex(NULL), sort(NULL), fft(NULL) {
     }
@@ -501,7 +493,6 @@ private:
     bool hasInitializedKernel;
     OpenCLArray<mm_float2>* sigmaEpsilon;
     OpenCLArray<mm_float4>* exceptionParams;
-    OpenCLArray<mm_int4>* exceptionIndices;
     OpenCLArray<mm_float2>* cosSinSums;
     OpenCLArray<mm_float2>* pmeGrid;
     OpenCLArray<mm_float2>* pmeGrid2;
@@ -513,7 +504,6 @@ private:
     OpenCLArray<mm_int2>* pmeAtomGridIndex;
     OpenCLSort<mm_int2>* sort;
     OpenCLFFT3D* fft;
-    cl::Kernel exceptionsKernel;
     cl::Kernel ewaldSumsKernel;
     cl::Kernel ewaldForcesKernel;
     cl::Kernel pmeGridIndexKernel;

platforms/opencl/src/OpenCLNonbondedUtilities.h

@@ -142,6 +142,12 @@ public:
     double getCutoffDistance() {
         return cutoff;
     }
+    /**
+     * Get whether any interactions have been added.
+     */
+    bool getHasInteractions() {
+        return cutoff != -1.0;
+    }
     /**
      * Prepare to compute interactions.  This updates the neighbor list.
      */

platforms/opencl/src/kernels/harmonicAngleForce.cl

-/**
- * Evaluate the forces due to harmonic angles.
- */
-
-__kernel void calcHarmonicAngleForce(int numAtoms, int numAngles, __global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global float2* params, __global int8* indices) {
-    int index = get_global_id(0);
-    float energy = 0.0f;
-    while (index < numAngles) {
-        // Look up the data for this angle.
-
-        int8 atoms = indices[index];
-        float2 angleParams = params[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*RSQRT(r21*r23), -1.0f, 1.0f);
-        float deltaIdeal = acos(cosine)-angleParams.x;
-        energy += 0.5f*angleParams.y*deltaIdeal*deltaIdeal;
-        float dEdR = angleParams.y*deltaIdeal;
-        float4 c21 = cross(v0, cp)*(dEdR/(r21*rp));
-        float4 c23 = cross(cp, v1)*(dEdR/(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;
-        index += get_global_size(0);
-    }
-    energyBuffer[get_global_id(0)] += energy;
-}
+float2 angleParams = PARAMS[index];
+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*RSQRT(r21*r23), -1.0f, 1.0f);
+float deltaIdeal = acos(cosine)-angleParams.x;
+energy += 0.5f*angleParams.y*deltaIdeal*deltaIdeal;
+float dEdR = angleParams.y*deltaIdeal;
+float4 force1 = cross(v0, cp)*(dEdR/(r21*rp));
+float4 force3 = cross(cp, v1)*(dEdR/(r23*rp));
+float4 force2 = -(force1+force3);

platforms/opencl/src/kernels/harmonicBondForce.cl

-/**
- * Evaluate the forces due to harmonic bonds.
- */
-
-__kernel void calcHarmonicBondForce(int numAtoms, int numBonds, __global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global float2* params, __global int4* indices) {
-    int index = get_global_id(0);
-    float energy = 0.0f;
-    while (index < numBonds) {
-        // Look up the data for this bonds.
-
-        int4 atoms = indices[index];
+float4 delta = pos2-pos1;
+float2 bondParams = PARAMS[index];
+float r = SQRT(delta.x*delta.x + delta.y*delta.y + delta.z*delta.z);
+float deltaIdeal = r-bondParams.x;
+energy += 0.5f * bondParams.y*deltaIdeal*deltaIdeal;
+float dEdR = bondParams.y * deltaIdeal;
+dEdR = (r > 0.0f) ? (dEdR / r) : 0.0f;
+delta.xyz *= dEdR;
+float4 force1 = delta;
+float4 force2 = -delta;
\ No newline at end of file
-        float4 delta = posq[atoms.y]-posq[atoms.x];
-        float2 bondParams = params[index];
-
-        // Compute the force.
-
-        float r = SQRT(delta.x*delta.x + delta.y*delta.y + delta.z*delta.z);
-        float deltaIdeal = r-bondParams.x;
-        energy += 0.5f * bondParams.y*deltaIdeal*deltaIdeal;
-        float dEdR = bondParams.y * deltaIdeal;
-        dEdR = (r > 0.0f) ? (dEdR / r) : 0.0f;
-        delta.xyz *= dEdR;
-
-        // 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/src/kernels/nonbondedExceptions.cl

-/**
- * Compute nonbonded exceptions.
- */
-
-__kernel void computeNonbondedExceptions(__global float4* forceBuffers, __global float* energyBuffer,
-        __global float4* posq, __global float4* params, __global int4* indices) {
-    int index = get_global_id(0);
-    float energy = 0.0f;
-    while (index < NUM_EXCEPTIONS) {
-        // Look up the data for this bonds.
-
-        int4 atoms = indices[index];
-        float4 exceptionParams = params[index];
-        float4 delta = posq[atoms.y]-posq[atoms.x];
-
-        // Compute the force.
-
-        float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
-        float invR = RSQRT(r2);
-        float sig2 = invR*exceptionParams.y;
-        sig2 *= sig2;
-        float sig6 = sig2*sig2*sig2;
-        float dEdR = exceptionParams.z*(12.0f*sig6-6.0f)*sig6;
-        float tempEnergy = exceptionParams.z*(sig6-1.0f)*sig6;
-        dEdR += exceptionParams.x*invR;
-        dEdR *= invR*invR;
-        tempEnergy += exceptionParams.x*invR;
-        energy += tempEnergy;
-        delta.xyz *= dEdR;
-
-        // Record the force on each of the two atoms.
-
-        unsigned int offsetA = atoms.x+atoms.z*NUM_ATOMS;
-        unsigned int offsetB = atoms.y+atoms.w*NUM_ATOMS;
-        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;
-}
+float4 exceptionParams = PARAMS[index];
+float4 delta = pos2-pos1;
+float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
+float invR = RSQRT(r2);
+float sig2 = invR*exceptionParams.y;
+sig2 *= sig2;
+float sig6 = sig2*sig2*sig2;
+float dEdR = exceptionParams.z*(12.0f*sig6-6.0f)*sig6;
+float tempEnergy = exceptionParams.z*(sig6-1.0f)*sig6;
+dEdR += exceptionParams.x*invR;
+dEdR *= invR*invR;
+tempEnergy += exceptionParams.x*invR;
+energy += tempEnergy;
+delta.xyz *= dEdR;
+float4 force1 = -delta;
+float4 force2 = delta;

platforms/opencl/src/kernels/periodicTorsionForce.cl

-/**
- * Evaluate the forces due to periodic torsions.
- */
+const float PI = 3.14159265358979323846f;
+float4 torsionParams = PARAMS[index];
+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 dihedralAngle;
+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 calcPeriodicTorsionForce(int numAtoms, int numTorsions, __global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global float4* params, __global int8* indices) {
-    int index = get_global_id(0);
-    float energy = 0.0f;
-    const float PI = 3.14159265358979323846f;
-    while (index < numTorsions) {
-        // Look up the data for this torsion.
-
-        int8 atoms = indices[index];
-        float4 torsionParams = params[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 dihedralAngle;
-        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);
-            dihedralAngle = asin(sqrt(dot(cross_prod, cross_prod)/scale));
-            if (cosangle < 0.0f)
-                dihedralAngle = PI-dihedralAngle;
-        }
-        else
-           dihedralAngle = acos(cosangle);
-        dihedralAngle = (dot(v0, cp1) >= 0 ? dihedralAngle : -dihedralAngle);
-        float deltaAngle = torsionParams.z*dihedralAngle-torsionParams.y;
-        energy += torsionParams.x*(1.0f+cos(deltaAngle));
-        float sinDeltaAngle = sin(deltaAngle);
-        float dEdAngle = -torsionParams.x*torsionParams.z*sinDeltaAngle;
-        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;
-        index += get_global_size(0);
-    }
-    energyBuffer[get_global_id(0)] += energy;
+    float4 cross_prod = cross(cp0, cp1);
+    float scale = dot(cp0, cp0)*dot(cp1, cp1);
+    dihedralAngle = asin(sqrt(dot(cross_prod, cross_prod)/scale));
+    if (cosangle < 0.0f)
+        dihedralAngle = PI-dihedralAngle;
 }
+else
+   dihedralAngle = acos(cosangle);
+dihedralAngle = (dot(v0, cp1) >= 0 ? dihedralAngle : -dihedralAngle);
+float deltaAngle = torsionParams.z*dihedralAngle-torsionParams.y;
+energy += torsionParams.x*(1.0f+cos(deltaAngle));
+float sinDeltaAngle = sin(deltaAngle);
+float dEdAngle = -torsionParams.x*torsionParams.z*sinDeltaAngle;
+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;
\ No newline at end of file

platforms/opencl/src/kernels/rbTorsionForce.cl

-/**
- * Evaluate the forces due to Ryckaert-Bellemans torsions.
- */
+const float PI = 3.14159265358979323846f;
+float8 torsionParams = PARAMS[index];
+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 dihedralAngle;
+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 calcRBTorsionForce(int numAtoms, int numTorsions, __global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global float8* params, __global int8* indices) {
-    int index = get_global_id(0);
-    float energy = 0.0f;
-    const float PI = 3.14159265358979323846f;
-    while (index < numTorsions) {
-        // Look up the data for this torsion.
-
-        int8 atoms = indices[index];
-        float8 torsionParams = params[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 dihedralAngle;
-        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);
-            dihedralAngle = asin(sqrt(dot(cross_prod, cross_prod)/scale));
-            if (cosangle < 0.0f)
-                dihedralAngle = PI-dihedralAngle;
-        }
-        else
-           dihedralAngle = acos(cosangle);
-        dihedralAngle = (dot(v0, cp1) >= 0 ? dihedralAngle : -dihedralAngle);
-        if (dihedralAngle < 0.0f)
-            dihedralAngle += PI;
-        else
-            dihedralAngle -= PI;
-        cosangle = -cosangle;
-        float cosFactor = cosangle;
-        float dEdAngle = -torsionParams.s1;
-        float rbEnergy = torsionParams.s0;
-        rbEnergy += torsionParams.s1*cosFactor;
-        dEdAngle -= 2.0f*torsionParams.s2*cosFactor;
-        cosFactor *= cosangle;
-        dEdAngle -= 3.0f*torsionParams.s3*cosFactor;
-        rbEnergy += torsionParams.s2*cosFactor;
-        cosFactor *= cosangle;
-        dEdAngle -= 4.0f*torsionParams.s4*cosFactor;
-        rbEnergy += torsionParams.s3*cosFactor;
-        cosFactor *= cosangle;
-        dEdAngle -= 5.0f*torsionParams.s5*cosFactor;
-        rbEnergy += torsionParams.s4*cosFactor;
-        rbEnergy += torsionParams.s5*cosFactor*cosangle;
-        energy += rbEnergy;
-        dEdAngle *= sin(dihedralAngle);
-        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;
-        index += get_global_size(0);
-    }
-    energyBuffer[get_global_id(0)] += energy;
+    float4 cross_prod = cross(cp0, cp1);
+    float scale = dot(cp0, cp0)*dot(cp1, cp1);
+    dihedralAngle = asin(sqrt(dot(cross_prod, cross_prod)/scale));
+    if (cosangle < 0.0f)
+        dihedralAngle = PI-dihedralAngle;
 }
+else
+   dihedralAngle = acos(cosangle);
+dihedralAngle = (dot(v0, cp1) >= 0 ? dihedralAngle : -dihedralAngle);
+if (dihedralAngle < 0.0f)
+    dihedralAngle += PI;
+else
+    dihedralAngle -= PI;
+cosangle = -cosangle;
+float cosFactor = cosangle;
+float dEdAngle = -torsionParams.s1;
+float rbEnergy = torsionParams.s0;
+rbEnergy += torsionParams.s1*cosFactor;
+dEdAngle -= 2.0f*torsionParams.s2*cosFactor;
+cosFactor *= cosangle;
+dEdAngle -= 3.0f*torsionParams.s3*cosFactor;
+rbEnergy += torsionParams.s2*cosFactor;
+cosFactor *= cosangle;
+dEdAngle -= 4.0f*torsionParams.s4*cosFactor;
+rbEnergy += torsionParams.s3*cosFactor;
+cosFactor *= cosangle;
+dEdAngle -= 5.0f*torsionParams.s5*cosFactor;
+rbEnergy += torsionParams.s4*cosFactor;
+rbEnergy += torsionParams.s5*cosFactor*cosangle;
+energy += rbEnergy;
+dEdAngle *= sin(dihedralAngle);
+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;

platforms/opencl/tests/TestOpenCLMultipleForces.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) 2011 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 a system with multiple forces, to make sure OpenCLBondedUtilities is
+ * processing them correctly.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "OpenCLPlatform.h"
+#include "ReferencePlatform.h"
+#include "openmm/HarmonicAngleForce.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/PeriodicTorsionForce.h"
+#include "openmm/RBTorsionForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testForces() {
+    const int numParticles = 100;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    for (int i = 0; i < numParticles-1; i++)
+        bonds->addBond(i, i+1, 1.0, 1.5);
+    system.addForce(bonds);
+    HarmonicAngleForce* angles = new HarmonicAngleForce();
+    for (int i = 0; i < numParticles-2; i++)
+        angles->addAngle(i, i+1, i+2, 2.0, 1.5);
+    system.addForce(angles);
+    PeriodicTorsionForce* periodic = new PeriodicTorsionForce();
+    for (int i = 0; i < numParticles-3; i++)
+        periodic->addTorsion(i, i+1, i+2, i+3, 2, PI_M/3, 1.1);
+    system.addForce(periodic);
+    RBTorsionForce* rb = new RBTorsionForce();
+    for (int i = 0; i < numParticles-3; i += 3)
+        rb->addTorsion(i, i+1, i+2, i+3, 1.0, 1.1, 1.2, 0.3, 0.4, 0.5);
+    system.addForce(rb);
+    ReferencePlatform ref;
+    VerletIntegrator integrator1(0.01);
+    Context context1(system, integrator1, ref);
+    OpenCLPlatform cl;
+    VerletIntegrator integrator2(0.01);
+    Context context2(system, integrator2, cl);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        positions[i] = Vec3(i, genrand_real2(sfmt), genrand_real2(sfmt));
+    context1.setPositions(positions);
+    context2.setPositions(positions);
+    State state1 = context1.getState(State::Forces | State::Energy);
+    State state2 = context1.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces1 = state1.getForces();
+    const vector<Vec3>& forces2 = state2.getForces();
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(forces1[i], forces2[i], TOL);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), TOL);
+}
+
+int main() {
+    try {
+        testForces();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+