ApplyConstraintsKernel did not work correctly

platforms/opencl/src/OpenCLKernels.cpp

@@ -221,8 +221,20 @@ void OpenCLApplyConstraintsKernel::initialize(const System& system) {
 }
 
 void OpenCLApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
-    cl.getIntegrationUtilities().applyConstraints(tol);
-    cl.getIntegrationUtilities().computeVirtualSites();
+    if (!hasInitializedKernel) {
+        hasInitializedKernel = true;
+        map<string, string> defines;
+        defines["NUM_ATOMS"] = intToString(cl.getNumAtoms());
+        cl::Program program = cl.createProgram(OpenCLKernelSources::constraints, defines);
+        applyDeltasKernel = cl::Kernel(program, "applyPositionDeltas");
+        applyDeltasKernel.setArg<cl::Buffer>(0, cl.getPosq().getDeviceBuffer());
+        applyDeltasKernel.setArg<cl::Buffer>(1, cl.getIntegrationUtilities().getPosDelta().getDeviceBuffer());
+    }
+    OpenCLIntegrationUtilities& integration = cl.getIntegrationUtilities();
+    cl.clearBuffer(integration.getPosDelta());
+    integration.applyConstraints(tol);
+    cl.executeKernel(applyDeltasKernel, cl.getNumAtoms());
+    integration.computeVirtualSites();
 }
 
 void OpenCLVirtualSitesKernel::initialize(const System& system) {

platforms/opencl/src/OpenCLKernels.h

@@ -159,7 +159,8 @@ private:
  */
 class OpenCLApplyConstraintsKernel : public ApplyConstraintsKernel {
 public:
-    OpenCLApplyConstraintsKernel(std::string name, const Platform& platform, OpenCLContext& cl) : ApplyConstraintsKernel(name, platform), cl(cl) {
+    OpenCLApplyConstraintsKernel(std::string name, const Platform& platform, OpenCLContext& cl) : ApplyConstraintsKernel(name, platform),
+            cl(cl), hasInitializedKernel(false) {
     }
     /**
      * Initialize the kernel.
@@ -176,6 +177,8 @@ public:
     void apply(ContextImpl& context, double tol);
 private:
     OpenCLContext& cl;
+    bool hasInitializedKernel;
+    cl::Kernel applyDeltasKernel;
 };
 
 /**

platforms/opencl/src/kernels/constraints.cl

+__kernel void applyPositionDeltas(__global float4* restrict posq, __global float4* restrict posDelta) {
+    for (unsigned int index = get_global_id(0); index < NUM_ATOMS; index += get_global_size(0)) {
+        float4 position = posq[index];
+        position.xyz += posDelta[index].xyz;
+        posq[index] = position;
+    }
+}