Created OpenCL implementations of VerletIntegrator and SHAKE.

platforms/opencl/src/OpenCLArray.h

@@ -29,6 +29,8 @@
 
 #include "OpenCLContext.h"
 #include "openmm/OpenMMException.h"
+#include <iostream>
+#include <sstream>
 #include <vector>
 
 namespace OpenMM {
@@ -117,13 +119,27 @@ public:
      * Copy the values in a vector to the Buffer.
      */
     void upload(std::vector<T>& data) {
-        context.getQueue().enqueueWriteBuffer(*buffer, CL_TRUE, 0, size*sizeof(T), &data[0]);
+        try {
+            context.getQueue().enqueueWriteBuffer(*buffer, CL_TRUE, 0, size*sizeof(T), &data[0]);
+        }
+        catch (cl::Error err) {
+            std::stringstream str;
+            str<<"Error uploading array "<<name<<": "<<err.what()<<" ("<<err.err()<<")";
+            throw OpenMMException(str.str());
+        }
     }
     /**
      * Copy the values in the Buffer to a vector.
      */
     void download(std::vector<T>& data) const {
-        context.getQueue().enqueueReadBuffer(*buffer, CL_TRUE, 0, size*sizeof(T), &data[0]);
+        try {
+            context.getQueue().enqueueReadBuffer(*buffer, CL_TRUE, 0, size*sizeof(T), &data[0]);
+        }
+        catch (cl::Error err) {
+            std::stringstream str;
+            str<<"Error downloading array "<<name<<": "<<err.what()<<" ("<<err.err()<<")";
+            throw OpenMMException(str.str());
+        }
     }
     /**
      * Copy the values in the host buffer to the OpenCL Buffer.

platforms/opencl/src/OpenCLContext.cpp

@@ -27,6 +27,7 @@
 #include "OpenCLContext.h"
 #include "OpenCLArray.h"
 #include "OpenCLForceInfo.h"
+#include "OpenCLIntegrationUtilities.h"
 #include "openmm/Platform.h"
 #include "openmm/System.h"
 #include <fstream>
@@ -36,7 +37,7 @@
 using namespace OpenMM;
 using namespace std;
 
-OpenCLContext::OpenCLContext(int numParticles, int deviceIndex) {
+OpenCLContext::OpenCLContext(int numParticles, int deviceIndex) : time(0.0), stepCount(0) {
     context = cl::Context(CL_DEVICE_TYPE_ALL);
     vector<cl::Device> devices = context.getInfo<CL_CONTEXT_DEVICES>();
     const int minThreadBlockSize = 32;
@@ -79,6 +80,7 @@ OpenCLContext::~OpenCLContext() {
     delete forceBuffers;
     delete energyBuffer;
     delete atomIndex;
+    delete integration;
 }
 
 void OpenCLContext::initialize(const System& system) {
@@ -93,7 +95,7 @@ void OpenCLContext::initialize(const System& system) {
     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)[i].w = (float) (1.0/system.getParticleMass(i));
     velm->upload();
     numForceBuffers = 1;
     for (int i = 0; i < (int) forces.size(); i++)
@@ -105,6 +107,7 @@ void OpenCLContext::initialize(const System& system) {
     for (int i = 0; i < paddedNumAtoms; ++i)
         (*atomIndex)[i] = i;
     atomIndex->upload();
+    integration = new OpenCLIntegrationUtilities(*this, system);
 }
 
 void OpenCLContext::addForce(OpenCLForceInfo* force) {
@@ -144,8 +147,7 @@ void OpenCLContext::executeKernel(cl::Kernel& kernel, int workUnits) {
     }
     catch (cl::Error err) {
         stringstream str;
-        str<<"Error invoking kernel ";
-        str<<kernel.getInfo<CL_KERNEL_FUNCTION_NAME>()<<": "<<err.what()<<" ("<<err.err()<<")";
+        str<<"Error invoking kernel "<<kernel.getInfo<CL_KERNEL_FUNCTION_NAME>()<<": "<<err.what()<<" ("<<err.err()<<")";
         throw OpenMMException(str.str());
     }
 }

platforms/opencl/src/OpenCLContext.h

@@ -35,6 +35,7 @@ namespace OpenMM {
 template <class T>
 class OpenCLArray;
 class OpenCLForceInfo;
+class OpenCLIntegrationUtilities;
 class System;
 
 /**
@@ -100,7 +101,7 @@ public:
         return *posq;
     }
     /**
-     * Get the array which contains the velocity and massof each atom.
+     * Get the array which contains the velocity and inverse mass of each atom.
      */
     OpenCLArray<mm_float4>& getVelm() {
         return *velm;
@@ -160,6 +161,30 @@ public:
      * @param numBuffers  the number of buffers packed into the array
      */
     void reduceBuffer(OpenCLArray<mm_float4>& array, int numBuffers);
+    /**
+     * Get the current simulation time.
+     */
+    double getTime() {
+        return time;
+    }
+    /**
+     * Set the current simulation time.
+     */
+    void setTime(double t) {
+        time = t;
+    }
+    /**
+     * Get the number of integration steps that have been taken.
+     */
+    int getStepCount() {
+        return stepCount;
+    }
+    /**
+     * Set the number of integration steps that have been taken.
+     */
+    void setStepCount(int steps) {
+        stepCount = steps;;
+    }
     /**
      * Get the number of atoms.
      */
@@ -197,7 +222,15 @@ public:
     int getNumForceBuffers() const {
         return numForceBuffers;
     }
+    /**
+     * Get the OpenCLIntegrationUtilities for this context.
+     */
+    OpenCLIntegrationUtilities& getIntegrationUtilties() {
+        return *integration;
+    }
 private:
+    double time;
+    int stepCount;
     int numAtoms;
     int paddedNumAtoms;
     int numAtomBlocks;
@@ -217,6 +250,7 @@ private:
     OpenCLArray<mm_float4>* forceBuffers;
     OpenCLArray<cl_float>* energyBuffer;
     OpenCLArray<cl_int>* atomIndex;
+    OpenCLIntegrationUtilities* integration;
 };
 
 } // namespace OpenMM

