Continuing implementation of OpenCL platform.

platforms/opencl/src/OpenCLArray.h

@@ -83,6 +83,12 @@ public:
     int getSize() {
         return size;
     }
+    /**
+     * Get the name of the array.
+     */
+    const std::string& getName() {
+        return name;
+    }
     /**
      * Get the OpenCL Buffer object.
      */
@@ -111,13 +117,13 @@ public:
      * Copy the values in a vector to the Buffer.
      */
     void upload(std::vector<T>& data) {
-        context.getQueue().enqueueWriteBuffer(*buffer, true, 0, size*sizeof(T), &data[0]);
+        context.getQueue().enqueueWriteBuffer(*buffer, CL_TRUE, 0, size*sizeof(T), &data[0]);
     }
     /**
      * Copy the values in the Buffer to a vector.
      */
     void download(std::vector<T>& data) const {
-        context.getQueue().enqueueReadBuffer(*buffer, true, 0, size*sizeof(T), &data[0]);
+        context.getQueue().enqueueReadBuffer(*buffer, CL_TRUE, 0, size*sizeof(T), &data[0]);
     }
     /**
      * Copy the values in the host buffer to the OpenCL Buffer.

platforms/opencl/src/OpenCLContext.cpp

@@ -26,16 +26,19 @@
 
 #include "OpenCLContext.h"
 #include "OpenCLArray.h"
+#include "OpenCLForceInfo.h"
 #include "openmm/Platform.h"
+#include "openmm/System.h"
 #include <fstream>
 #include <iostream>
+#include <sstream>
 
 using namespace OpenMM;
 using namespace std;
 
 OpenCLContext::OpenCLContext(int numParticles, int platformIndex, int deviceIndex) {
     // TODO Select the platform and device correctly
-    context = cl::Context(CL_DEVICE_TYPE_CPU);
+    context = cl::Context(CL_DEVICE_TYPE_GPU);
     device = context.getInfo<CL_CONTEXT_DEVICES>()[0];
     queue = cl::CommandQueue(context, device);
     numAtoms = numParticles;
@@ -43,36 +46,55 @@ OpenCLContext::OpenCLContext(int numParticles, int platformIndex, int deviceInde
     numAtomBlocks = (paddedNumAtoms+(TileSize-1))/TileSize;
     numTiles = numAtomBlocks*(numAtomBlocks+1)/2;
     numThreadBlocks = 8*device.getInfo<CL_DEVICE_MAX_COMPUTE_UNITS>();
-    forceBufferPerWarp = true;
-    numForceBuffers = numThreadBlocks*ThreadBlockSize/TileSize;
-    if (numForceBuffers >= numAtomBlocks) {
-        // For small systems, it is more efficient to have one force buffer per block of 32 atoms instead of one per warp.
-
-        forceBufferPerWarp = false;
-        numForceBuffers = numAtomBlocks;
-    }
-    posq = new OpenCLArray<mm_float4>(*this, paddedNumAtoms, "posq", true);
-    velm = new OpenCLArray<mm_float4>(*this, paddedNumAtoms, "velm", true);
-    forceBuffers = new OpenCLArray<mm_float4>(*this, paddedNumAtoms*numForceBuffers, "forceBuffers", false);
-    force = new OpenCLArray<mm_float4>(*this, &forceBuffers->getDeviceBuffer(), paddedNumAtoms, "force", true);
-    atomIndex = new OpenCLArray<cl_int>(*this, paddedNumAtoms, "atomIndex", true);
-    for (int i = 0; i < paddedNumAtoms; ++i)
-        atomIndex->set(i, i);
-    atomIndex->upload();
 
     // Create utility kernels that are used in multiple places.
 
     utilities = createProgram(loadSourceFromFile("utilities.cl"));
     clearBufferKernel = cl::Kernel(utilities, "clearBuffer");
+    reduceFloat4Kernel = cl::Kernel(utilities, "reduceFloat4Buffer");
 }
 
 OpenCLContext::~OpenCLContext() {
+    for (int i = 0; i < (int) forces.size(); i++)
+        delete forces[i];
     delete posq;
     delete velm;
     delete force;
+    delete forceBuffers;
+    delete energyBuffer;
     delete atomIndex;
 }
 
+void OpenCLContext::initialize(const System& system) {
+//    forceBufferPerWarp = true;
+//    numForceBuffers = numThreadBlocks*ThreadBlockSize/TileSize;
+//    if (numForceBuffers >= numAtomBlocks) {
+//        // For small systems, it is more efficient to have one force buffer per block of 32 atoms instead of one per warp.
+//
+//        forceBufferPerWarp = false;
+//        numForceBuffers = numAtomBlocks;
+//    }
+    posq = new OpenCLArray<mm_float4>(*this, paddedNumAtoms, "posq", true);
+    velm = new OpenCLArray<mm_float4>(*this, paddedNumAtoms, "velm", true);
+    for (int i = 0; i < numAtoms; i++)
+        (*velm)[i].w = system.getParticleMass(i);
+    velm->upload();
+    numForceBuffers = 1;
+    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);
+    force = new OpenCLArray<mm_float4>(*this, &forceBuffers->getDeviceBuffer(), paddedNumAtoms, "force", true);
+    energyBuffer = new OpenCLArray<cl_float>(*this, numThreadBlocks*ThreadBlockSize, "energyBuffer", true);
+    atomIndex = new OpenCLArray<cl_int>(*this, paddedNumAtoms, "atomIndex", true);
+    for (int i = 0; i < paddedNumAtoms; ++i)
+        (*atomIndex)[i] = i;
+    atomIndex->upload();
+}
+
+void OpenCLContext::addForce(OpenCLForceInfo* force) {
+    forces.push_back(force);
+}
+
 string OpenCLContext::loadSourceFromFile(const string& filename) const {
     ifstream file((Platform::getDefaultPluginsDirectory()+"/opencl/"+filename).c_str());
     if (!file.is_open())
@@ -99,14 +121,35 @@ cl::Program OpenCLContext::createProgram(const std::string source) {
     return program;
 }
 
+void OpenCLContext::executeKernel(cl::Kernel& kernel, int workUnits) {
+    int size = std::min((workUnits+ThreadBlockSize-1)/ThreadBlockSize, numThreadBlocks)*ThreadBlockSize;
+    try {
+        queue.enqueueNDRangeKernel(kernel, cl::NullRange, cl::NDRange(size), cl::NDRange(ThreadBlockSize));
+    }
+    catch (cl::Error err) {
+        stringstream str;
+        str<<"Error invoking kernel ";
+        str<<kernel.getInfo<CL_KERNEL_FUNCTION_NAME>()<<": "<<err.what()<<" ("<<err.err()<<")";
+        throw OpenMMException(str.str());
+    }
+}
+
 void OpenCLContext::clearBuffer(OpenCLArray<float>& array) {
     clearBufferKernel.setArg<cl::Buffer>(0, array.getDeviceBuffer());
     clearBufferKernel.setArg<cl_int>(1, array.getSize());
-    queue.enqueueNDRangeKernel(clearBufferKernel, cl::NullRange, cl::NDRange(numThreadBlocks*ThreadBlockSize), cl::NDRange(ThreadBlockSize));
+    executeKernel(clearBufferKernel, array.getSize());
 }
 
 void OpenCLContext::clearBuffer(OpenCLArray<mm_float4>& array) {
     clearBufferKernel.setArg<cl::Buffer>(0, array.getDeviceBuffer());
     clearBufferKernel.setArg<cl_int>(1, array.getSize()*4);
-    queue.enqueueNDRangeKernel(clearBufferKernel, cl::NullRange, cl::NDRange(numThreadBlocks*ThreadBlockSize), cl::NDRange(ThreadBlockSize));
+    executeKernel(clearBufferKernel, array.getSize()*4);
+}
+
+void OpenCLContext::reduceBuffer(OpenCLArray<mm_float4>& array, int numBuffers) {
+    int bufferSize = array.getSize()/numBuffers;
+    reduceFloat4Kernel.setArg<cl::Buffer>(0, array.getDeviceBuffer());
+    reduceFloat4Kernel.setArg<cl_int>(1, bufferSize);
+    reduceFloat4Kernel.setArg<cl_int>(2, numBuffers);
+    executeKernel(reduceFloat4Kernel, bufferSize);
 }

platforms/opencl/src/OpenCLContext.h

@@ -34,16 +34,27 @@ namespace OpenMM {
 
 template <class T>
 class OpenCLArray;
+class OpenCLForceInfo;
+class System;
 
 /**
- * We can't use cl_float4, since different OpenCL implementations currently define it in
- * incompatible ways.  Hopefully that will be fixed in the future.  In the mean time, we
- * define our own type to represent float4 on the host.
+ * We can't use predefined vector types like cl_float4, since different OpenCL implementations currently define
+ * them in incompatible ways.  Hopefully that will be fixed in the future.  In the mean time, we define our own
+ * types to represent them on the host.
  */
 
