Continuing implementation of nonbonded forces

platforms/opencl/src/OpenCLArray.h

@@ -54,7 +54,14 @@ public:
      */
     OpenCLArray(OpenCLContext& context, int size, const std::string& name, bool createHostBuffer = false) :
             context(context), size(size), name(name), local(createHostBuffer ? size : 0), ownsBuffer(true) {
-        buffer = new cl::Buffer(context.getContext(), CL_MEM_READ_WRITE, size*sizeof(T));
+        try {
+            buffer = new cl::Buffer(context.getContext(), CL_MEM_READ_WRITE, size*sizeof(T));
+        }
+        catch (cl::Error err) {
+            std::stringstream str;
+            str<<"Error creating array "<<name<<": "<<err.what()<<" ("<<err.err()<<")";
+            throw OpenMMException(str.str());
+        }
     }
     /**
      * Create an OpenCLArray object the uses a preexisting Buffer.

platforms/opencl/src/OpenCLForceInfo.cpp

@@ -29,7 +29,7 @@
 using namespace OpenMM;
 using namespace std;
 
-bool OpenCLForceInfo::areParticlesIdentical() {
+bool OpenCLForceInfo::areParticlesIdentical(int particle1, int particle2) {
     return true;
 }
 

platforms/opencl/src/OpenCLForceInfo.h

@@ -38,8 +38,7 @@ namespace OpenMM {
 
 class OpenCLForceInfo {
 public:
-    OpenCLForceInfo(int requiredForceBuffers, bool usesNeighborList, double cutoffDistance) :
-        requiredForceBuffers(requiredForceBuffers), usesNeighborList(usesNeighborList), cutoffDistance(cutoffDistance) {
+    OpenCLForceInfo(int requiredForceBuffers) : requiredForceBuffers(requiredForceBuffers) {
     }
     /**
      * Get the number of force buffers this force requires.
@@ -47,23 +46,10 @@ public:
     int getRequiredForceBuffers() {
         return requiredForceBuffers;
     }
-    /**
-     * Get whether or not this force makes use of the neighbor list.
-     */
-    bool getUsesNeighborList() {
-        return usesNeighborList;
-    }
-    /**
-     * Get the cutoff distance used by this force.  If usesNeighborList() returns false,
-     * the return value from this method is ignored.
-     */
-    double getCutoffDistance() {
-        return cutoffDistance;
-    }
     /**
      * Get whether or not two particles have identical force field parameters.
      */
-    virtual bool areParticlesIdentical();
+    virtual bool areParticlesIdentical(int particle1, int particle2);
     /**
      * Get the number of particle groups defined by this force.
      */
@@ -78,8 +64,6 @@ public:
     virtual bool areGroupsIdentical(int group1, int group2);
 private:
     int requiredForceBuffers;
-    bool usesNeighborList;
-    double cutoffDistance;
 };
 
 } // namespace OpenMM

platforms/opencl/src/OpenCLKernels.cpp