platforms/opencl/src/OpenCLIntegrationUtilities.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 "OpenCLIntegrationUtilities.h"
+#include "OpenCLArray.h"
+#include <map>
+
+using namespace OpenMM;
+using namespace std;
+
+struct OpenCLIntegrationUtilities::ShakeCluster {
+    int centralID;
+    int peripheralID[3];
+    int size;
+    bool valid;
+    double distance;
+    double centralInvMass, peripheralInvMass;
+    ShakeCluster() : valid(true) {
+    }
+    ShakeCluster(int centralID, double invMass) : centralID(centralID), centralInvMass(invMass), size(0), valid(true) {
+    }
+    void addAtom(int id, double dist, double invMass) {
+        if (size == 3 || (size > 0 && dist != distance) || (size > 0 && invMass != peripheralInvMass))
+            valid = false;
+        else {
+            peripheralID[size++] = id;
+            distance = dist;
+            peripheralInvMass = invMass;
+        }
+    }
+    void markInvalid(map<int, ShakeCluster>& allClusters, vector<bool>& invalidForShake)
+    {
+        valid = false;
+        invalidForShake[centralID] = true;
+        for (int i = 0; i < size; i++) {
+            invalidForShake[peripheralID[i]] = true;
+            map<int, ShakeCluster>::iterator otherCluster = allClusters.find(peripheralID[i]);
+            if (otherCluster != allClusters.end() && otherCluster->second.valid)
+                otherCluster->second.markInvalid(allClusters, invalidForShake);
+        }
+    }
+};
+
+OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, const System& system) : context(context),
+        oldPos(NULL), posDelta(NULL), settleAtoms(NULL), settleParams(NULL), shakeAtoms(NULL), shakeParams(NULL) {
+    // Create workspace arrays.
+
+    posDelta = new OpenCLArray<mm_float4>(context, context.getPaddedNumAtoms(), "posDelta");
+    oldPos = new OpenCLArray<mm_float4>(context, context.getPaddedNumAtoms(), "oldPos");
+
+    // Create kernels for enforcing constraints.
+
+    cl::Program program = context.createProgram(context.loadSourceFromFile("shakeHydrogens.cl"));
+    shakeKernel = cl::Kernel(program, "applyShakeToHydrogens");
+
+    // Record the set of constraints and how many constraints each atom is involved in.
+
+    int numConstraints = system.getNumConstraints();
+    vector<int> atom1(numConstraints);
+    vector<int> atom2(numConstraints);
+    vector<double> distance(numConstraints);
+    vector<int> constraintCount(context.getNumAtoms(), 0);
+    for (int i = 0; i < numConstraints; i++) {
+        system.getConstraintParameters(i, atom1[i], atom2[i], distance[i]);
+        constraintCount[atom1[i]]++;
+        constraintCount[atom2[i]]++;
+    }
+
+    // Identify clusters of three atoms that can be treated with SETTLE.  First, for every
+    // atom that might be part of such a cluster, make a list of the two other atoms it is
+    // connected to.
+
+    vector<map<int, float> > settleConstraints(system.getNumParticles());
+    for (int i = 0; i < (int)atom1.size(); i++) {
+        if (constraintCount[atom1[i]] == 2 && constraintCount[atom2[i]] == 2) {
+            settleConstraints[atom1[i]][atom2[i]] = distance[i];
+            settleConstraints[atom2[i]][atom1[i]] = distance[i];
+        }
+    }
+
+    // Now remove the ones that don't actually form closed loops of three atoms.
+
+    vector<int> settleClusters;
+    for (int i = 0; i < (int)settleConstraints.size(); i++) {
+        if (settleConstraints[i].size() == 2) {
+            int partner1 = settleConstraints[i].begin()->first;
+            int partner2 = (++settleConstraints[i].begin())->first;
+            if (settleConstraints[partner1].size() != 2 || settleConstraints[partner2].size() != 2 ||
+                    settleConstraints[partner1].find(partner2) == settleConstraints[partner1].end())
+                settleConstraints[i].clear();
+            else if (i < partner1 && i < partner2)
+                settleClusters.push_back(i);
+        }
+        else
+            settleConstraints[i].clear();
+    }
+
+    // Record the SETTLE clusters.
+
+    vector<bool> isShakeAtom(system.getNumParticles(), false);
+    if (settleClusters.size() > 0) {
+        vector<mm_int4> atoms(settleAtoms->getSize());
+        vector<mm_float2> params(settleParams->getSize());
+        for (int i = 0; i < (int) settleClusters.size(); i++) {
+            int atom1 = settleClusters[i];
+            int atom2 = settleConstraints[atom1].begin()->first;
+            int atom3 = (++settleConstraints[atom1].begin())->first;
+            float dist12 = settleConstraints[atom1].find(atom2)->second;
+            float dist13 = settleConstraints[atom1].find(atom3)->second;
+            float dist23 = settleConstraints[atom2].find(atom3)->second;
+            if (dist12 == dist13) {
+                // atom1 is the central atom
+                atoms.push_back((mm_int4) {atom1, atom2, atom3, 0});
+                params.push_back((mm_float2) {dist12, dist23});
+            }
+            else if (dist12 == dist23) {
+                // atom2 is the central atom
+                atoms.push_back((mm_int4) {atom2, atom1, atom3, 0});
+                params.push_back((mm_float2) {dist12, dist13});
+            }
+            else if (dist13 == dist23) {
+                // atom3 is the central atom
+                atoms.push_back((mm_int4) {atom3, atom1, atom2, 0});
+                params.push_back((mm_float2) {dist13, dist12});
+            }
+            else
+                throw OpenMMException("Two of the three distances constrained with SETTLE must be the same.");
+            isShakeAtom[atom1] = true;
+            isShakeAtom[atom2] = true;
+            isShakeAtom[atom3] = true;
+        }
+        settleAtoms = new OpenCLArray<mm_int4>(context, atoms.size(), "settleAtoms");
+        settleParams = new OpenCLArray<mm_float2>(context, params.size(), "settleParams");
+        settleAtoms->upload(atoms);
+        settleParams->upload(params);
+    }
+
+    // Find clusters consisting of a central atom with up to three peripheral atoms.
+
+    map<int, ShakeCluster> clusters;
+    vector<bool> invalidForShake(system.getNumParticles(), false);
+    for (int i = 0; i < (int) atom1.size(); i++) {
+        if (isShakeAtom[atom1[i]])
+            continue; // This is being taken care of with SETTLE.
+
+        // Determine which is the central atom.
+
+        bool firstIsCentral;
+        if (constraintCount[atom1[i]] > 1)
+            firstIsCentral = true;
+        else if (constraintCount[atom2[i]] > 1)
+            firstIsCentral = false;
+        else if (atom1[i] < atom2[i])
+            firstIsCentral = true;
+        else
+            firstIsCentral = false;
+        int centralID, peripheralID;
+        if (firstIsCentral) {
+            centralID = atom1[i];
+            peripheralID = atom2[i];
+        }
+        else {
+            centralID = atom2[i];
+            peripheralID = atom1[i];
+        }
+
+        // Add it to the cluster.
+
+        if (clusters.find(centralID) == clusters.end()) {
+            clusters[centralID] = ShakeCluster(centralID, 1.0/system.getParticleMass(centralID));
+        }
+        ShakeCluster& cluster = clusters[centralID];
+        cluster.addAtom(peripheralID, distance[i], 1.0/system.getParticleMass(peripheralID));
+        if (constraintCount[peripheralID] != 1 || invalidForShake[atom1[i]] || invalidForShake[atom2[i]]) {
+            cluster.markInvalid(clusters, invalidForShake);
+            map<int, ShakeCluster>::iterator otherCluster = clusters.find(peripheralID);
+            if (otherCluster != clusters.end() && otherCluster->second.valid)
+                otherCluster->second.markInvalid(clusters, invalidForShake);
+        }
+    }
+    int validShakeClusters = 0;
+    for (map<int, ShakeCluster>::iterator iter = clusters.begin(); iter != clusters.end(); ++iter) {
+        ShakeCluster& cluster = iter->second;
+        if (cluster.valid) {
+            cluster.valid = !invalidForShake[cluster.centralID];
+            for (int i = 0; i < cluster.size; i++)
+                if (invalidForShake[cluster.peripheralID[i]])
+                    cluster.valid = false;
+            if (cluster.valid)
+                ++validShakeClusters;
+        }
+    }
+
+    // Record the SHAKE clusters.
+
+    if (validShakeClusters > 0) {
+        vector<mm_int4> atoms(shakeAtoms->getSize());
+        vector<mm_float4> params(shakeParams->getSize());
+        int index = 0;
+        for (map<int, ShakeCluster>::const_iterator iter = clusters.begin(); iter != clusters.end(); ++iter) {
+            const ShakeCluster& cluster = iter->second;
+            if (!cluster.valid)
+                continue;
+            atoms.push_back((mm_int4) {cluster.centralID, cluster.peripheralID[0], (cluster.size > 1 ? cluster.peripheralID[1] : -1), (cluster.size > 2 ? cluster.peripheralID[2] : -1)});
+            params.push_back((mm_float4) {cluster.centralInvMass, 0.5f/(cluster.centralInvMass+cluster.peripheralInvMass), cluster.distance*cluster.distance, cluster.peripheralInvMass});
+            isShakeAtom[cluster.centralID] = true;
+            isShakeAtom[cluster.peripheralID[0]] = true;
+            if (cluster.size > 1)
+                isShakeAtom[cluster.peripheralID[1]] = true;
+            if (cluster.size > 2)
+                isShakeAtom[cluster.peripheralID[2]] = true;
+            ++index;
+        }
+        shakeAtoms = new OpenCLArray<mm_int4>(context, atoms.size(), "shakeAtoms");
+        shakeParams = new OpenCLArray<mm_float4>(context, params.size(), "shakeParams");
+        shakeAtoms->upload(atoms);
+        shakeParams->upload(params);
+    }
+}
+
+OpenCLIntegrationUtilities::~OpenCLIntegrationUtilities() {
+    if (posDelta != NULL)
+        delete posDelta;
+    if (oldPos != NULL)
+        delete oldPos;
+    if (settleAtoms != NULL)
+        delete settleAtoms;
+    if (settleParams != NULL)
+        delete settleParams;
+    if (shakeAtoms != NULL)
+        delete shakeAtoms;
+    if (shakeParams != NULL)
+        delete shakeParams;
+}
+
+void OpenCLIntegrationUtilities::applyConstraints(double tol, OpenCLArray<mm_float4>& oldPositions, OpenCLArray<mm_float4>& positionDeltas, OpenCLArray<mm_float4>& newDeltas) {
+    if (shakeAtoms != NULL) {
+        shakeKernel.setArg<cl_int>(0, shakeAtoms->getSize());
+        shakeKernel.setArg<cl_float>(1, tol);
+        shakeKernel.setArg<cl::Buffer>(2, oldPositions.getDeviceBuffer());
+        shakeKernel.setArg<cl::Buffer>(3, positionDeltas.getDeviceBuffer());
+        shakeKernel.setArg<cl::Buffer>(4, newDeltas.getDeviceBuffer());
+        shakeKernel.setArg<cl::Buffer>(5, shakeAtoms->getDeviceBuffer());
+        shakeKernel.setArg<cl::Buffer>(6, shakeParams->getDeviceBuffer());
+        context.executeKernel(shakeKernel, shakeAtoms->getSize());
+    }
+}
+//
+//void OpenCLIntegrationUtilities::clearBuffer(OpenCLArray<mm_float4>& array) {
+//    clearBufferKernel.setArg<cl::Buffer>(0, array.getDeviceBuffer());
+//    clearBufferKernel.setArg<cl_int>(1, array.getSize()*4);
+//    executeKernel(clearBufferKernel, array.getSize()*4);
+//}
+//
+//void OpenCLIntegrationUtilities::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/OpenCLIntegrationUtilities.h

