Implemented cutoffs and periodic boundary conditions

platforms/opencl/src/OpenCLKernels.cpp

@@ -530,6 +530,16 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
     sigmaEpsilon->upload(sigmaEpsilonVector);
     bool useCutoff = (force.getNonbondedMethod() != NonbondedForce::NoCutoff);
     bool usePeriodic = (force.getNonbondedMethod() != NonbondedForce::NoCutoff && force.getNonbondedMethod() != NonbondedForce::CutoffNonPeriodic);
+    map<string, string> defines;
+    if (useCutoff) {
+        double reactionFieldK = pow(force.getCutoffDistance(), -3.0)*(force.getReactionFieldDielectric()-1.0)/(2.0*force.getReactionFieldDielectric()+1.0);
+        double reactionFieldC = (1.0 / force.getCutoffDistance())*(3.0*force.getReactionFieldDielectric())/(2.0*force.getReactionFieldDielectric()+1.0);
+        char k[50], c[50];
+        sprintf(k, "%.8ef", reactionFieldK);
+        sprintf(c, "%.8ef", reactionFieldC);
+        defines["REACTION_FIELD_K"] = string(k);
+        defines["REACTION_FIELD_C"] = string(c);
+    }
 //    if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
 //        method = CUTOFF;
 //    }
@@ -566,7 +576,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
 //    }
 //    data.nonbondedMethod = method;
 //    gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList, method);
-    cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, force.getCutoffDistance(), exclusionList);
+    cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, force.getCutoffDistance(), exclusionList, defines);
     cl.getNonbondedUtilities().addParameter("sigmaEpsilon", "float2", 8, sigmaEpsilon->getDeviceBuffer());
     cutoffSquared = force.getCutoffDistance()*force.getCutoffDistance();
 
@@ -602,7 +612,12 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
             maxBuffers = max(maxBuffers, forceBufferCounter[i]);
     }
     cl.addForce(new OpenCLNonbondedForceInfo(maxBuffers, force));
-    cl::Program program = cl.createProgram(cl.loadSourceFromFile("nonbondedExceptions.cl"));
+    if (useCutoff) {
+        defines["USE_CUTOFF"] = "1";
+    }
+    if (usePeriodic)
+        defines["USE_PERIODIC"] = "1";
+    cl::Program program = cl.createProgram(cl.loadSourceFromFile("nonbondedExceptions.cl"), defines);
     exceptionsKernel = cl::Kernel(program, "computeNonbondedExceptions");
 }
 

platforms/opencl/src/OpenCLNonbondedUtilities.cpp

@@ -68,7 +68,7 @@ OpenCLNonbondedUtilities::~OpenCLNonbondedUtilities() {
         delete compact;
 }
 