+typedef struct {
+    cl_float x, y;
+} mm_float2;
 typedef struct {
     cl_float x, y, z, w;
 } mm_float4;
+typedef struct {
+    cl_int x, y;
+} mm_int2;
+typedef struct {
+    cl_int x, y, z, w;
+} mm_int4;
 
 /**
  * This class contains the information associated with a Context by the OpenCL Platform.
@@ -55,6 +66,15 @@ public:
     static const int TileSize = 32;
     OpenCLContext(int numParticles, int platformIndex, int deviceIndex);
     ~OpenCLContext();
+    /**
+     * This is called to initialize internal data structures after all Forces in the system
+     * have been initialized.
+     */
+    void initialize(const System& system);
+    /**
+     * Add an OpenCLForce to this context.
+     */
+    void addForce(OpenCLForceInfo* force);
     /**
      * Get the cl::Context associated with this object.
      */
@@ -91,6 +111,12 @@ public:
     OpenCLArray<mm_float4>& getForceBuffers() {
         return *forceBuffers;
     }
+    /**
+     * Get the array which contains the buffer in which energy is computed.
+     */
+    OpenCLArray<cl_float>& getEnergyBuffer() {
+        return *energyBuffer;
+    }
     /**
      * Get the array which contains the index of each atom.
      */