+#ifndef OPENMM_OPENCLINTEGRATIONUTILITIES_H_
+#define OPENMM_OPENCLINTEGRATIONUTILITIES_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 "OpenMM/System.h"
+#include "OpenCLContext.h"
+
+namespace OpenMM {
+
+/**
+ * This class implements features that are used by many different integrators, including
+ * common workspace arrays and enforcing constraints.
+ */
+
+class OpenCLIntegrationUtilities {
+public:
+    OpenCLIntegrationUtilities(OpenCLContext& context, const System& system);
+    ~OpenCLIntegrationUtilities();
+    /**
+     * Get the array which contains position deltas.
+     */
+    OpenCLArray<mm_float4>& getPosDelta() {
+        return *posDelta;
+    }
+    /**
+     * Get the array which contains positions from before the current time step.
+     */
+    OpenCLArray<mm_float4>& getOldPos() {
+        return *oldPos;
+    }
+    /**
+     * Apply constraints to the atom positions.
+     *
+     * @param tol             the constraint tolerance
+     * @param oldPositions    the atom positions from before the current time step
+     * @param positionDeltas  the offsets being added to atom positions
+     * @param newDeltas       the array to store constrained deltas into
+     */
+    void applyConstraints(double tol, OpenCLArray<mm_float4>& oldPositions, OpenCLArray<mm_float4>& positionDeltas, OpenCLArray<mm_float4>& newDeltas);
+private:
+    OpenCLContext& context;
+//    cl::Kernel clearBufferKernel;
+    cl::Kernel shakeKernel;
+    OpenCLArray<mm_float4>* posDelta;
+    OpenCLArray<mm_float4>* oldPos;
+    OpenCLArray<mm_int4>* settleAtoms;
+    OpenCLArray<mm_float2>* settleParams;
+    OpenCLArray<mm_int4>* shakeAtoms;
+    OpenCLArray<mm_float4>* shakeParams;
+    struct ShakeCluster;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_OPENCLINTEGRATIONUTILITIES_H_*/

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -33,39 +33,40 @@ using namespace OpenMM;
 
 KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
     OpenCLPlatform::PlatformData& data = *static_cast<OpenCLPlatform::PlatformData*>(context.getPlatformData());
+    OpenCLContext& cl = *data.context;
     if (name == CalcForcesAndEnergyKernel::Name())
-        return new OpenCLCalcForcesAndEnergyKernel(name, platform, data);
+        return new OpenCLCalcForcesAndEnergyKernel(name, platform, cl);
     if (name == UpdateStateDataKernel::Name())
-        return new OpenCLUpdateStateDataKernel(name, platform, data);
+        return new OpenCLUpdateStateDataKernel(name, platform, cl);
     if (name == CalcHarmonicBondForceKernel::Name())
-        return new OpenCLCalcHarmonicBondForceKernel(name, platform, data, context.getSystem());
+        return new OpenCLCalcHarmonicBondForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcHarmonicAngleForceKernel::Name())
-        return new OpenCLCalcHarmonicAngleForceKernel(name, platform, data, context.getSystem());
+        return new OpenCLCalcHarmonicAngleForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcPeriodicTorsionForceKernel::Name())
-        return new OpenCLCalcPeriodicTorsionForceKernel(name, platform, data, context.getSystem());
+        return new OpenCLCalcPeriodicTorsionForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcRBTorsionForceKernel::Name())
-        return new OpenCLCalcRBTorsionForceKernel(name, platform, data, context.getSystem());
+        return new OpenCLCalcRBTorsionForceKernel(name, platform, cl, context.getSystem());
 //    if (name == CalcNonbondedForceKernel::Name())
-//        return new OpenCLCalcNonbondedForceKernel(name, platform, data, context.getSystem());
+//        return new OpenCLCalcNonbondedForceKernel(name, platform, cl, context.getSystem());
 //    if (name == CalcCustomNonbondedForceKernel::Name())
-//        return new OpenCLCalcCustomNonbondedForceKernel(name, platform, data, context.getSystem());
+//        return new OpenCLCalcCustomNonbondedForceKernel(name, platform, cl, context.getSystem());
 //    if (name == CalcGBSAOBCForceKernel::Name())
-//        return new OpenCLCalcGBSAOBCForceKernel(name, platform, data);
+//        return new OpenCLCalcGBSAOBCForceKernel(name, platform, cl);
     if (name == IntegrateVerletStepKernel::Name())
-        return new OpenCLIntegrateVerletStepKernel(name, platform, data);
+        return new OpenCLIntegrateVerletStepKernel(name, platform, cl);
 //    if (name == IntegrateLangevinStepKernel::Name())
-//        return new OpenCLIntegrateLangevinStepKernel(name, platform, data);
+//        return new OpenCLIntegrateLangevinStepKernel(name, platform, cl);
 //    if (name == IntegrateBrownianStepKernel::Name())
-//        return new OpenCLIntegrateBrownianStepKernel(name, platform, data);
+//        return new OpenCLIntegrateBrownianStepKernel(name, platform, cl);
 //    if (name == IntegrateVariableVerletStepKernel::Name())
-//        return new OpenCLIntegrateVariableVerletStepKernel(name, platform, data);
+//        return new OpenCLIntegrateVariableVerletStepKernel(name, platform, cl);
 //    if (name == IntegrateVariableLangevinStepKernel::Name())
-//        return new OpenCLIntegrateVariableLangevinStepKernel(name, platform, data);
+//        return new OpenCLIntegrateVariableLangevinStepKernel(name, platform, cl);
 //    if (name == ApplyAndersenThermostatKernel::Name())
-//        return new OpenCLApplyAndersenThermostatKernel(name, platform, data);
+//        return new OpenCLApplyAndersenThermostatKernel(name, platform, cl);
     if (name == CalcKineticEnergyKernel::Name())
-        return new OpenCLCalcKineticEnergyKernel(name, platform, data);
+        return new OpenCLCalcKineticEnergyKernel(name, platform, cl);
 //    if (name == RemoveCMMotionKernel::Name())
