Fixed race condition

platforms/cuda/include/CudaKernels.h

@@ -709,7 +709,7 @@ private:
     CudaFFT3D* dispersionFft;
     cufftHandle dispersionFftForward;
     cufftHandle dispersionFftBackward;
-    CUfunction computeParamsKernel;
+    CUfunction computeParamsKernel, computeExclusionParamsKernel;
     CUfunction ewaldSumsKernel;
     CUfunction ewaldForcesKernel;
     CUfunction pmeGridIndexKernel;

platforms/cuda/src/CudaKernels.cpp

@@ -2110,6 +2110,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
     
     CUmodule module = cu.createModule(CudaKernelSources::nonbondedParameters, paramsDefines);
     computeParamsKernel = cu.getKernel(module, "computeParameters");
+    computeExclusionParamsKernel = cu.getKernel(module, "computeExclusionParameters");
     info = new ForceInfo(force);
     cu.addForce(info);
 }
@@ -2136,7 +2137,7 @@ double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeF
         vector<void*> paramsArgs = {&cu.getEnergyBuffer().getDevicePointer(), &computeSelfEnergy, &globalParams.getDevicePointer(), &numAtoms,
                 &baseParticleParams.getDevicePointer(), &cu.getPosq().getDevicePointer(), &charges.getDevicePointer(), &sigmaEpsilon.getDevicePointer(),
                 &particleParamOffsets.getDevicePointer(), &particleOffsetIndices.getDevicePointer()};
-        int numExceptions, numExclusions;
+        int numExceptions;
         if (exceptionParams.isInitialized()) {
             numExceptions = exceptionParams.getSize();
             paramsArgs.push_back(&numExceptions);
@@ -2145,13 +2146,13 @@ double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeF
             paramsArgs.push_back(&exceptionParamOffsets.getDevicePointer());
             paramsArgs.push_back(&exceptionOffsetIndices.getDevicePointer());
         }
+        cu.executeKernel(computeParamsKernel, &paramsArgs[0], cu.getPaddedNumAtoms());
         if (exclusionParams.isInitialized()) {
-            numExclusions = exclusionParams.getSize();
-            paramsArgs.push_back(&numExclusions);
-            paramsArgs.push_back(&exclusionAtoms.getDevicePointer());
-            paramsArgs.push_back(&exclusionParams.getDevicePointer());
+            int numExclusions = exclusionParams.getSize();
+            vector<void*> exclusionParamsArgs = {&cu.getPosq().getDevicePointer(), &charges.getDevicePointer(), &sigmaEpsilon.getDevicePointer(),
+                    &numExclusions, &exclusionAtoms.getDevicePointer(), &exclusionParams.getDevicePointer()};
+            cu.executeKernel(computeExclusionParamsKernel, &exclusionParamsArgs[0], numExclusions);
         }
-        cu.executeKernel(computeParamsKernel, &paramsArgs[0], cu.getPaddedNumAtoms());
         if (usePmeStream) {
             cuEventRecord(paramsSyncEvent, cu.getCurrentStream());
             cuStreamWaitEvent(pmeStream, paramsSyncEvent, 0);

platforms/cuda/src/kernels/nonbondedParameters.cu

@@ -7,9 +7,6 @@ extern "C" __global__ void computeParameters(mixed* __restrict__ energyBuffer, b
 #ifdef HAS_EXCEPTIONS
         , int numExceptions, const float4* __restrict__ baseExceptionParams, float4* __restrict__ exceptionParams,
         float4* __restrict__ exceptionParamOffsets, int* __restrict__ exceptionOffsetIndices
-#endif
-#ifdef HAS_EXCLUSIONS
-        , int numExclusions, const int2* __restrict__ exclusionAtoms, float4* __restrict__ exclusionParams
 #endif
         ) {
     mixed energy = 0;
@@ -63,10 +60,15 @@ extern "C" __global__ void computeParameters(mixed* __restrict__ energyBuffer, b
         exceptionParams[i] = make_float4((float) (138.935456f*params.x), (float) params.y, (float) (4*params.z), 0);
     }
 #endif
+    if (includeSelfEnergy)
+        energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += energy;
+}
 
-    // Compute parameters for subtracting the reciprocal part of excluded interactions.
-
-#ifdef HAS_EXCLUSIONS
+/**
+ * Compute parameters for subtracting the reciprocal part of excluded interactions.
+ */
+extern "C" __global__ void computeExclusionParameters(real4* __restrict__ posq, real* __restrict__ charge, float2* __restrict__ sigmaEpsilon,
+        int numExclusions, const int2* __restrict__ exclusionAtoms, float4* __restrict__ exclusionParams) {
     for (int i = blockIdx.x*blockDim.x+threadIdx.x; i < numExclusions; i += blockDim.x*gridDim.x) {
         int2 atoms = exclusionAtoms[i];
 #ifdef USE_POSQ_CHARGES
@@ -85,7 +87,4 @@ extern "C" __global__ void computeParameters(mixed* __restrict__ energyBuffer, b
 #endif
         exclusionParams[i] = make_float4((float) (138.935456f*chargeProd), sigma, epsilon, 0);
     }
-#endif
-    if (includeSelfEnergy)
-        energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += energy;
-}
+}
\ No newline at end of file

platforms/opencl/include/OpenCLKernels.h

@@ -685,7 +685,7 @@ private:
     Kernel cpuPme;
     PmeIO* pmeio;
     SyncQueuePostComputation* syncQueue;
-    cl::Kernel computeParamsKernel;
+    cl::Kernel computeParamsKernel, computeExclusionParamsKernel;
     cl::Kernel ewaldSumsKernel;
     cl::Kernel ewaldForcesKernel;
     cl::Kernel pmeAtomRangeKernel;

platforms/opencl/src/OpenCLKernels.cpp

@@ -2036,6 +2036,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
     
     cl::Program program = cl.createProgram(OpenCLKernelSources::nonbondedParameters, paramsDefines);
     computeParamsKernel = cl::Kernel(program, "computeParameters");
+    computeExclusionParamsKernel = cl::Kernel(program, "computeExclusionParameters");
     info = new ForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force);
     cl.addForce(info);
 }
