Bug fixes and an optimization to CustomNonbondedForce with interaction groups

platforms/cuda/src/CudaKernels.cpp

@@ -46,6 +46,7 @@
 #include "lepton/ParsedExpression.h"
 #include "SimTKOpenMMRealType.h"
 #include "SimTKOpenMMUtilities.h"
+#include <algorithm>
 #include <cmath>
 #include <set>
 
@@ -2290,6 +2291,8 @@ double CudaCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool in
             interactionGroupArgs.push_back(&cu.getEnergyBuffer().getDevicePointer());
             interactionGroupArgs.push_back(&cu.getPosq().getDevicePointer());
             interactionGroupArgs.push_back(&interactionGroupData->getDevicePointer());
+            interactionGroupArgs.push_back(cu.getPeriodicBoxSizePointer());
+            interactionGroupArgs.push_back(cu.getInvPeriodicBoxSizePointer());
             for (int i = 0; i < (int) params->getBuffers().size(); i++)
                 interactionGroupArgs.push_back(&params->getBuffers()[i].getMemory());
             if (globals != NULL)

platforms/cuda/src/kernels/customNonbondedGroups.cu

@@ -11,7 +11,8 @@ typedef struct {
 } AtomData;
 
 extern "C" __global__ void computeInteractionGroups(
-        unsigned long long* __restrict__ forceBuffers, real* __restrict__ energyBuffer, const real4* __restrict__ posq, const int4* __restrict__ groupData
+        unsigned long long* __restrict__ forceBuffers, real* __restrict__ energyBuffer, const real4* __restrict__ posq, const int4* __restrict__ groupData,
+        real4 periodicBoxSize, real4 invPeriodicBoxSize
         PARAMETER_ARGUMENTS) {
     const unsigned int totalWarps = (blockDim.x*gridDim.x)/TILE_SIZE;
     const unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/TILE_SIZE; // global warpIndex
@@ -53,6 +54,9 @@ extern "C" __global__ void computeInteractionGroups(
             delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
             real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
+#ifdef USE_CUTOFF
+            if (!isExcluded && r2 < CUTOFF_SQUARED) {
+#endif
                 real invR = RSQRT(r2);
                 real r = RECIP(invR);
                 LOAD_ATOM2_PARAMETERS
@@ -67,6 +71,9 @@ extern "C" __global__ void computeInteractionGroups(
                 localData[localIndex].fx += delta.x;
                 localData[localIndex].fy += delta.y;
                 localData[localIndex].fz += delta.z;
+#ifdef USE_CUTOFF
+            }
+#endif
             tj = (tj == rangeEnd-1 ? rangeStart : tj+1);
         }
         if (exclusions != 0) {