-//        return new OpenCLRemoveCMMotionKernel(name, platform, data);
+//        return new OpenCLRemoveCMMotionKernel(name, platform, cl);
     throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());
 }

platforms/opencl/src/OpenCLKernels.cpp

@@ -29,6 +29,7 @@
 #include "openmm/LangevinIntegrator.h"
 #include "openmm/Context.h"
 #include "openmm/internal/ContextImpl.h"
+#include "OpenCLIntegrationUtilities.h"
 #include <cmath>
 
 using namespace OpenMM;
@@ -38,19 +39,19 @@ void OpenCLCalcForcesAndEnergyKernel::initialize(const System& system) {
 }
 
 void OpenCLCalcForcesAndEnergyKernel::beginForceComputation(ContextImpl& context) {
-    data.context->clearBuffer(data.context->getForceBuffers());
+    cl.clearBuffer(cl.getForceBuffers());
 }
 
 void OpenCLCalcForcesAndEnergyKernel::finishForceComputation(ContextImpl& context) {
-    data.context->reduceBuffer(data.context->getForceBuffers(), data.context->getNumForceBuffers());
+    cl.reduceBuffer(cl.getForceBuffers(), cl.getNumForceBuffers());
 }
 
 void OpenCLCalcForcesAndEnergyKernel::beginEnergyComputation(ContextImpl& context) {
-    data.context->clearBuffer(data.context->getEnergyBuffer());
+    cl.clearBuffer(cl.getEnergyBuffer());
 }
 
 double OpenCLCalcForcesAndEnergyKernel::finishEnergyComputation(ContextImpl& context) {
-    OpenCLArray<cl_float>& energy = data.context->getEnergyBuffer();
+    OpenCLArray<cl_float>& energy = cl.getEnergyBuffer();
     energy.download();
     double sum = 0.0f;
     for (int i = 0; i < energy.getSize(); i++)
@@ -62,17 +63,17 @@ void OpenCLUpdateStateDataKernel::initialize(const System& system) {
 }
 
 double OpenCLUpdateStateDataKernel::getTime(const ContextImpl& context) const {
-    return data.time;
+    return cl.getTime();
 }
 
 void OpenCLUpdateStateDataKernel::setTime(ContextImpl& context, double time) {
-    data.time = time;
+    cl.setTime(time);
 }
 
 void OpenCLUpdateStateDataKernel::getPositions(ContextImpl& context, std::vector<Vec3>& positions) {
-    OpenCLArray<mm_float4>& posq = data.context->getPosq();
+    OpenCLArray<mm_float4>& posq = cl.getPosq();
     posq.download();
-    OpenCLArray<cl_int>& order = data.context->getAtomIndex();
+    OpenCLArray<cl_int>& order = cl.getAtomIndex();
     int numParticles = context.getSystem().getNumParticles();
     positions.resize(numParticles);
     for (int i = 0; i < numParticles; ++i) {
@@ -84,8 +85,8 @@ void OpenCLUpdateStateDataKernel::getPositions(ContextImpl& context, std::vector
 }
 
 void OpenCLUpdateStateDataKernel::setPositions(ContextImpl& context, const std::vector<Vec3>& positions) {
-    OpenCLArray<mm_float4>& posq = data.context->getPosq();
-    OpenCLArray<cl_int>& order = data.context->getAtomIndex();
+    OpenCLArray<mm_float4>& posq = cl.getPosq();
+    OpenCLArray<cl_int>& order = cl.getAtomIndex();
     int numParticles = context.getSystem().getNumParticles();
     for (int i = 0; i < numParticles; ++i) {
         mm_float4& pos = posq[i];
@@ -100,9 +101,9 @@ void OpenCLUpdateStateDataKernel::setPositions(ContextImpl& context, const std::
 }
 
 void OpenCLUpdateStateDataKernel::getVelocities(ContextImpl& context, std::vector<Vec3>& velocities) {
-    OpenCLArray<mm_float4>& velm = data.context->getVelm();
+    OpenCLArray<mm_float4>& velm = cl.getVelm();
     velm.download();
-    OpenCLArray<cl_int>& order = data.context->getAtomIndex();
+    OpenCLArray<cl_int>& order = cl.getAtomIndex();
     int numParticles = context.getSystem().getNumParticles();
     velocities.resize(numParticles);
     for (int i = 0; i < numParticles; ++i) {
@@ -112,8 +113,8 @@ void OpenCLUpdateStateDataKernel::getVelocities(ContextImpl& context, std::vecto
 }
 
 void OpenCLUpdateStateDataKernel::setVelocities(ContextImpl& context, const std::vector<Vec3>& velocities) {
-    OpenCLArray<mm_float4>& velm = data.context->getVelm();
-    OpenCLArray<cl_int>& order = data.context->getAtomIndex();
+    OpenCLArray<mm_float4>& velm = cl.getVelm();
+    OpenCLArray<cl_int>& order = cl.getAtomIndex();
     int numParticles = context.getSystem().getNumParticles();
     for (int i = 0; i < numParticles; ++i) {
         mm_float4& vel = velm[i];
@@ -126,9 +127,9 @@ void OpenCLUpdateStateDataKernel::setVelocities(ContextImpl& context, const std:
 }
 
 void OpenCLUpdateStateDataKernel::getForces(ContextImpl& context, std::vector<Vec3>& forces) {
-    OpenCLArray<mm_float4>& force = data.context->getForce();
+    OpenCLArray<mm_float4>& force = cl.getForce();
     force.download();
-    OpenCLArray<cl_int>& order = data.context->getAtomIndex();
+    OpenCLArray<cl_int>& order = cl.getAtomIndex();
     int numParticles = context.getSystem().getNumParticles();
     forces.resize(numParticles);
     for (int i = 0; i < numParticles; ++i) {
@@ -172,8 +173,8 @@ OpenCLCalcHarmonicBondForceKernel::~OpenCLCalcHarmonicBondForceKernel() {
 
 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");
+    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);
@@ -191,20 +192,20 @@ void OpenCLCalcHarmonicBondForceKernel::initialize(const System& system, const H
     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"));
+    cl.addForce(new OpenCLBondForceInfo(maxBuffers, force));
+    cl::Program program = cl.createProgram(cl.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>(0, cl.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>(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());
-    data.context->executeKernel(kernel, numBonds);
+    cl.executeKernel(kernel, numBonds);
 }
 
 double OpenCLCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
@@ -248,8 +249,8 @@ OpenCLCalcHarmonicAngleForceKernel::~OpenCLCalcHarmonicAngleForceKernel() {
 
 void OpenCLCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
     numAngles = force.getNumAngles();
-    params = new OpenCLArray<mm_float2>(*data.context, numAngles, "angleParams");
-    indices = new OpenCLArray<mm_int8>(*data.context, numAngles, "angleIndices");
+    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);
@@ -268,20 +269,20 @@ void OpenCLCalcHarmonicAngleForceKernel::initialize(const System& system, const
     for (int i = 0; i < forceBufferCounter.size(); i++) {
         maxBuffers = max(maxBuffers, forceBufferCounter[i]);
     }
-    data.context->addForce(new OpenCLAngleForceInfo(maxBuffers, force));
-    cl::Program program = data.context->createProgram(data.context->loadSourceFromFile("harmonicAngleForce.cl"));
+    cl.addForce(new OpenCLAngleForceInfo(maxBuffers, force));
+    cl::Program program = cl.createProgram(cl.loadSourceFromFile("harmonicAngleForce.cl"));
     kernel = cl::Kernel(program, "calcHarmonicAngleForce");
 }
 
 void OpenCLCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
-    kernel.setArg<cl_int>(0, data.context->getPaddedNumAtoms());
+    kernel.setArg<cl_int>(0, cl.getPaddedNumAtoms());
     kernel.setArg<cl_int>(1, numAngles);
-    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>(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());
-    data.context->executeKernel(kernel, numAngles);
+    cl.executeKernel(kernel, numAngles);
 }
 
 double OpenCLCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
@@ -326,8 +327,8 @@ OpenCLCalcPeriodicTorsionForceKernel::~OpenCLCalcPeriodicTorsionForceKernel() {
 
 void OpenCLCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
     numTorsions = force.getNumTorsions();
-    params = new OpenCLArray<mm_float4>(*data.context, numTorsions, "periodicTorsionParams");
-    indices = new OpenCLArray<mm_int8>(*data.context, numTorsions, "periodicTorsionIndices");
+    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);
@@ -346,20 +347,20 @@ void OpenCLCalcPeriodicTorsionForceKernel::initialize(const System& system, cons
     for (int i = 0; i < forceBufferCounter.size(); i++) {
         maxBuffers = max(maxBuffers, forceBufferCounter[i]);
     }
-    data.context->addForce(new OpenCLPeriodicTorsionForceInfo(maxBuffers, force));
-    cl::Program program = data.context->createProgram(data.context->loadSourceFromFile("periodicTorsionForce.cl"));
+    cl.addForce(new OpenCLPeriodicTorsionForceInfo(maxBuffers, force));
+    cl::Program program = cl.createProgram(cl.loadSourceFromFile("periodicTorsionForce.cl"));
     kernel = cl::Kernel(program, "calcPeriodicTorsionForce");
 }
 
 void OpenCLCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
-    kernel.setArg<cl_int>(0, data.context->getPaddedNumAtoms());
+    kernel.setArg<cl_int>(0, cl.getPaddedNumAtoms());
     kernel.setArg<cl_int>(1, numTorsions);
-    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>(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());
-    data.context->executeKernel(kernel, numTorsions);
+    cl.executeKernel(kernel, numTorsions);
 }
 
 double OpenCLCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context) {
@@ -404,8 +405,8 @@ OpenCLCalcRBTorsionForceKernel::~OpenCLCalcRBTorsionForceKernel() {
 
 void OpenCLCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
     numTorsions = force.getNumTorsions();
-    params = new OpenCLArray<mm_float8>(*data.context, numTorsions, "rbTorsionParams");
-    indices = new OpenCLArray<mm_int8>(*data.context, numTorsions, "rbTorsionIndices");
+    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);
@@ -424,20 +425,20 @@ void OpenCLCalcRBTorsionForceKernel::initialize(const System& system, const RBTo
     for (int i = 0; i < forceBufferCounter.size(); i++) {
         maxBuffers = max(maxBuffers, forceBufferCounter[i]);
     }
-    data.context->addForce(new OpenCLRBTorsionForceInfo(maxBuffers, force));
-    cl::Program program = data.context->createProgram(data.context->loadSourceFromFile("rbTorsionForce.cl"));
+    cl.addForce(new OpenCLRBTorsionForceInfo(maxBuffers, force));
+    cl::Program program = cl.createProgram(cl.loadSourceFromFile("rbTorsionForce.cl"));
     kernel = cl::Kernel(program, "calcRBTorsionForce");
 }
 
 void OpenCLCalcRBTorsionForceKernel::executeForces(ContextImpl& context) {
-    kernel.setArg<cl_int>(0, data.context->getPaddedNumAtoms());
+    kernel.setArg<cl_int>(0, cl.getPaddedNumAtoms());
     kernel.setArg<cl_int>(1, numTorsions);
-    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>(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());
-    data.context->executeKernel(kernel, numTorsions);
+    cl.executeKernel(kernel, numTorsions);
 }
 
 double OpenCLCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
@@ -794,31 +795,45 @@ OpenCLIntegrateVerletStepKernel::~OpenCLIntegrateVerletStepKernel() {
 }
 
 void OpenCLIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
-//    initializeIntegration(system, data, integrator);
-    data.context->initialize(system);
-    prevStepSize = -1.0;
+    cl.initialize(system);
+    cl::Program program = cl.createProgram(cl.loadSourceFromFile("verlet.cl"));
+    kernel1 = cl::Kernel(program, "integrateVerletPart1");
+    kernel2 = cl::Kernel(program, "integrateVerletPart2");
 }
 
 void OpenCLIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
-//    _gpuContext* gpu = data.gpu;
-//    double stepSize = integrator.getStepSize();
-//    if (stepSize != prevStepSize) {
-//        // Initialize the GPU parameters.
-//
-//        gpuSetVerletIntegrationParameters(gpu, (float) stepSize, 0.0f);
-//        gpuSetConstants(gpu);
-//        prevStepSize = stepSize;
-//    }
-//    kVerletUpdatePart1(gpu);
-//    kApplyFirstShake(gpu);
-//    kApplyFirstSettle(gpu);
-//    kApplyFirstCCMA(gpu);
-//    if (data.removeCM)
-//        if (data.stepCount%data.cmMotionFrequency == 0)
-//            gpu->bCalculateCM = true;
-//    kVerletUpdatePart2(gpu);
-//    data.time += stepSize;
-//    data.stepCount++;
+    OpenCLIntegrationUtilities& integration = cl.getIntegrationUtilties();
+    int numAtoms = cl.getNumAtoms();
+    double dt = integrator.getStepSize();
+
+    // Call the first integration kernel.
+
+    kernel1.setArg<cl_int>(0, numAtoms);
+    kernel1.setArg<cl_float>(1, dt);
+    kernel1.setArg<cl::Buffer>(2, cl.getPosq().getDeviceBuffer());
+    kernel1.setArg<cl::Buffer>(3, cl.getVelm().getDeviceBuffer());
+    kernel1.setArg<cl::Buffer>(4, cl.getForce().getDeviceBuffer());
+    kernel1.setArg<cl::Buffer>(5, integration.getPosDelta().getDeviceBuffer());
+    kernel1.setArg<cl::Buffer>(6, integration.getOldPos().getDeviceBuffer());
+    cl.executeKernel(kernel1, numAtoms);
+
+    // Apply constraints.
+
+    integration.applyConstraints(integrator.getConstraintTolerance(), integration.getOldPos(), integration.getPosDelta(), integration.getPosDelta());
+
+    // Call the second integration kernel.
+
+    kernel2.setArg<cl_int>(0, numAtoms);
+    kernel2.setArg<cl_float>(1, dt);
+    kernel2.setArg<cl::Buffer>(2, cl.getPosq().getDeviceBuffer());
+    kernel2.setArg<cl::Buffer>(3, cl.getVelm().getDeviceBuffer());
+    kernel2.setArg<cl::Buffer>(4, integration.getPosDelta().getDeviceBuffer());
+    cl.executeKernel(kernel2, numAtoms);
+
+    // Update the time and step count.
+
+    cl.setTime(cl.getTime()+dt);
+    cl.setStepCount(cl.getStepCount()+1);
 }
 
 //OpenCLIntegrateLangevinStepKernel::~OpenCLIntegrateLangevinStepKernel() {
@@ -1023,7 +1038,7 @@ double OpenCLCalcKineticEnergyKernel::execute(ContextImpl& context) {
     // We don't currently have a GPU kernel to do this, so we retrieve the velocities and calculate the energy
     // on the CPU.
 
-    OpenCLArray<mm_float4>& velm = data.context->getVelm();
+    OpenCLArray<mm_float4>& velm = cl.getVelm();
     double energy = 0.0;
     for (size_t i = 0; i < masses.size(); ++i) {
         mm_float4 v = velm[i];

platforms/opencl/src/OpenCLKernels.h

@@ -42,7 +42,7 @@ namespace OpenMM {
  */
 class OpenCLCalcForcesAndEnergyKernel : public CalcForcesAndEnergyKernel {
 public:
-    OpenCLCalcForcesAndEnergyKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data) : CalcForcesAndEnergyKernel(name, platform), data(data) {
+    OpenCLCalcForcesAndEnergyKernel(std::string name, const Platform& platform, OpenCLContext& cl) : CalcForcesAndEnergyKernel(name, platform), cl(cl) {
     }
     /**
      * Initialize the kernel.
@@ -82,7 +82,7 @@ public:
      */
     double finishEnergyComputation(ContextImpl& context);
 private:
-   OpenCLPlatform::PlatformData& data;
+   OpenCLContext& cl;
 };
 
 /**
@@ -91,7 +91,7 @@ private:
  */
 class OpenCLUpdateStateDataKernel : public UpdateStateDataKernel {
 public:
-    OpenCLUpdateStateDataKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data) : UpdateStateDataKernel(name, platform), data(data) {
+    OpenCLUpdateStateDataKernel(std::string name, const Platform& platform, OpenCLContext& cl) : UpdateStateDataKernel(name, platform), cl(cl) {
     }
     /**
      * Initialize the kernel.
@@ -142,7 +142,7 @@ public:
      */
     void getForces(ContextImpl& context, std::vector<Vec3>& forces);
 private:
-    OpenCLPlatform::PlatformData& data;
+    OpenCLContext& cl;
 };
 
 /**
@@ -150,8 +150,8 @@ private:
  */
 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(std::string name, const Platform& platform, OpenCLContext& cl, System& system) :
+        CalcHarmonicBondForceKernel(name, platform), cl(cl), system(system), params(NULL), indices(NULL) {
     }
     ~OpenCLCalcHarmonicBondForceKernel();
     /**
@@ -176,7 +176,7 @@ public:
     double executeEnergy(ContextImpl& context);
 private:
     int numBonds;
-    OpenCLPlatform::PlatformData& data;
+    OpenCLContext& cl;
     System& system;
     OpenCLArray<mm_float2>* params;
     OpenCLArray<mm_int4>* indices;
@@ -188,7 +188,7 @@ private:
  */
 class OpenCLCalcHarmonicAngleForceKernel : public CalcHarmonicAngleForceKernel {
 public:
-    OpenCLCalcHarmonicAngleForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcHarmonicAngleForceKernel(name, platform), data(data), system(system) {
+    OpenCLCalcHarmonicAngleForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcHarmonicAngleForceKernel(name, platform), cl(cl), system(system) {
     }
     ~OpenCLCalcHarmonicAngleForceKernel();
     /**
@@ -213,7 +213,7 @@ public:
     double executeEnergy(ContextImpl& context);
 private:
     int numAngles;
-    OpenCLPlatform::PlatformData& data;
+    OpenCLContext& cl;
     System& system;
     OpenCLArray<mm_float2>* params;
     OpenCLArray<mm_int8>* indices;
@@ -225,7 +225,7 @@ private:
  */
 class OpenCLCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
 public:
-    OpenCLCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcPeriodicTorsionForceKernel(name, platform), data(data), system(system) {
+    OpenCLCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcPeriodicTorsionForceKernel(name, platform), cl(cl), system(system) {
     }
     ~OpenCLCalcPeriodicTorsionForceKernel();
     /**
@@ -250,7 +250,7 @@ public:
     double executeEnergy(ContextImpl& context);
 private:
     int numTorsions;
-    OpenCLPlatform::PlatformData& data;
+    OpenCLContext& cl;
     System& system;
     OpenCLArray<mm_float4>* params;
     OpenCLArray<mm_int8>* indices;
@@ -262,7 +262,7 @@ private:
  */
 class OpenCLCalcRBTorsionForceKernel : public CalcRBTorsionForceKernel {
 public:
-    OpenCLCalcRBTorsionForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcRBTorsionForceKernel(name, platform), data(data), system(system) {
+    OpenCLCalcRBTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcRBTorsionForceKernel(name, platform), cl(cl), system(system) {
     }
     ~OpenCLCalcRBTorsionForceKernel();
     /**
@@ -287,7 +287,7 @@ public:
     double executeEnergy(ContextImpl& context);
 private:
     int numTorsions;
-    OpenCLPlatform::PlatformData& data;
+    OpenCLContext& cl;
     System& system;
     OpenCLArray<mm_float8>* params;
     OpenCLArray<mm_int8>* indices;
@@ -299,7 +299,7 @@ private:
 // */
 //class OpenCLCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
 //public:
-//    OpenCLCalcNonbondedForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcNonbondedForceKernel(name, platform), data(data), system(system) {
+//    OpenCLCalcNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcNonbondedForceKernel(name, platform), cl(cl), system(system) {
 //    }
 //    ~OpenCLCalcNonbondedForceKernel();
 //    /**
@@ -323,7 +323,7 @@ private:
 //     */
 //    double executeEnergy(ContextImpl& context);
 //private:
-//    OpenCLPlatform::PlatformData& data;
+//    OpenCLContext& cl;
 //    int numParticles;
 //    System& system;
 //};
@@ -333,7 +333,7 @@ private:
 // */
 //class OpenCLCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
 //public:
-//    OpenCLCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcCustomNonbondedForceKernel(name, platform), data(data), system(system) {
+//    OpenCLCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomNonbondedForceKernel(name, platform), cl(cl), system(system) {
 //    }
 //    ~OpenCLCalcCustomNonbondedForceKernel();
 //    /**
@@ -358,7 +358,7 @@ private:
 //    double executeEnergy(ContextImpl& context);
 //private:
 //    void updateGlobalParams(ContextImpl& context);
-//    OpenCLPlatform::PlatformData& data;
+//    OpenCLContext& cl;
 //    int numParticles;
 //    std::vector<std::string> globalParamNames;
 //    std::vector<float> globalParamValues;
@@ -370,7 +370,7 @@ private:
 // */
 //class OpenCLCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
 //public:
-//    OpenCLCalcGBSAOBCForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data) : CalcGBSAOBCForceKernel(name, platform), data(data) {
+//    OpenCLCalcGBSAOBCForceKernel(std::string name, const Platform& platform, OpenCLContext& cl) : CalcGBSAOBCForceKernel(name, platform), cl(cl) {
 //    }
 //    ~OpenCLCalcGBSAOBCForceKernel();
 //    /**
@@ -394,7 +394,7 @@ private:
 //     */
 //    double executeEnergy(ContextImpl& context);
 //private:
-//    OpenCLPlatform::PlatformData& data;
+//    OpenCLContext& cl;
 //};
 
 /**
@@ -402,7 +402,7 @@ private:
  */
 class OpenCLIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
 public:
-    OpenCLIntegrateVerletStepKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data) : IntegrateVerletStepKernel(name, platform), data(data) {
+    OpenCLIntegrateVerletStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateVerletStepKernel(name, platform), cl(cl) {
     }
     ~OpenCLIntegrateVerletStepKernel();
     /**
@@ -420,8 +420,9 @@ public:
      */
     void execute(ContextImpl& context, const VerletIntegrator& integrator);
 private:
-    OpenCLPlatform::PlatformData& data;
-    double prevStepSize;
+    OpenCLContext& cl;
+    cl::Kernel kernel1;
+    cl::Kernel kernel2;
 };
 
 ///**
@@ -429,7 +430,7 @@ private:
 // */
 //class OpenCLIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
 //public:
-//    OpenCLIntegrateLangevinStepKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data) : IntegrateLangevinStepKernel(name, platform), data(data) {
+//    OpenCLIntegrateLangevinStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateLangevinStepKernel(name, platform), cl(cl) {
 //    }
 //    ~OpenCLIntegrateLangevinStepKernel();
 //    /**
@@ -447,7 +448,7 @@ private:
 //     */
 //    void execute(ContextImpl& context, const LangevinIntegrator& integrator);
 //private:
-//    OpenCLPlatform::PlatformData& data;
+//    OpenCLContext& cl;
 //    double prevTemp, prevFriction, prevStepSize;
 //};
 //
@@ -456,7 +457,7 @@ private:
 // */
 //class OpenCLIntegrateBrownianStepKernel : public IntegrateBrownianStepKernel {
 //public:
-//    OpenCLIntegrateBrownianStepKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data) : IntegrateBrownianStepKernel(name, platform), data(data) {
+//    OpenCLIntegrateBrownianStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateBrownianStepKernel(name, platform), cl(cl) {
 //    }
 //    ~OpenCLIntegrateBrownianStepKernel();
 //    /**
@@ -474,7 +475,7 @@ private:
 //     */
 //    void execute(ContextImpl& context, const BrownianIntegrator& integrator);
 //private:
-//    OpenCLPlatform::PlatformData& data;
+//    OpenCLContext& cl;
 //    double prevTemp, prevFriction, prevStepSize;
 //};
 //
@@ -483,7 +484,7 @@ private:
 // */
 //class OpenCLIntegrateVariableVerletStepKernel : public IntegrateVariableVerletStepKernel {
 //public:
-//    OpenCLIntegrateVariableVerletStepKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data) : IntegrateVariableVerletStepKernel(name, platform), data(data) {
+//    OpenCLIntegrateVariableVerletStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateVariableVerletStepKernel(name, platform), cl(cl) {
 //    }
 //    ~OpenCLIntegrateVariableVerletStepKernel();
 //    /**
@@ -502,7 +503,7 @@ private:
 //     */
 //    void execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
 //private:
-//    OpenCLPlatform::PlatformData& data;
+//    OpenCLContext& cl;
 //    double prevErrorTol;
 //};
 //
@@ -511,7 +512,7 @@ private:
 // */
 //class OpenCLIntegrateVariableLangevinStepKernel : public IntegrateVariableLangevinStepKernel {
 //public:
-//    OpenCLIntegrateVariableLangevinStepKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data) : IntegrateVariableLangevinStepKernel(name, platform), data(data) {
+//    OpenCLIntegrateVariableLangevinStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateVariableLangevinStepKernel(name, platform), cl(cl) {
 //    }
 //    ~OpenCLIntegrateVariableLangevinStepKernel();
 //    /**
@@ -530,7 +531,7 @@ private:
 //     */
 //    void execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
 //private:
-//    OpenCLPlatform::PlatformData& data;
+//    OpenCLContext& cl;
 //    double prevTemp, prevFriction, prevErrorTol;
 //};
 //
@@ -539,7 +540,7 @@ private:
 // */
 //class OpenCLApplyAndersenThermostatKernel : public ApplyAndersenThermostatKernel {
 //public:
-//    OpenCLApplyAndersenThermostatKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data) : ApplyAndersenThermostatKernel(name, platform), data(data) {
+//    OpenCLApplyAndersenThermostatKernel(std::string name, const Platform& platform, OpenCLContext& cl) : ApplyAndersenThermostatKernel(name, platform), cl(cl) {
 //    }
 //    ~OpenCLApplyAndersenThermostatKernel();
 //    /**
@@ -556,7 +557,7 @@ private:
 //     */
 //    void execute(ContextImpl& context);
 //private:
-//    OpenCLPlatform::PlatformData& data;
+//    OpenCLContext& cl;
 //    double prevTemp, prevFrequency, prevStepSize;
 //};
 
@@ -565,7 +566,7 @@ private:
  */
 class OpenCLCalcKineticEnergyKernel : public CalcKineticEnergyKernel {
 public:
-    OpenCLCalcKineticEnergyKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data) : CalcKineticEnergyKernel(name, platform), data(data) {
+    OpenCLCalcKineticEnergyKernel(std::string name, const Platform& platform, OpenCLContext& cl) : CalcKineticEnergyKernel(name, platform), cl(cl) {
     }
     /**
      * Initialize the kernel.
@@ -580,7 +581,7 @@ public:
      */
     double execute(ContextImpl& context);
 private:
-    OpenCLPlatform::PlatformData& data;
+    OpenCLContext& cl;
     std::vector<double> masses;
 };
 
@@ -589,7 +590,7 @@ private:
 // */
 //class OpenCLRemoveCMMotionKernel : public RemoveCMMotionKernel {
 //public:
-//    OpenCLRemoveCMMotionKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data) : RemoveCMMotionKernel(name, platform), data(data) {
+//    OpenCLRemoveCMMotionKernel(std::string name, const Platform& platform, OpenCLContext& cl) : RemoveCMMotionKernel(name, platform), cl(cl) {
 //    }
 //    /**
 //     * Initialize the kernel, setting up the particle masses.
@@ -605,7 +606,7 @@ private:
 //     */
 //    void execute(ContextImpl& context);
 //private:
-//    OpenCLPlatform::PlatformData& data;
+//    OpenCLContext& cl;
 //};
 
 } // namespace OpenMM

platforms/opencl/src/kernels/shakeHydrogens.cl

+/**
+ * Enforce constraints on SHAKE clusters
+ */
+
+__kernel void applyShakeToHydrogens(int numClusters, float tol, __global float4* oldPos, __global float4* posDelta, __global float4* newDelta, __global int4* clusterAtoms, __global float4* clusterParams) {
+    int index = get_global_id(0);
+    while (index < numClusters) {
+        // Load the data for this cluster.
+
+        int4 atoms = clusterAtoms[index];
+        float4 params = clusterParams[index];
+        float4 pos = oldPos[atoms.x];
+        float4 xpi = posDelta[atoms.x];
+        float4 pos1 = oldPos[atoms.y];
+        float4 xpj1 = posDelta[atoms.y];
+        float4 pos2 = {0.0f, 0.0f, 0.0f, 0.0f};
+        float4 xpj2 = {0.0f, 0.0f, 0.0f, 0.0f};
+        float invMassCentral = params.x;
+        float avgMass = params.y;
+        float d2 = params.z;
+        float invMassPeripheral = params.w;
+        if (atoms.z != -1) {
+            pos2 = oldPos[atoms.z];
+            xpj2 = posDelta[atoms.z];
+        }
+        float4 pos3 = {0.0f, 0.0f, 0.0f, 0.0f};
+        float4 xpj3 = {0.0f, 0.0f, 0.0f, 0.0f};
+        if (atoms.w != -1) {
+            pos3 = oldPos[atoms.w];
+            xpj3 = posDelta[atoms.w];
+        }
+
+        // Precompute quantities.
+
+        float4 rij1 = pos-pos1;
+        float4 rij2 = pos-pos2;
+        float4 rij3 = pos-pos3;
+        float rij1sq = rij1.x*rij1.x + rij1.y*rij1.y + rij1.z*rij1.z;
+        float rij2sq = rij2.x*rij2.x + rij2.y*rij2.y + rij2.z*rij2.z;
+        float rij3sq = rij3.x*rij3.x + rij3.y*rij3.y + rij3.z*rij3.z;
+        float ld1 = d2-rij1sq;
+        float ld2 = d2-rij2sq;
+        float ld3 = d2-rij3sq;
+
+        // Iterate until convergence.
+
+        bool converged = false;
+        int iteration = 0;
+        while (iteration < 15 && !converged) {
+            converged = true;
+            float4 rpij = xpi-xpj1;
+            float rpsqij = rpij.x*rpij.x + rpij.y*rpij.y + rpij.z*rpij.z;
+            float rrpr = rij1.x*rpij.x + rij1.y*rpij.y + rij1.z*rpij.z;
+            float diff = fabs(ld1-2.0f*rrpr-rpsqij) / (d2*tol);
+            if (diff >= 1.0f) {
+                float acor  = (ld1-2.0f*rrpr-rpsqij)*avgMass / (rrpr+rij1sq);
+                float4 dr = rij1*acor;
+                xpi.xyz += dr.xyz*invMassCentral;
+                xpj1.xyz -= dr.xyz*invMassPeripheral;
+                converged = false;
+            }
+            if (atoms.z != -1) {
+                rpij.xyz = xpi.xyz-xpj2.xyz;
+                rpsqij = rpij.x*rpij.x + rpij.y*rpij.y + rpij.z*rpij.z;
+                rrpr = rij2.x*rpij.x + rij2.y*rpij.y + rij2.z*rpij.z;
+                diff = fabs(ld2-2.0f*rrpr-rpsqij) / (d2*tol);
+                if (diff >= 1.0f) {
+                    float acor  = (ld2 - 2.0f*rrpr - rpsqij)*avgMass / (rrpr + rij2sq);
+                    float4 dr = rij2*acor;
+                    xpi.xyz += dr.xyz*invMassCentral;
+                    xpj2.xyz -= dr.xyz*invMassPeripheral;
+                    converged = false;
+                }
+            }
+            if (atoms.w != -1) {
+                rpij.xyz = xpi.xyz-xpj3.xyz;
+                rpsqij = rpij.x*rpij.x + rpij.y*rpij.y + rpij.z*rpij.z;
+                rrpr = rij3.x*rpij.x + rij3.y*rpij.y + rij3.z*rpij.z;
+                diff = fabs(ld3 - 2.0f*rrpr - rpsqij) / (d2*tol);
+                if (diff >= 1.0f) {
+                    float acor  = (ld3-2.0f*rrpr-rpsqij)*avgMass / (rrpr+rij3sq);
+                    float4 dr = rij3*acor;
+                    xpi.xyz += dr.xyz*invMassCentral;
+                    xpj3.xyz -= dr.xyz*invMassPeripheral;
+                    converged = false;
+                }
+            }
+            iteration++;
+        }
+
+        // Record the new positions.
+
+        newDelta[atoms.x] = xpi;
+        newDelta[atoms.y] = xpj1;
+        if (atoms.z != -1)
+            newDelta[atoms.z] = xpj2;
+        if (atoms.w != -1)
+            newDelta[atoms.w] = xpj3;
+        index += get_global_size(0);
+    }
+}

platforms/opencl/src/kernels/verlet.cl

+/**
+ * Perform the first step of verlet integration.
+ */
+
+__kernel void integrateVerletPart1(int numAtoms, float dt, __global float4* posq, __global float4* velm, __global float4* force, __global float4* posDelta, __global float4* oldPos) {
+    int index = get_global_id(0);
+    while (index < numAtoms) {
+        float4 pos = posq[index];
+        float4 velocity = velm[index];
+        oldPos[index] = pos;
+        velocity.xyz += force[index].xyz*dt*velocity.w;
+        pos.xyz = velocity.xyz*dt;
+        posDelta[index] = pos;
+        velm[index] = velocity;
+        index += get_global_size(0);
+    }
+}
+
+/**
+ * Perform the second step of verlet integration.
+ */
+
+__kernel void integrateVerletPart2(int numAtoms, float dt, __global float4* posq, __global float4* velm, __global float4* posDelta) {
+    int index = get_global_id(0);
+    float oneOverDt = 1.0f/dt;
+    while (index < numAtoms) {
+        float4 pos = posq[index];
+        float4 delta = posDelta[index];
+        float4 velocity = velm[index];
+        pos.xyz += delta.xyz;
+        velocity.xyz = delta.xyz*oneOverDt;
+        posq[index] = pos;
+        velm[index] = velocity;
+        index += get_global_size(0);
+    }
+}

platforms/opencl/tests/TestOpenCLVerletIntegrator.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 VerletIntegrator.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "OpenCLPlatform.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "../src/sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testSingleBond() {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(2.0);
+    system.addParticle(2.0);
+    VerletIntegrator integrator(0.01);
+    HarmonicBondForce* forceField = new HarmonicBondForce();
+    forceField->addBond(0, 1, 1.5, 1);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    context.setPositions(positions);
+
+    // This is simply a harmonic oscillator, so compare it to the analytical solution.
+
+    const double freq = 1.0;;
+    State state = context.getState(State::Energy);
+    const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
+    for (int i = 0; i < 1000; ++i) {
+        state = context.getState(State::Positions | State::Velocities | State::Energy);
+        double time = state.getTime();
+        double expectedDist = 1.5+0.5*std::cos(freq*time);
+        ASSERT_EQUAL_VEC(Vec3(-0.5*expectedDist, 0, 0), state.getPositions()[0], 0.02);
+        ASSERT_EQUAL_VEC(Vec3(0.5*expectedDist, 0, 0), state.getPositions()[1], 0.02);
+        double expectedSpeed = -0.5*freq*std::sin(freq*time);
+        ASSERT_EQUAL_VEC(Vec3(-0.5*expectedSpeed, 0, 0), state.getVelocities()[0], 0.02);
+        ASSERT_EQUAL_VEC(Vec3(0.5*expectedSpeed, 0, 0), state.getVelocities()[1], 0.02);
+        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
+        integrator.step(1);
+    }
+}
+
+void testConstraints() {
+    const int numParticles = 8;
+    const int numConstraints = 5;
+    const double temp = 100.0;
+    OpenCLPlatform platform;
+    System system;
+    VerletIntegrator integrator(0.001);
+    integrator.setConstraintTolerance(1e-5);
+    NonbondedForce* forceField = new NonbondedForce();
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(10.0);
+        forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
+    }
+    system.addConstraint(0, 1, 1.0);
+    system.addConstraint(1, 2, 1.0);
+    system.addConstraint(2, 3, 1.0);
+    system.addConstraint(4, 5, 1.0);
+    system.addConstraint(6, 7, 1.0);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    init_gen_rand(0);
+    for (int i = 0; i < numParticles; ++i) {
+        positions[i] = Vec3(i/2, (i+1)/2, 0);
+        velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
+    }
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+
+    // Simulate it and see whether the constraints remain satisfied.
+
+    double initialEnergy = 0.0;
+    for (int i = 0; i < 1000; ++i) {
+        State state = context.getState(State::Positions | State::Energy);
+        for (int j = 0; j < numConstraints; ++j) {
+            int particle1, particle2;
+            double distance;
+            system.getConstraintParameters(j, particle1, particle2, distance);
+            Vec3 p1 = state.getPositions()[particle1];
+            Vec3 p2 = state.getPositions()[particle2];
+            double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
+            ASSERT_EQUAL_TOL(distance, dist, 1e-4);
+        }
+        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        if (i == 1)
+            initialEnergy = energy;
+        else if (i > 1)
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
+        integrator.step(1);
+    }
+}
+
+void testConstrainedClusters() {
+    const int numParticles = 7;
+    const double temp = 500.0;
+    OpenCLPlatform platform;
+    System system;
+    VerletIntegrator integrator(0.002);
+    integrator.setConstraintTolerance(1e-5);
+    NonbondedForce* forceField = new NonbondedForce();
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(i > 1 ? 1.0 : 10.0);
+        forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
+    }
+    system.addConstraint(0, 1, 1.0);
+    system.addConstraint(0, 2, 1.0);
+    system.addConstraint(0, 3, 1.0);
+    system.addConstraint(0, 4, 1.0);
+    system.addConstraint(1, 5, 1.0);
+    system.addConstraint(1, 6, 1.0);
+    system.addConstraint(2, 3, sqrt(2.0));
+    system.addConstraint(2, 4, sqrt(2.0));
+    system.addConstraint(3, 4, sqrt(2.0));
+    system.addConstraint(5, 6, sqrt(2.0));
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(numParticles);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    positions[2] = Vec3(-1, 0, 0);
+    positions[3] = Vec3(0, 1, 0);
+    positions[4] = Vec3(0, 0, 1);
+    positions[5] = Vec3(2, 0, 0);
+    positions[6] = Vec3(1, 1, 0);
+    vector<Vec3> velocities(numParticles);
+    init_gen_rand(0);
+    for (int i = 0; i < numParticles; ++i)
+        velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+
+    // Simulate it and see whether the constraints remain satisfied.
+
+    double initialEnergy = 0.0;
+    for (int i = 0; i < 1000; ++i) {
+        State state = context.getState(State::Positions | State::Energy);
+        for (int j = 0; j < system.getNumConstraints(); ++j) {
+            int particle1, particle2;
+            double distance;
+            system.getConstraintParameters(j, particle1, particle2, distance);
+            Vec3 p1 = state.getPositions()[particle1];
+            Vec3 p2 = state.getPositions()[particle2];
+            double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
+            ASSERT_EQUAL_TOL(distance, dist, 2e-5);
+        }
+        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        if (i == 1)
+            initialEnergy = energy;
+        else if (i > 1)
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
+        integrator.step(1);
+    }
+}
+
+int main() {
+    try {
+        testSingleBond();
+//        testConstraints();
+//        testConstrainedClusters();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+