@@ -2063,9 +2064,12 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
             computeParamsKernel.setArg<cl::Buffer>(index++, exceptionOffsetIndices.getDeviceBuffer());
         }
         if (exclusionParams.isInitialized()) {
-            computeParamsKernel.setArg<cl_int>(index++, exclusionParams.getSize());
-            computeParamsKernel.setArg<cl::Buffer>(index++, exclusionAtoms.getDeviceBuffer());
-            computeParamsKernel.setArg<cl::Buffer>(index++, exclusionParams.getDeviceBuffer());
+            computeExclusionParamsKernel.setArg<cl::Buffer>(0, cl.getPosq().getDeviceBuffer());
+            computeExclusionParamsKernel.setArg<cl::Buffer>(1, charges.getDeviceBuffer());
+            computeExclusionParamsKernel.setArg<cl::Buffer>(2, sigmaEpsilon.getDeviceBuffer());
+            computeExclusionParamsKernel.setArg<cl_int>(3, exclusionParams.getSize());
+            computeExclusionParamsKernel.setArg<cl::Buffer>(4, exclusionAtoms.getDeviceBuffer());
+            computeExclusionParamsKernel.setArg<cl::Buffer>(5, exclusionParams.getDeviceBuffer());
         }
         if (cosSinSums.isInitialized()) {
             ewaldSumsKernel.setArg<cl::Buffer>(0, cl.getEnergyBuffer().getDeviceBuffer());
@@ -2215,6 +2219,8 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
     if (recomputeParams || hasOffsets) {
         computeParamsKernel.setArg<cl_int>(1, includeEnergy && includeReciprocal);
         cl.executeKernel(computeParamsKernel, cl.getPaddedNumAtoms());
+        if (exclusionParams.isInitialized())
+            cl.executeKernel(computeExclusionParamsKernel, exclusionParams.getSize());
         if (usePmeQueue) {
             vector<cl::Event> events(1);
             cl.getQueue().enqueueMarker(&events[0]);

platforms/opencl/src/kernels/nonbondedParameters.cl

@@ -7,9 +7,6 @@ __kernel void computeParameters(__global mixed* restrict energyBuffer, int inclu
 #ifdef HAS_EXCEPTIONS
         , int numExceptions, __global const float4* restrict baseExceptionParams, __global float4* restrict exceptionParams,
         __global float4* restrict exceptionParamOffsets, __global int* restrict exceptionOffsetIndices
-#endif
-#ifdef HAS_EXCLUSIONS
-        , int numExclusions, __global const int2* restrict exclusionAtoms, __global float4* restrict exclusionParams
 #endif
         ) {
     mixed energy = 0;
@@ -63,10 +60,15 @@ __kernel void computeParameters(__global mixed* restrict energyBuffer, int inclu
         exceptionParams[i] = (float4) ((float) (138.935456f*params.x), (float) params.y, (float) (4*params.z), 0);
     }
 #endif
+    if (includeSelfEnergy)
+        energyBuffer[get_global_id(0)] += energy;
+}
 
-    // Compute parameters for subtracting the reciprocal part of excluded interactions.
-
-#ifdef HAS_EXCLUSIONS
+/**
+ * Compute parameters for subtracting the reciprocal part of excluded interactions.
+ */
+__kernel void computeExclusionParameters(__global real4* restrict posq, __global real* restrict charge, __global float2* restrict sigmaEpsilon,
+        int numExclusions, __global const int2* restrict exclusionAtoms, __global float4* restrict exclusionParams) {
     for (int i = get_global_id(0); i < numExclusions; i += get_global_size(0)) {
         int2 atoms = exclusionAtoms[i];
 #ifdef USE_POSQ_CHARGES
@@ -85,7 +87,4 @@ __kernel void computeParameters(__global mixed* restrict energyBuffer, int inclu
 #endif
         exclusionParams[i] = (float4) ((float) (138.935456f*chargeProd), sigma, epsilon, 0);
     }
-#endif
-    if (includeSelfEnergy)
-        energyBuffer[get_global_id(0)] += energy;
 }