@@ -105,6 +131,13 @@ public:
      * Create an OpenCL Program from source code.
      */
     cl::Program createProgram(const std::string source);
+    /**
+     * Execute a kernel.
+     *
+     * @param kernel    the kernel to execute
+     * @param workUnits the maximum number of work units that should be used
+     */
+    void executeKernel(cl::Kernel& kernel, int workUnits);
     /**
      * Set all elements of an array to 0.
      */
@@ -113,23 +146,70 @@ public:
      * Set all elements of an array to 0.
      */
     void clearBuffer(OpenCLArray<mm_float4>& array);
+    /**
+     * Given a collection of buffers packed into an array, sum them and store
+     * the sum in the first buffer.
+     *
+     * @param array       the array containing the buffers to reduce
+     * @param numBuffers  the number of buffers packed into the array
+     */
+    void reduceBuffer(OpenCLArray<mm_float4>& array, int numBuffers);
+    /**
+     * Get the number of atoms.
+     */
+    int getNumAtoms() const {
+        return numAtoms;
+    }
+    /**
+     * Get the number of atoms, rounded up to a multiple of TileSize.  This is the actual size of
+     * most arrays with one element per atom.
+     */
+    int getPaddedNumAtoms() const {
+        return paddedNumAtoms;
+    }
+    /**
+     * Get the number of blocks of TileSize atoms.
+     */
+    int getNumAtomBlocks() const {
+        return numAtomBlocks;
+    }
+    /**
+     * Get the standard number of thread blocks to use when executing kernels.
+     */
+    int getNumThreadBlocks() const {
+        return numThreadBlocks;
+    }
+    /**
+     * Get the total number of tiles used for nonbonded computation.
+     */
+    int getNumTiles() const {
+        return numTiles;
+    }
+    /**
+     * Get the number of force buffers.
+     */
+    int getNumForceBuffers() const {
+        return numForceBuffers;
+    }
+private:
     int numAtoms;
     int paddedNumAtoms;
     int numAtomBlocks;
     int numTiles;
     int numThreadBlocks;
     int numForceBuffers;