-void OpenCLNonbondedUtilities::addInteraction(bool usesCutoff, bool usesPeriodic, double cutoffDistance, const std::vector<std::vector<int> >& exclusionList) {
+void OpenCLNonbondedUtilities::addInteraction(bool usesCutoff, bool usesPeriodic, double cutoffDistance, const vector<vector<int> >& exclusionList, const map<string, string>& defines) {
     if (cutoff != -1.0) {
         if (usesCutoff != useCutoff)
             throw OpenMMException("All Forces must agree on whether to use a cutoff");
@@ -86,12 +86,18 @@ void OpenCLNonbondedUtilities::addInteraction(bool usesCutoff, bool usesPeriodic
         }
         if (!sameExclusions)
             throw OpenMMException("All Forces must have identical exceptions");
+        for (map<string, string>::const_iterator iter = defines.begin(); iter != defines.end(); ++iter) {
+            map<string, string>::const_iterator existing = defines.find(iter->first);
+            if (existing != defines.end() && existing->second != iter->second)
+                throw OpenMMException("Two Forces define different values for "+iter->first);
+        }
     }
     else {
         useCutoff = usesCutoff;
         usePeriodic = usesPeriodic;
         cutoff = cutoffDistance;
         atomExclusions = exclusionList;
+        kernelDefines.insert(defines.begin(), defines.end());
     }
 }
 
@@ -116,9 +122,13 @@ void OpenCLNonbondedUtilities::initialize(const System& system) {
 
     // Create kernels.
 
-    map<string, string> defines;
+    map<string, string> defines = kernelDefines;
     if (forceBufferPerAtomBlock)
-        defines["USE_OUTPUT_BUFFER_PER_BLOCK"] = "true";
+        defines["USE_OUTPUT_BUFFER_PER_BLOCK"] = "1";
+    if (useCutoff)
+        defines["USE_CUTOFF"] = "1";
+    if (usePeriodic)
+        defines["USE_PERIODIC"] = "1";
     cl::Program forceProgram = context.createProgram(context.loadSourceFromFile("nonbonded.cl"), defines);
     forceKernel = cl::Kernel(forceProgram, "computeNonbonded");
     cl::Program interactingBlocksProgram = context.createProgram(context.loadSourceFromFile("findInteractingBlocks.cl"), defines);
@@ -264,19 +274,25 @@ void OpenCLNonbondedUtilities::computeInteractions() {
     hasComputedInteractions = true;
     forceKernel.setArg<cl_int>(0, tiles->getSize());
     forceKernel.setArg<cl_int>(1, context.getPaddedNumAtoms());
-    forceKernel.setArg<cl_float>(2, cutoff*cutoff);
-    forceKernel.setArg<mm_float4>(3, periodicBoxSize);
-    forceKernel.setArg<cl::Buffer>(4, context.getForceBuffers().getDeviceBuffer());
-    forceKernel.setArg<cl::Buffer>(5, context.getEnergyBuffer().getDeviceBuffer());
-    forceKernel.setArg<cl::Buffer>(6, context.getPosq().getDeviceBuffer());
-    forceKernel.setArg<cl::Buffer>(7, tiles->getDeviceBuffer());
-    forceKernel.setArg<cl::Buffer>(8, exclusions->getDeviceBuffer());
-    forceKernel.setArg<cl::Buffer>(9, exclusionIndex->getDeviceBuffer());
-    forceKernel.setArg(10, OpenCLContext::ThreadBlockSize*sizeof(cl_float4), NULL);
-    forceKernel.setArg(11, OpenCLContext::ThreadBlockSize*sizeof(cl_float4), NULL);
+    forceKernel.setArg<cl::Buffer>(2, context.getForceBuffers().getDeviceBuffer());
+    forceKernel.setArg<cl::Buffer>(3, context.getEnergyBuffer().getDeviceBuffer());
+    forceKernel.setArg<cl::Buffer>(4, context.getPosq().getDeviceBuffer());
+    forceKernel.setArg<cl::Buffer>(5, tiles->getDeviceBuffer());
+    forceKernel.setArg<cl::Buffer>(6, exclusions->getDeviceBuffer());
+    forceKernel.setArg<cl::Buffer>(7, exclusionIndex->getDeviceBuffer());
+    forceKernel.setArg(8, OpenCLContext::ThreadBlockSize*sizeof(cl_float4), NULL);
+    forceKernel.setArg(9, OpenCLContext::ThreadBlockSize*sizeof(cl_float4), NULL);
+    int paramBase = 10;
+    if (useCutoff) {
+        paramBase = 14;
+        forceKernel.setArg<cl_float>(10, cutoff*cutoff);
+        forceKernel.setArg<mm_float4>(11, periodicBoxSize);
+        forceKernel.setArg<cl::Buffer>(12, interactionFlags->getDeviceBuffer());
+        forceKernel.setArg(13, OpenCLContext::ThreadBlockSize*sizeof(cl_float4), NULL);
+    }
     for (int i = 0; i < (int) parameters.size(); i++) {
-        forceKernel.setArg<cl::Buffer>(i*2+12, *parameters[i].buffer);
-        forceKernel.setArg(i*2+13, OpenCLContext::ThreadBlockSize*parameters[i].size, NULL);
+        forceKernel.setArg<cl::Buffer>(i*2+paramBase, *parameters[i].buffer);
+        forceKernel.setArg(i*2+paramBase+1, OpenCLContext::ThreadBlockSize*parameters[i].size, NULL);
     }
     context.executeKernel(forceKernel, tiles->getSize()*OpenCLContext::TileSize);
 }

platforms/opencl/src/OpenCLNonbondedUtilities.h

@@ -52,8 +52,9 @@ public:
      * @param usesPeriodic   specifies whether periodic boundary conditions should be applied to this interaction
      * @param cutoffDistance the cutoff distance for this interaction (ignored if usesCutoff is false)
      * @param exclusionList  for each atom, specifies the list of other atoms whose interactions should be excluded
+     * @param defines        preprocessor macros to define when compiling the kernel
      */
-    void addInteraction(bool usesCutoff, bool usesPeriodic, double cutoffDistance, const std::vector<std::vector<int> >& exclusionList);
+    void addInteraction(bool usesCutoff, bool usesPeriodic, double cutoffDistance, const std::vector<std::vector<int> >& exclusionList, const std::map<std::string, std::string>& defines);
     /**
      * Add a per-atom parameter that interactions may depend on.
      *
@@ -142,6 +143,7 @@ private:
     std::vector<std::vector<int> > atomExclusions;
     std::vector<ParameterInfo> parameters;
     OpenCLCompact* compact;
+    std::map<std::string, std::string> kernelDefines;
     double cutoff;
     bool useCutoff, usePeriodic, forceBufferPerAtomBlock, hasComputedInteractions;
     int numForceBuffers;

platforms/opencl/src/kernels/findInteractingBlocks.cl

@@ -118,6 +118,7 @@ __kernel void findInteractionsWithinBlocks(float cutoffSquared, float4 periodicB
 
             // Sum the flags.
 
+#ifdef WARPS_ARE_ATOMIC
             if (index % 2 == 0)
                 flags[thread] += flags[thread+1];
             if (index % 4 == 0)
@@ -126,6 +127,21 @@ __kernel void findInteractionsWithinBlocks(float cutoffSquared, float4 periodicB
                 flags[thread] += flags[thread+4];
             if (index % 16 == 0)
                 flags[thread] += flags[thread+8];
+#else
+            barrier(CLK_LOCAL_MEM_FENCE);
+            if (index % 2 == 0)
+                flags[thread] += flags[thread+1];
+            barrier(CLK_LOCAL_MEM_FENCE);
+            if (index % 4 == 0)
+                flags[thread] += flags[thread+2];
+            barrier(CLK_LOCAL_MEM_FENCE);
+            if (index % 8 == 0)
+                flags[thread] += flags[thread+4];
+            barrier(CLK_LOCAL_MEM_FENCE);
+            if (index % 16 == 0)
+                flags[thread] += flags[thread+8];
+            barrier(CLK_LOCAL_MEM_FENCE);
+#endif
             if (index == 0)
             {
                 unsigned int allFlags = flags[thread] + flags[thread+16];

platforms/opencl/src/kernels/nonbonded.cl

@@ -5,18 +5,18 @@ const float EpsilonFactor = 138.935485f;
  * Compute nonbonded interactions.
  */
 
-__kernel void computeNonbonded(int numTiles, int paddedNumAtoms, float cutoffSquared, float4 periodicBoxSize,
+__kernel void computeNonbonded(int numTiles, int paddedNumAtoms,
         __global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global unsigned int* tiles,
         __global unsigned int* exclusions,  __global unsigned int* exclusionIndices, __local float4* local_posq, __local float4* local_force,
+#ifdef USE_CUTOFF
+        float cutoffSquared, float4 periodicBoxSize, __global unsigned int* interactionFlags, __local float4* tempBuffer,
+#endif
         __global float2* sigmaEpsilon, __local float2* local_sigmaEpsilon) {
     unsigned int totalWarps = get_global_size(0)/TileSize;
     unsigned int warp = get_global_id(0)/TileSize;
     unsigned int pos = warp*numTiles/totalWarps;
     unsigned int end = (warp+1)*numTiles/totalWarps;
     float energy = 0.0f;
-#ifdef USE_CUTOFF
-    float3* tempBuffer = (float3*) &sA[cSim.nonbond_threads_per_block];
-#endif
     unsigned int lasty = 0xFFFFFFFF;
 
     while (pos < end) {
@@ -61,8 +61,8 @@ __kernel void computeNonbonded(int numTiles, int paddedNumAtoms, float cutoffSqu
                 float dEdR = eps * (12.0f * sig6 - 6.0f) * sig6;
                 float tempEnergy = eps * (sig6 - 1.0f) * sig6;
 #ifdef USE_CUTOFF
-                dEdR += apos.w * local_posq[tbx+j].w * (invR - 2.0f * cSim.reactionFieldK * r2);
-                tempEnergy += apos.w * local_posq[tbx+j].w * (invR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
+                dEdR += apos.w * local_posq[tbx+j].w * (invR - 2.0f * REACTION_FIELD_K * r2);
+                tempEnergy += apos.w * local_posq[tbx+j].w * (invR + REACTION_FIELD_K * r2 - REACTION_FIELD_C);
 #else
                 dEdR += apos.w * local_posq[tbx+j].w * invR;
                 tempEnergy += apos.w * local_posq[tbx+j].w * invR;
@@ -108,7 +108,7 @@ __kernel void computeNonbonded(int numTiles, int paddedNumAtoms, float cutoffSqu
             local_force[get_local_id(0)] = 0.0f;
             apos.w *= EpsilonFactor;
 #ifdef USE_CUTOFF
-            unsigned int flags = cSim.pInteractionFlag[pos];
+            unsigned int flags = interactionFlags[pos];
             if (!hasExclusions && flags != 0xFFFFFFFF) {
                 if (flags == 0) {
                     // No interactions in this tile.
@@ -131,12 +131,12 @@ __kernel void computeNonbonded(int numTiles, int paddedNumAtoms, float cutoffSqu
                             float sig2 = invR * sig;
                             sig2 *= sig2;
                             float sig6 = sig2 * sig2 * sig2;
-                            float eps = a.y * local_sigmaEpsilon[tbx+j].Y;
+                            float eps = a.y * local_sigmaEpsilon[tbx+j].y;
                             float dEdR = eps * (12.0f * sig6 - 6.0f) * sig6;
 			    float tempEnergy = eps * (sig6 - 1.0f) * sig6;
 #ifdef USE_CUTOFF
-                            dEdR += apos.w * local_posq[tbx+j].w * (invR - 2.0f * cSim.reactionFieldK * r2);
-                            tempEnergy += apos.w * local_posq[tbx+j].w * (invR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
+                            dEdR += apos.w * local_posq[tbx+j].w * (invR - 2.0f * REACTION_FIELD_K * r2);
+                            tempEnergy += apos.w * local_posq[tbx+j].w * (invR + REACTION_FIELD_K * r2 - REACTION_FIELD_C);
 #else
                             dEdR += apos.w * local_posq[tbx+j].w * invR;
                             tempEnergy += apos.w * local_posq[tbx+j].w * invR;
@@ -188,19 +188,19 @@ __kernel void computeNonbonded(int numTiles, int paddedNumAtoms, float cutoffSqu
 #endif
                     float r2 = delta.x * delta.x + delta.y * delta.y + delta.z * delta.z;
                     float invR = 1.0f / sqrt(r2);
-                    float sig = a.x + local_sigmaEpsilon[tbx+j].x;
+                    float sig = a.x + local_sigmaEpsilon[tbx+tj].x;
                     float sig2 = invR * sig;
                     sig2 *= sig2;
                     float sig6 = sig2 * sig2 * sig2;
-                    float eps = a.y * local_sigmaEpsilon[tbx+j].y;
+                    float eps = a.y * local_sigmaEpsilon[tbx+tj].y;
                     float dEdR = eps * (12.0f * sig6 - 6.0f) * sig6;
 		    float tempEnergy = eps * (sig6 - 1.0f) * sig6;
 #ifdef USE_CUTOFF
-                    dEdR += apos.w * local_posq[tbx+j].w * (invR - 2.0f * cSim.reactionFieldK * r2);
-		    tempEnergy += apos.w * local_posq[tbx+j].w * (invR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
+                    dEdR += apos.w * local_posq[tbx+tj].w * (invR - 2.0f * REACTION_FIELD_K * r2);
+		    tempEnergy += apos.w * local_posq[tbx+tj].w * (invR + REACTION_FIELD_K * r2 - REACTION_FIELD_C);
 #else
-                    dEdR += apos.w * local_posq[tbx+j].w * invR;
-                    tempEnergy += apos.w * local_posq[tbx+j].w * invR;
+                    dEdR += apos.w * local_posq[tbx+tj].w * invR;
+                    tempEnergy += apos.w * local_posq[tbx+tj].w * invR;
 #endif
                     dEdR *= invR * invR;
 #ifdef USE_CUTOFF

platforms/opencl/src/kernels/nonbondedExceptions.cl

@@ -28,8 +28,8 @@ __kernel void computeNonbondedExceptions(int numAtoms, int numExceptions, float
         float dEdR = exceptionParams.z*(12.0f*sig6-6.0f)*sig6;
         float tempEnergy = exceptionParams.z*(sig6-1.0f)*sig6;
 #ifdef USE_CUTOFF
-        dEdR += exceptionParams.x*(invR-2.0f*cSim.reactionFieldK*r2);
-        tempEnergy += exceptionParams.x*(invR+cSim.reactionFieldK*r2-cSim.reactionFieldC);
+        dEdR += exceptionParams.x*(invR-2.0f*REACTION_FIELD_K*r2);
+        tempEnergy += exceptionParams.x*(invR+REACTION_FIELD_K*r2-REACTION_FIELD_C);
 #else
         dEdR += exceptionParams.x*invR;
         tempEnergy += exceptionParams.x*invR;

platforms/opencl/tests/TestOpenCLNonbondedForce.cpp

@@ -525,7 +525,6 @@ void testBlockInteractions(bool periodic) {
     const unsigned int atoms = clcontext.getPaddedNumAtoms();
     const unsigned int grid = OpenCLContext::TileSize;
     const unsigned int dim = clcontext.getNumAtomBlocks();
-    printf("%d of %d\n", numWithInteractions, hasInteractions.size());
     for (int i = 0; i < numWithInteractions; i++) {
         unsigned int tile = interactingTiles[i];
         unsigned int x = (tile >> 17);
@@ -615,12 +614,12 @@ int main() {
         testCoulomb();
         testLJ();
         testExclusionsAnd14();
-//        testCutoff();
-//        testCutoff14();
-//        testPeriodic();
-//        testLargeSystem();
+        testCutoff();
+        testCutoff14();
+        testPeriodic();
+        testLargeSystem();
         testBlockInteractions(false);
-//        testBlockInteractions(true);
+        testBlockInteractions(true);
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;