Merge pull request #1127 from peastman/centroidbugs

Fixed various bugs in CustomCentroidBondForce

platforms/cuda/src/CudaKernels.cpp

@@ -4634,6 +4634,8 @@ void CudaCalcCustomCentroidBondForceKernel::initialize(const System& system, con
 }
 
 double CudaCalcCustomCentroidBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    if (numBonds == 0)
+        return 0.0;
     if (globals != NULL) {
         bool changed = false;
         for (int i = 0; i < (int) globalParamNames.size(); i++) {
@@ -4658,10 +4660,7 @@ double CudaCalcCustomCentroidBondForceKernel::execute(ContextImpl& context, bool
 
 void CudaCalcCustomCentroidBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force) {
     cu.setAsCurrent();
-    int numContexts = cu.getPlatformData().contexts.size();
-    int startIndex = cu.getContextIndex()*force.getNumBonds()/numContexts;
-    int endIndex = (cu.getContextIndex()+1)*force.getNumBonds()/numContexts;
-    if (numBonds != endIndex-startIndex)
+    if (numBonds != force.getNumBonds())
         throw OpenMMException("updateParametersInContext: The number of bonds has changed");
     if (numBonds == 0)
         return;
@@ -4672,7 +4671,7 @@ void CudaCalcCustomCentroidBondForceKernel::copyParametersToContext(ContextImpl&
     vector<int> particles;
     vector<double> parameters;
     for (int i = 0; i < numBonds; i++) {
-        force.getBondParameters(startIndex+i, particles, parameters);
+        force.getBondParameters(i, particles, parameters);
         paramVector[i].resize(parameters.size());
         for (int j = 0; j < (int) parameters.size(); j++)
             paramVector[i][j] = (float) parameters[j];

platforms/opencl/src/OpenCLKernels.cpp

@@ -4837,6 +4837,8 @@ void OpenCLCalcCustomCentroidBondForceKernel::initialize(const System& system, c
 }
 
 double OpenCLCalcCustomCentroidBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    if (numBonds == 0)
+        return 0.0;
     if (globals != NULL) {
         bool changed = false;
         for (int i = 0; i < (int) globalParamNames.size(); i++) {
@@ -4857,10 +4859,7 @@ double OpenCLCalcCustomCentroidBondForceKernel::execute(ContextImpl& context, bo
 }
 
 void OpenCLCalcCustomCentroidBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force) {
-    int numContexts = cl.getPlatformData().contexts.size();
-    int startIndex = cl.getContextIndex()*force.getNumBonds()/numContexts;
-    int endIndex = (cl.getContextIndex()+1)*force.getNumBonds()/numContexts;
-    if (numBonds != endIndex-startIndex)
+    if (numBonds != force.getNumBonds())
         throw OpenMMException("updateParametersInContext: The number of bonds has changed");
     if (numBonds == 0)
         return;
@@ -4871,7 +4870,7 @@ void OpenCLCalcCustomCentroidBondForceKernel::copyParametersToContext(ContextImp
     vector<int> particles;
     vector<double> parameters;
     for (int i = 0; i < numBonds; i++) {
-        force.getBondParameters(startIndex+i, particles, parameters);
+        force.getBondParameters(i, particles, parameters);
         paramVector[i].resize(parameters.size());
         for (int j = 0; j < (int) parameters.size(); j++)
             paramVector[i][j] = (float) parameters[j];

platforms/opencl/src/kernels/customCentroidBond.cl

@@ -27,30 +27,35 @@ __kernel void computeGroupCenters(__global const real4* restrict posq, __global
         temp[thread].x = center.x;
         temp[thread].y = center.y;
         temp[thread].z = center.z;
-        __syncthreads();
+        barrier(CLK_LOCAL_MEM_FENCE);
         if (thread < 32) {
             temp[thread].x += temp[thread+32].x;
             temp[thread].y += temp[thread+32].y;
             temp[thread].z += temp[thread+32].z;
+            SYNC_WARPS;
             if (thread < 16) {
                 temp[thread].x += temp[thread+16].x;
                 temp[thread].y += temp[thread+16].y;
                 temp[thread].z += temp[thread+16].z;
+                SYNC_WARPS;
             }
             if (thread < 8) {
                 temp[thread].x += temp[thread+8].x;
                 temp[thread].y += temp[thread+8].y;
                 temp[thread].z += temp[thread+8].z;
+                SYNC_WARPS;
             }
             if (thread < 4) {
                 temp[thread].x += temp[thread+4].x;
                 temp[thread].y += temp[thread+4].y;
                 temp[thread].z += temp[thread+4].z;
+                SYNC_WARPS;
             }
             if (thread < 2) {
                 temp[thread].x += temp[thread+2].x;
                 temp[thread].y += temp[thread+2].y;
                 temp[thread].z += temp[thread+2].z;
+                SYNC_WARPS;
             }
         }
         if (thread == 0)