-    bool forceBufferPerWarp;
-private:
     cl::Context context;
     cl::Device device;
     cl::CommandQueue queue;
     cl::Program utilities;
     cl::Kernel clearBufferKernel;
+    cl::Kernel reduceFloat4Kernel;
+    std::vector<OpenCLForceInfo*> forces;
     OpenCLArray<mm_float4>* posq;
     OpenCLArray<mm_float4>* velm;
     OpenCLArray<mm_float4>* force;
     OpenCLArray<mm_float4>* forceBuffers;
+    OpenCLArray<cl_float>* energyBuffer;
     OpenCLArray<cl_int>* atomIndex;
 };
 

platforms/opencl/src/OpenCLForceInfo.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2009 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "OpenCLForceInfo.h"
+
+using namespace OpenMM;
+using namespace std;
+
+bool OpenCLForceInfo::areParticlesIdentical() {
+    return true;
+}
+
+int OpenCLForceInfo::getNumParticleGroups() {
+    return 0;
+}
+
+void OpenCLForceInfo::getParticlesInGroup(int index, vector<int>& particles) {
+    return;
+}
+
+bool OpenCLForceInfo::areGroupsIdentical(int group1, int group2) {
+    return true;
+}

platforms/opencl/src/OpenCLForceInfo.h