@@ -57,6 +57,7 @@ void OpenCLCalcForcesAndEnergyKernel::finishForceComputation(ContextImpl& contex
 
 void OpenCLCalcForcesAndEnergyKernel::beginEnergyComputation(ContextImpl& context) {
     cl.clearBuffer(cl.getEnergyBuffer());
+    cl.getNonbondedUtilties().prepareInteractions();
 }
 
 double OpenCLCalcForcesAndEnergyKernel::finishEnergyComputation(ContextImpl& context) {
@@ -149,7 +150,7 @@ void OpenCLUpdateStateDataKernel::getForces(ContextImpl& context, std::vector<Ve
 
 class OpenCLBondForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLBondForceInfo(int requiredBuffers, const HarmonicBondForce& force) : OpenCLForceInfo(requiredBuffers, false, 0.0), force(force) {
+    OpenCLBondForceInfo(int requiredBuffers, const HarmonicBondForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumBonds();
@@ -198,9 +199,8 @@ void OpenCLCalcHarmonicBondForceKernel::initialize(const System& system, const H
     params->upload(paramVector);
     indices->upload(indicesVector);
     int maxBuffers = 1;
-    for (int i = 0; i < forceBufferCounter.size(); i++) {
+    for (int i = 0; i < forceBufferCounter.size(); i++)
         maxBuffers = max(maxBuffers, forceBufferCounter[i]);
-    }
     cl.addForce(new OpenCLBondForceInfo(maxBuffers, force));
     cl::Program program = cl.createProgram(cl.loadSourceFromFile("harmonicBondForce.cl"));
     kernel = cl::Kernel(program, "calcHarmonicBondForce");
@@ -224,7 +224,7 @@ double OpenCLCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
 
 class OpenCLAngleForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLAngleForceInfo(int requiredBuffers, const HarmonicAngleForce& force) : OpenCLForceInfo(requiredBuffers, false, 0.0), force(force) {
+    OpenCLAngleForceInfo(int requiredBuffers, const HarmonicAngleForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumAngles();
@@ -275,9 +275,8 @@ void OpenCLCalcHarmonicAngleForceKernel::initialize(const System& system, const
     params->upload(paramVector);
     indices->upload(indicesVector);
     int maxBuffers = 1;
-    for (int i = 0; i < forceBufferCounter.size(); i++) {
+    for (int i = 0; i < forceBufferCounter.size(); i++)
         maxBuffers = max(maxBuffers, forceBufferCounter[i]);
-    }
     cl.addForce(new OpenCLAngleForceInfo(maxBuffers, force));
     cl::Program program = cl.createProgram(cl.loadSourceFromFile("harmonicAngleForce.cl"));
     kernel = cl::Kernel(program, "calcHarmonicAngleForce");
@@ -301,7 +300,7 @@ double OpenCLCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
 
 class OpenCLPeriodicTorsionForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLPeriodicTorsionForceInfo(int requiredBuffers, const PeriodicTorsionForce& force) : OpenCLForceInfo(requiredBuffers, false, 0.0), force(force) {
+    OpenCLPeriodicTorsionForceInfo(int requiredBuffers, const PeriodicTorsionForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumTorsions();
@@ -353,9 +352,8 @@ void OpenCLCalcPeriodicTorsionForceKernel::initialize(const System& system, cons
     params->upload(paramVector);
     indices->upload(indicesVector);
     int maxBuffers = 1;
-    for (int i = 0; i < forceBufferCounter.size(); i++) {
+    for (int i = 0; i < forceBufferCounter.size(); i++)
         maxBuffers = max(maxBuffers, forceBufferCounter[i]);
-    }
     cl.addForce(new OpenCLPeriodicTorsionForceInfo(maxBuffers, force));
     cl::Program program = cl.createProgram(cl.loadSourceFromFile("periodicTorsionForce.cl"));
     kernel = cl::Kernel(program, "calcPeriodicTorsionForce");
@@ -379,7 +377,7 @@ double OpenCLCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context)
 
 class OpenCLRBTorsionForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLRBTorsionForceInfo(int requiredBuffers, const RBTorsionForce& force) : OpenCLForceInfo(requiredBuffers, false, 0.0), force(force) {
+    OpenCLRBTorsionForceInfo(int requiredBuffers, const RBTorsionForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumTorsions();
@@ -431,9 +429,8 @@ void OpenCLCalcRBTorsionForceKernel::initialize(const System& system, const RBTo
     params->upload(paramVector);
     indices->upload(indicesVector);
     int maxBuffers = 1;
-    for (int i = 0; i < forceBufferCounter.size(); i++) {
+    for (int i = 0; i < forceBufferCounter.size(); i++)
         maxBuffers = max(maxBuffers, forceBufferCounter[i]);
-    }
     cl.addForce(new OpenCLRBTorsionForceInfo(maxBuffers, force));
     cl::Program program = cl.createProgram(cl.loadSourceFromFile("rbTorsionForce.cl"));
     kernel = cl::Kernel(program, "calcRBTorsionForce");
@@ -455,9 +452,45 @@ double OpenCLCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
     return 0.0;
 }
 
+class OpenCLNonbondedForceInfo : public OpenCLForceInfo {
+public:
+    OpenCLNonbondedForceInfo(int requiredBuffers, const NonbondedForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
+    }
+    bool areParticlesIdentical(int particle1, int particle2) {
+        double charge1, charge2, sigma1, sigma2, epsilon1, epsilon2;
+        force.getParticleParameters(particle1, charge1, sigma1, epsilon1);
+        force.getParticleParameters(particle2, charge2, sigma2, epsilon2);
+        return (charge1 == charge2 && sigma1 == sigma2 && epsilon1 == epsilon2);
+    }
+    int getNumParticleGroups() {
+        return force.getNumExceptions();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2;
+        double chargeProd, sigma, epsilon;
+        force.getExceptionParameters(index, particle1, particle2, chargeProd, sigma, epsilon);
+        particles.resize(2);
+        particles[0] = particle1;
+        particles[1] = particle2;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2;
+        double chargeProd1, chargeProd2, sigma1, sigma2, epsilon1, epsilon2;
+        force.getExceptionParameters(group1, particle1, particle2, chargeProd1, sigma1, epsilon1);
+        force.getExceptionParameters(group2, particle1, particle2, chargeProd2, sigma2, epsilon2);
+        return (chargeProd1 == chargeProd2 && sigma1 == sigma2 && epsilon1 == epsilon2);
+    }
+private:
+    const NonbondedForce& force;
+};
+
 OpenCLCalcNonbondedForceKernel::~OpenCLCalcNonbondedForceKernel() {
     if (sigmaEpsilon != NULL)
         delete sigmaEpsilon;
+    if (exceptionParams != NULL)
+        delete exceptionParams;
+    if (exceptionIndices != NULL)
+        delete exceptionIndices;
 }
 
 void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {
@@ -535,6 +568,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
 //    gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList, method);
     cl.getNonbondedUtilties().addInteraction(useCutoff, usePeriodic, force.getCutoffDistance(), exclusionList);
     cl.getNonbondedUtilties().addParameter("sigmaEpsilon", "float2", 8, sigmaEpsilon->getDeviceBuffer());
+    cutoffSquared = force.getCutoffDistance()*force.getCutoffDistance();
 
     // Compute the Ewald self energy.
 
@@ -545,30 +579,48 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
 //                ewaldSelfEnergy -= selfEnergyScale*q[i]*q[i];
     }
 
-    // Initialize 1-4 nonbonded interactions.
+    // Initialize the exceptions.
 
-    {
-        int numExceptions = exceptions.size();
-        vector<int> particle1(numExceptions);
-        vector<int> particle2(numExceptions);
-        vector<float> c6(numExceptions);
-        vector<float> c12(numExceptions);
-        vector<float> q1(numExceptions);
-        vector<float> q2(numExceptions);
+    int numExceptions = exceptions.size();
+    int maxBuffers = 1;
+    if (numExceptions > 0) {
+        exceptionParams = new OpenCLArray<mm_float4>(cl, numExceptions, "exceptionParams");
+        exceptionIndices = new OpenCLArray<mm_int4>(cl, numExceptions, "exceptionIndices");
+        vector<mm_float4> exceptionParamsVector(numExceptions);
+        vector<mm_int4> exceptionIndicesVector(numExceptions);
+        vector<int> forceBufferCounter(system.getNumParticles(), 0);
         for (int i = 0; i < numExceptions; i++) {
-            double charge, sig, eps;
-            force.getExceptionParameters(exceptions[i], particle1[i], particle2[i], charge, sig, eps);
-            c6[i] = (float) (4*eps*pow(sig, 6.0));
-            c12[i] = (float) (4*eps*pow(sig, 12.0));
-            q1[i] = (float) charge;
-            q2[i] = 1.0f;
+            int particle1, particle2;
+            double chargeProd, sigma, epsilon;
+            force.getExceptionParameters(exceptions[i], particle1, particle2, chargeProd, sigma, epsilon);
+            exceptionParamsVector[i] = (mm_float4) {(float) (138.935485*chargeProd), (float) sigma, (float) (4.0*epsilon), 0.0f};
+            exceptionIndicesVector[i] = (mm_int4) {particle1, particle2, forceBufferCounter[particle1]++, forceBufferCounter[particle2]++};
         }
-//        gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, particle1, particle2, c6, c12, q1, q2);
+        exceptionParams->upload(exceptionParamsVector);
+        exceptionIndices->upload(exceptionIndicesVector);
+        for (int i = 0; i < forceBufferCounter.size(); i++)
+            maxBuffers = max(maxBuffers, forceBufferCounter[i]);
     }
+    cl.addForce(new OpenCLNonbondedForceInfo(maxBuffers, force));
+    cl::Program program = cl.createProgram(cl.loadSourceFromFile("nonbondedExceptions.cl"));
+    exceptionsKernel = cl::Kernel(program, "computeNonbondedExceptions");
 }
 
 void OpenCLCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
     cl.getNonbondedUtilties().computeInteractions();
+    if (exceptionIndices != NULL) {
+        int numExceptions = exceptionIndices->getSize();
+        exceptionsKernel.setArg<cl_int>(0, cl.getPaddedNumAtoms());
+        exceptionsKernel.setArg<cl_int>(1, numExceptions);
+        exceptionsKernel.setArg<cl_float>(2, cutoffSquared);
+        exceptionsKernel.setArg<mm_float4>(3, cl.getNonbondedUtilties().getPeriodicBoxSize());
+        exceptionsKernel.setArg<cl::Buffer>(4, cl.getForceBuffers().getDeviceBuffer());
+        exceptionsKernel.setArg<cl::Buffer>(5, cl.getEnergyBuffer().getDeviceBuffer());
+        exceptionsKernel.setArg<cl::Buffer>(6, cl.getPosq().getDeviceBuffer());
+        exceptionsKernel.setArg<cl::Buffer>(7, exceptionParams->getDeviceBuffer());
+        exceptionsKernel.setArg<cl::Buffer>(8, exceptionIndices->getDeviceBuffer());
+        cl.executeKernel(exceptionsKernel, numExceptions);
+    }
 }
 
 double OpenCLCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {

platforms/opencl/src/OpenCLKernels.h

@@ -299,7 +299,8 @@ private:
  */
 class OpenCLCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
 public:
-    OpenCLCalcNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcNonbondedForceKernel(name, platform), cl(cl) {
+    OpenCLCalcNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcNonbondedForceKernel(name, platform), cl(cl),
+            exceptionParams(NULL), exceptionIndices(NULL) {
     }
     ~OpenCLCalcNonbondedForceKernel();
     /**
@@ -325,7 +326,10 @@ public:
 private:
     OpenCLContext& cl;
     OpenCLArray<mm_float2>* sigmaEpsilon;
-    double ewaldSelfEnergy;
+    OpenCLArray<mm_float4>* exceptionParams;
+    OpenCLArray<mm_int4>* exceptionIndices;
+    cl::Kernel exceptionsKernel;
+    double cutoffSquared, ewaldSelfEnergy;
 };
 
 ///**

platforms/opencl/src/OpenCLNonbondedUtilities.h

@@ -71,6 +71,12 @@ public:
     int getNumForceBuffers() {
         return numForceBuffers;
     }
+    /**
+     * Get the periodic box size.
+     */
+    mm_float4 getPeriodicBoxSize() {
+        return periodicBoxSize;
+    }
     /**
      * Prepare to compute interactions.  This updates the neighbor list.
      */

platforms/opencl/src/kernels/nonbonded.cl

@@ -39,7 +39,6 @@ __kernel void computeNonbonded(int numTiles, int paddedNumAtoms, float cutoffSqu
 
             local_posq[get_local_id(0)] = apos;
             local_sigmaEpsilon[get_local_id(0)] = a;
-            barrier(CLK_LOCAL_MEM_FENCE);
             apos.w *= EpsilonFactor;
             unsigned int xi = x/TileSize;
             unsigned int tile = xi+xi*paddedNumAtoms/TileSize-xi*(xi+1)/2;
@@ -107,7 +106,6 @@ __kernel void computeNonbonded(int numTiles, int paddedNumAtoms, float cutoffSqu
                 local_sigmaEpsilon[get_local_id(0)] = sigmaEpsilon[j];
             }
             local_force[get_local_id(0)] = 0.0f;
-            barrier(CLK_LOCAL_MEM_FENCE);
             apos.w *= EpsilonFactor;
 #ifdef USE_CUTOFF
             unsigned int flags = cSim.pInteractionFlag[pos];
@@ -159,19 +157,14 @@ __kernel void computeNonbonded(int numTiles, int paddedNumAtoms, float cutoffSqu
 
                             if (tgx % 2 == 0)
                                 tempBuffer[get_local_id(0)].xyz += tempBuffer[get_local_id(0)+1].xyz;
-                            barrier(CLK_LOCAL_MEM_FENCE);
                             if (tgx % 4 == 0)
                                 tempBuffer[get_local_id(0)].xyz += tempBuffer[get_local_id(0)+2].xyz;
-                            barrier(CLK_LOCAL_MEM_FENCE);
                             if (tgx % 8 == 0)
                                 tempBuffer[get_local_id(0)].xyz += tempBuffer[get_local_id(0)+4].xyz;
-                            barrier(CLK_LOCAL_MEM_FENCE);
                             if (tgx % 16 == 0)
                                 tempBuffer[get_local_id(0)].xyz += tempBuffer[get_local_id(0)+8].xyz;
-                            barrier(CLK_LOCAL_MEM_FENCE);
                             if (tgx == 0)
                                 local_force[tbx+j].xyz += tempBuffer[get_local_id(0)].xyz + tempBuffer[get_local_id(0)+16].xyz;
-                            barrier(CLK_LOCAL_MEM_FENCE);
                         }
                     }
                 }
@@ -222,7 +215,6 @@ __kernel void computeNonbonded(int numTiles, int paddedNumAtoms, float cutoffSqu
                     delta.xyz *= dEdR;
                     af.xyz -= delta.xyz;
                     local_force[tbx+tj].xyz += delta.xyz;
-                    barrier(CLK_LOCAL_MEM_FENCE);
                     excl >>= 1;
                     tj = (tj + 1) & (TileSize - 1);
                 }

platforms/opencl/src/kernels/nonbondedExceptions.cl

+/**
+ * Compute nonbonded exceptions.
+ */
+
+__kernel void computeNonbondedExceptions(int numAtoms, int numExceptions, float cutoffSquared, float4 periodicBoxSize, __global float4* forceBuffers, __global float* energyBuffer,
+        __global float4* posq, __global float4* params, __global int4* indices) {
+    int index = get_global_id(0);
+    float energy = 0.0f;
+    while (index < numExceptions) {
+        // Look up the data for this bonds.
+
+        int4 atoms = indices[index];
+        float4 exceptionParams = params[index];
+        float4 delta = posq[atoms.y]-posq[atoms.x];
+#ifdef USE_PERIODIC
+        delta.x -= floor(delta.x/periodicBoxSize.x+0.5f)*periodicBoxSize.x;
+        delta.y -= floor(delta.y/periodicBoxSize.y+0.5f)*periodicBoxSize.y;
+        delta.z -= floor(delta.z/periodicBoxSize.z+0.5f)*periodicBoxSize.z;
+#endif
+
+        // Compute the force.
+
+        float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
+        float invR = 1.0f/sqrt(r2);
+        float sig2 = invR*exceptionParams.y;
+        sig2 *= sig2;
+        float sig6 = sig2*sig2*sig2;
+        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);
+#else
+        dEdR += exceptionParams.x*invR;
+        tempEnergy += exceptionParams.x*invR;
+#endif
+        dEdR *= invR*invR;
+#ifdef USE_CUTOFF
+        if (r2 > cutoffSquared) {
+            dEdR = 0.0f;
+            tempEnergy  = 0.0f;
+        }
+#endif
+        energy += tempEnergy;
+        delta.xyz *= dEdR;
+
+        // Record the force on each of the two atoms.
+
+        unsigned int offsetA = atoms.x+atoms.z*numAtoms;
+        unsigned int offsetB = atoms.y+atoms.w*numAtoms;
+        float4 forceA = forceBuffers[offsetA];
+        float4 forceB = forceBuffers[offsetB];
+        forceA.xyz -= delta.xyz;
+        forceB.xyz += delta.xyz;
+        forceBuffers[offsetA] = forceA;
+        forceBuffers[offsetB] = forceB;
+        index += get_global_size(0);
+    }
+    energyBuffer[get_global_id(0)] += energy;
+}

platforms/opencl/tests/TestOpenCLLangevinIntegrator.cpp

@@ -105,7 +105,7 @@ void testTemperature() {
         system.addParticle(2.0);
         forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
     }
-//    system.addForce(forceField);
+    system.addForce(forceField);
     Context context(system, integrator, platform);
     vector<Vec3> positions(numParticles);
     for (int i = 0; i < numParticles; ++i)
@@ -188,7 +188,7 @@ void testRandomSeed() {
         system.addParticle(2.0);
         forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
     }
-//    system.addForce(forceField);
+    system.addForce(forceField);
     vector<Vec3> positions(numParticles);
     vector<Vec3> velocities(numParticles);
     for (int i = 0; i < numParticles; ++i) {

platforms/opencl/tests/TestOpenCLNonbondedForce.cpp

@@ -604,7 +604,7 @@ int main() {
     try {
         testCoulomb();
         testLJ();
-//        testExclusionsAnd14();
+        testExclusionsAnd14();
 //        testCutoff();
 //        testCutoff14();
 //        testPeriodic();

platforms/opencl/tests/TestOpenCLSettle.cpp

@@ -68,7 +68,7 @@ void testConstraints() {
         system.addConstraint(i*3, i*3+2, 0.1);
         system.addConstraint(i*3+1, i*3+2, 0.163);
     }
-//    system.addForce(forceField);
+    system.addForce(forceField);
     Context context(system, integrator, platform);
     vector<Vec3> positions(numParticles);
     vector<Vec3> velocities(numParticles);