+#ifndef OPENMM_OPENCLFORCEINFO_H_
+#define OPENMM_OPENCLFORCEINFO_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2009 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include <vector>
+
+namespace OpenMM {
+
+/**
+ * This class is used by the OpenCL implementation of a Force class to convey information
+ * about the behavior and requirements of that force.
+ */
+
+class OpenCLForceInfo {
+public:
+    OpenCLForceInfo(int requiredForceBuffers, bool usesNeighborList, double cutoffDistance) :
+        requiredForceBuffers(requiredForceBuffers), usesNeighborList(usesNeighborList), cutoffDistance(cutoffDistance) {
+    }
+    /**
+     * Get the number of force buffers this force requires.
+     */
+    int getRequiredForceBuffers() {
+        return requiredForceBuffers;
+    }
+    /**
+     * Get whether or not this force makes use of the neighbor list.
+     */
+    bool getUsesNeighborList() {
+        return usesNeighborList;
+    }
+    /**
+     * Get the cutoff distance used by this force.  If usesNeighborList() returns false,
+     * the return value from this method is ignored.
+     */
+    double getCutoffDistance() {
+        return cutoffDistance;
+    }
+    /**
+     * Get whether or not two particles have identical force field parameters.
+     */
+    virtual bool areParticlesIdentical();
+    /**
+     * Get the number of particle groups defined by this force.
+     */
+    virtual int getNumParticleGroups();
+    /**
+     * Get the list of particles in a particular group.
+     */
+    virtual void getParticlesInGroup(int index, std::vector<int>& particles);
+    /**
+     * Get whether two particle groups are identical.
+     */
+    virtual bool areGroupsIdentical(int group1, int group2);
+private:
+    int requiredForceBuffers;
+    bool usesNeighborList;
+    double cutoffDistance;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_OPENCLFORCEINFO_H_*/

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -37,8 +37,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
         return new OpenCLCalcForcesAndEnergyKernel(name, platform, data);
     if (name == UpdateStateDataKernel::Name())
         return new OpenCLUpdateStateDataKernel(name, platform, data);
-//    if (name == CalcHarmonicBondForceKernel::Name())
-//        return new OpenCLCalcHarmonicBondForceKernel(name, platform, data, context.getSystem());
+    if (name == CalcHarmonicBondForceKernel::Name())
+        return new OpenCLCalcHarmonicBondForceKernel(name, platform, data, context.getSystem());
 //    if (name == CalcHarmonicAngleForceKernel::Name())
 //        return new OpenCLCalcHarmonicAngleForceKernel(name, platform, data, context.getSystem());
 //    if (name == CalcPeriodicTorsionForceKernel::Name())

platforms/opencl/src/OpenCLKernels.cpp

@@ -24,9 +24,8 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
  * -------------------------------------------------------------------------- */
 
-#include "OpenCLArray.h"
 #include "OpenCLKernels.h"
-#include "OpenCLContext.h"
+#include "OpenCLForceInfo.h"
 #include "openmm/LangevinIntegrator.h"
 #include "openmm/Context.h"
 #include "openmm/internal/ContextImpl.h"
@@ -43,13 +42,20 @@ void OpenCLCalcForcesAndEnergyKernel::beginForceComputation(ContextImpl& context
 }
 
 void OpenCLCalcForcesAndEnergyKernel::finishForceComputation(ContextImpl& context) {
+    data.context->reduceBuffer(data.context->getForceBuffers(), data.context->getNumForceBuffers());
 }
 
 void OpenCLCalcForcesAndEnergyKernel::beginEnergyComputation(ContextImpl& context) {
+    data.context->clearBuffer(data.context->getEnergyBuffer());
 }
 
 double OpenCLCalcForcesAndEnergyKernel::finishEnergyComputation(ContextImpl& context) {
-    return 0.0;
+    OpenCLArray<cl_float>& energy = data.context->getEnergyBuffer();
+    energy.download();
+    double sum = 0.0f;
+    for (int i = 0; i < energy.getSize(); i++)
+        sum += energy[i];
+    return sum;
 }
 
 void OpenCLUpdateStateDataKernel::initialize(const System& system) {
@@ -131,37 +137,80 @@ void OpenCLUpdateStateDataKernel::getForces(ContextImpl& context, std::vector<Ve
     }
 }
 
-//OpenCLCalcHarmonicBondForceKernel::~OpenCLCalcHarmonicBondForceKernel() {
-//}
-//
-//void OpenCLCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
-//    if (data.primaryKernel == NULL)
-//        data.primaryKernel = this;
-//    data.hasBonds = true;
-//    numBonds = force.getNumBonds();
-//    vector<int> particle1(numBonds);
-//    vector<int> particle2(numBonds);
-//    vector<float> length(numBonds);
-//    vector<float> k(numBonds);
-//    for (int i = 0; i < numBonds; i++) {
-//        double lengthValue, kValue;
-//        force.getBondParameters(i, particle1[i], particle2[i], lengthValue, kValue);
-//        length[i] = (float) lengthValue;
-//        k[i] = (float) kValue;
-//    }
-//    gpuSetBondParameters(data.gpu, particle1, particle2, length, k);
-//}
-//
-//void OpenCLCalcHarmonicBondForceKernel::executeForces(ContextImpl& context) {
-//    if (data.primaryKernel == this)
-//        calcForces(context, data);
-//}
-//
-//double OpenCLCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
-//    if (data.primaryKernel == this)
-//        return calcEnergy(context, data, system);
-//    return 0.0;
-//}
+OpenCLCalcHarmonicBondForceKernel::~OpenCLCalcHarmonicBondForceKernel() {
+    if (params != NULL)
+        delete params;
+    if (indices != NULL)
+        delete indices;
+}
+
+class OpenCLBondForceInfo : public OpenCLForceInfo {
+public:
+    OpenCLBondForceInfo(int requiredBuffers, const HarmonicBondForce& force) : OpenCLForceInfo(requiredBuffers, false, 0.0), force(force) {
+    }
+    int getNumParticleGroups() {
+        return force.getNumBonds();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2;
+        double length, k;
+        force.getBondParameters(index, particle1, particle2, length, k);
+        particles.resize(2);
+        particles[0] = particle1;
+        particles[1] = particle2;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2;
+        double length1, length2, k1, k2;
+        force.getBondParameters(group1, particle1, particle2, length1, k1);
+        force.getBondParameters(group2, particle1, particle2, length2, k2);
+        return (length1 == length2 && k1 == k2);
+    }
+private:
+    const HarmonicBondForce& force;
+};
+
+void OpenCLCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
+    numBonds = force.getNumBonds();
+    params = new OpenCLArray<mm_float2>(*data.context, numBonds, "BondParams");
+    indices = new OpenCLArray<mm_int4>(*data.context, 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(i, particle1, particle2, length, k);
+        paramVector[i] = (mm_float2) {length, 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 < forceBufferCounter.size(); i++) {
+        maxBuffers = max(maxBuffers, forceBufferCounter[i]);
+    }
+    data.context->addForce(new OpenCLBondForceInfo(maxBuffers, force));
+    cl::Program program = data.context->createProgram(data.context->loadSourceFromFile("harmonicBondForce.cl"));
+    kernel = cl::Kernel(program, "calcHarmonicBondForce");
+}
+
+void OpenCLCalcHarmonicBondForceKernel::executeForces(ContextImpl& context) {
+    kernel.setArg<cl_int>(0, data.context->getPaddedNumAtoms());
+    kernel.setArg<cl_int>(1, numBonds);
+    kernel.setArg<cl::Buffer>(2, data.context->getForceBuffers().getDeviceBuffer());
+    kernel.setArg<cl::Buffer>(3, data.context->getEnergyBuffer().getDeviceBuffer());
+    kernel.setArg<cl::Buffer>(4, data.context->getPosq().getDeviceBuffer());
+    kernel.setArg<cl::Buffer>(5, params->getDeviceBuffer());
+    kernel.setArg<cl::Buffer>(6, indices->getDeviceBuffer());
+    data.context->executeKernel(kernel, numBonds);
+}
+
+double OpenCLCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
+    executeForces(context);
+    return 0.0;
+}
 //
 //OpenCLCalcHarmonicAngleForceKernel::~OpenCLCalcHarmonicAngleForceKernel() {
 //}
@@ -625,6 +674,7 @@ OpenCLIntegrateVerletStepKernel::~OpenCLIntegrateVerletStepKernel() {
 
 void OpenCLIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
 //    initializeIntegration(system, data, integrator);
+    data.context->initialize(system);
     prevStepSize = -1.0;
 }
 

platforms/opencl/src/OpenCLKernels.h

@@ -28,6 +28,8 @@
  * -------------------------------------------------------------------------- */
 
 #include "OpenCLPlatform.h"
+#include "OpenCLArray.h"
+#include "OpenCLContext.h"
 #include "openmm/kernels.h"
 #include "openmm/System.h"
 
@@ -142,41 +144,45 @@ public:
 private:
     OpenCLPlatform::PlatformData& data;
 };
-//
-///**
-// * This kernel is invoked by HarmonicBondForce to calculate the forces acting on the system and the energy of the system.
-// */
-//class OpenCLCalcHarmonicBondForceKernel : public CalcHarmonicBondForceKernel {
-//public:
-//    OpenCLCalcHarmonicBondForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcHarmonicBondForceKernel(name, platform), data(data), system(system) {
-//    }
-//    ~OpenCLCalcHarmonicBondForceKernel();
-//    /**
-//     * Initialize the kernel.
-//     *
-//     * @param system     the System this kernel will be applied to
-//     * @param force      the HarmonicBondForce this kernel will be used for
-//     */
-//    void initialize(const System& system, const HarmonicBondForce& force);
-//    /**
-//     * Execute the kernel to calculate the forces.
-//     *
-//     * @param context    the context in which to execute this kernel
-//     */
-//    void executeForces(ContextImpl& context);
-//    /**
-//     * Execute the kernel to calculate the energy.
-//     *
-//     * @param context    the context in which to execute this kernel
-//     * @return the potential energy due to the HarmonicBondForce
-//     */
-//    double executeEnergy(ContextImpl& context);
-//private:
-//    int numBonds;
-//    OpenCLPlatform::PlatformData& data;
-//    System& system;
-//};
-//
+
+/**
+ * This kernel is invoked by HarmonicBondForce to calculate the forces acting on the system and the energy of the system.
+ */
+class OpenCLCalcHarmonicBondForceKernel : public CalcHarmonicBondForceKernel {
+public:
+    OpenCLCalcHarmonicBondForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) :
+        CalcHarmonicBondForceKernel(name, platform), data(data), system(system), params(NULL), indices(NULL) {
+    }
+    ~OpenCLCalcHarmonicBondForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the HarmonicBondForce this kernel will be used for
+     */
+    void initialize(const System& system, const HarmonicBondForce& force);
+    /**
+     * Execute the kernel to calculate the forces.
+     *
+     * @param context    the context in which to execute this kernel
+     */
+    void executeForces(ContextImpl& context);
+    /**
+     * Execute the kernel to calculate the energy.
+     *
+     * @param context    the context in which to execute this kernel
+     * @return the potential energy due to the HarmonicBondForce
+     */
+    double executeEnergy(ContextImpl& context);
+private:
+    int numBonds;
+    OpenCLPlatform::PlatformData& data;
+    System& system;
+    OpenCLArray<mm_float2>* params;
+    OpenCLArray<mm_int4>* indices;
+    cl::Kernel kernel;
+};
+
 ///**
 // * This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
 // */

platforms/opencl/src/OpenCLPlatform.cpp

@@ -47,7 +47,7 @@ OpenCLPlatform::OpenCLPlatform() {
     OpenCLKernelFactory* factory = new OpenCLKernelFactory();
     registerKernelFactory(CalcForcesAndEnergyKernel::Name(), factory);
     registerKernelFactory(UpdateStateDataKernel::Name(), factory);
-//    registerKernelFactory(CalcHarmonicBondForceKernel::Name(), factory);
+    registerKernelFactory(CalcHarmonicBondForceKernel::Name(), factory);
 //    registerKernelFactory(CalcHarmonicAngleForceKernel::Name(), factory);
 //    registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
 //    registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);

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) {
+        int4 atoms = indices[index];
+        float4 delta = posq[atoms.y]-posq[atoms.x];
+        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;
+        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/utilities.cl

+/**
+ * Fill a buffer with 0.
+ */
+
 __kernel void clearBuffer(__global float* buffer, int size) {
     int index = get_global_id(0);
     int step = get_global_size(0);
@@ -11,3 +15,20 @@ __kernel void clearBuffer(__global float* buffer, int size) {
         for (int i = sizeDiv4*4; i < size; i++)
             buffer[i] = 0.0f;
 }
+
+/**
+ * Sum a collection of buffers into the first one.
+ */
+
+__kernel void reduceFloat4Buffer(__global float4* buffer, int bufferSize, int numBuffers) {
+    int index = get_global_id(0);
+    int step = get_global_size(0);
+    int totalSize = bufferSize*numBuffers;
+    while (index < bufferSize) {
+        float4 sum = buffer[index];
+        for (int i = index+bufferSize; i < totalSize; i += bufferSize)
+            sum += buffer[i];
+        buffer[index] = sum;
+        index += step;
+    }
+}

platforms/opencl/tests/TestOpenCLHarmonicBondForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the OpenCL implementation of HarmonicBondForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "OpenCLPlatform.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testBonds() {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    HarmonicBondForce* forceField = new HarmonicBondForce();
+    forceField->addBond(0, 1, 1.5, 0.8);
+    forceField->addBond(1, 2, 1.2, 0.7);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 2, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(0, -0.8*0.5, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0.7*0.2, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-forces[0][0]-forces[2][0], -forces[0][1]-forces[2][1], -forces[0][2]-forces[2][2]), forces[1], TOL);
+    ASSERT_EQUAL_TOL(0.5*0.8*0.5*0.5 + 0.5*0.7*0.2*0.2, state.getPotentialEnergy(), TOL);
+}
+
+int main() {
+    try {
+        testBonds();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+