Continuing to implement OpenCL version of GayBerneForce

openmmapi/src/GayBerneForceImpl.cpp

@@ -115,6 +115,7 @@ void GayBerneForceImpl::initialize(ContextImpl& context) {
 double GayBerneForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
     if ((groups&(1<<owner.getForceGroup())) != 0)
         return kernel.getAs<CalcGayBerneForceKernel>().execute(context, includeForces, includeEnergy);
+    return 0.0;
 }
 
 std::vector<std::string> GayBerneForceImpl::getKernelNames() {

platforms/opencl/src/OpenCLKernels.cpp

@@ -6154,9 +6154,7 @@ void OpenCLCalcGayBerneForceKernel::initialize(const System& system, const GayBe
             defines["SWITCH_C5"] = cl.doubleToString(6/pow(force.getSwitchingDistance()-cutoff, 5.0));
         }
     }
-    if (!cl.getSupports64BitGlobalAtomics()) {
     defines["PADDED_NUM_ATOMS"] = cl.intToString(cl.getPaddedNumAtoms());
-    }
     cl::Program program = cl.createProgram(OpenCLKernelSources::gayBerne, defines);
     framesKernel = cl::Kernel(program, "computeEllipsoidFrames");
     blockBoundsKernel = cl::Kernel(program, "findBlockBounds");
@@ -6234,12 +6232,37 @@ double OpenCLCalcGayBerneForceKernel::execute(ContextImpl& context, bool include
     cl.executeKernel(framesKernel, numRealParticles);
     setPeriodicBoxArgs(cl, blockBoundsKernel, 1);
     cl.executeKernel(blockBoundsKernel, (numRealParticles+31)/32);
-    if (nonbondedMethod != GayBerneForce::NoCutoff) {
+    if (nonbondedMethod == GayBerneForce::NoCutoff) {
+        cl.executeKernel(forceKernel, cl.getNonbondedUtilities().getNumForceThreadBlocks()*cl.getNonbondedUtilities().getForceThreadBlockSize());
+    }
+    else {
+        while (true) {
             setPeriodicBoxArgs(cl, neighborsKernel, 2);
             cl.executeKernel(neighborsKernel, numRealParticles);
+            cl_int* count = (cl_int*) cl.getPinnedBuffer();
+            cl::Event event;
+            cl.getQueue().enqueueReadBuffer(neighborBlockCount->getDeviceBuffer(), CL_FALSE, 0, neighborBlockCount->getSize()*neighborBlockCount->getElementSize(), count, NULL, &event);
             setPeriodicBoxArgs(cl, forceKernel, 20);
-    }
             cl.executeKernel(forceKernel, cl.getNonbondedUtilities().getNumForceThreadBlocks()*cl.getNonbondedUtilities().getForceThreadBlockSize());
+            event.wait();
+            if (*count <= maxNeighborBlocks)
+                break;
+            
+            // There wasn't enough room for the neighbor list, so we need to recreate it.
+
+            delete neighbors;
+            neighbors = NULL;
+            delete neighborIndex;
+            neighborIndex = NULL;
+            maxNeighborBlocks = (int) ceil((*count)*1.1);
+            neighbors = OpenCLArray::create<cl_int>(cl, maxNeighborBlocks*32, "neighbors");
+            neighborIndex = OpenCLArray::create<cl_int>(cl, maxNeighborBlocks, "neighbors");
+            neighborsKernel.setArg<cl::Buffer>(10, neighbors->getDeviceBuffer());
+            neighborsKernel.setArg<cl::Buffer>(11, neighborIndex->getDeviceBuffer());
+            forceKernel.setArg<cl::Buffer>(17, neighbors->getDeviceBuffer());
+            forceKernel.setArg<cl::Buffer>(18, neighborIndex->getDeviceBuffer());
+        }
+    }
     cl.executeKernel(torqueKernel, numRealParticles);
     return 0.0;
 }
@@ -6266,20 +6289,15 @@ void OpenCLCalcGayBerneForceKernel::copyParametersToContext(ContextImpl& context
     vector<mm_float4> sigParamsVector(cl.getPaddedNumAtoms(), mm_float4(0, 0, 0, 0));
     vector<mm_float2> epsParamsVector(cl.getPaddedNumAtoms(), mm_float2(0, 0));
     vector<mm_float4> scaleVector(cl.getPaddedNumAtoms(), mm_float4(0, 0, 0, 0));
-    numRealParticles = 0;
     for (int i = 0; i < force.getNumParticles(); i++) {
         int xparticle, yparticle;
         double sigma, epsilon, sx, sy, sz, ex, ey, ez;
         force.getParticleParameters(i, sigma, epsilon, xparticle, yparticle, sx, sy, sz, ex, ey, ez);
-        sigParamsVector[i] = mm_float4((float) (0.5*sigma), (float) (0.5*sx), (float) (0.5*sy), (float) (0.5*sz));
+        sigParamsVector[i] = mm_float4((float) (0.5*sigma), (float) (0.25*sx*sx), (float) (0.25*sy*sy), (float) (0.25*sz*sz));
         epsParamsVector[i] = mm_float2((float) sqrt(epsilon), (float) (0.125*(sx*sy + sz*sz)*sqrt(sx*sy)));
         scaleVector[i] = mm_float4((float) (1/sqrt(ex)), (float) (1/sqrt(ey)), (float) (1/sqrt(ez)), 0);
-        if (epsilon != 0.0) {
-            numRealParticles++;
-            isRealParticle[i] = true;
-        }
-        else
-            isRealParticle[i] = false;
+        if (epsilon != 0.0 && !isRealParticle[i])
+            throw OpenMMException("updateParametersInContext: The set of ignored particles (ones with epsilon=0) has changed");
     }
     sigParams->upload(sigParamsVector);
     epsParams->upload(epsParamsVector);

platforms/opencl/src/kernels/gayBerne.cl

 #pragma OPENCL EXTENSION cl_khr_global_int32_base_atomics : enable
+#ifdef SUPPORTS_64_BIT_ATOMICS
+#pragma OPENCL EXTENSION cl_khr_int64_base_atomics : enable
+#endif
 
 #define TILE_SIZE 32
 #define NEIGHBOR_BLOCK_SIZE 32
@@ -55,14 +58,8 @@ __kernel void computeEllipsoidFrames(int numParticles, __global const real4* res
         float3 e2 = scale[originalIndex].xyz;
         for (int i = 0; i < 3; i++)
             for (int j = 0; j < 3; j++) {
-                real belem = 0;
-                real gelem = 0;
-                for (int k = 0; k < 3; k++) {
-                    belem += a[k][i]*e2[k]*a[k][j];
-                    gelem += a[k][i]*r2[k]*a[k][j];
-                }
-                b[i][j] = belem;
-                g[i][j] = gelem;
+                b[i][j] = a[0][i]*e2.x*a[0][j] + a[1][i]*e2.y*a[1][j] + a[2][i]*e2.z*a[2][j];
+                g[i][j] = a[0][i]*r2.x*a[0][j] + a[1][i]*r2.y*a[1][j] + a[2][i]*r2.z*a[2][j];
             }
     }
 }
@@ -201,15 +198,15 @@ void loadAtomData(AtomData* data, int sortedIndex, int originalIndex, __global c
 }
 
 real3 matrixVectorProduct(real (*m)[3], real3 v) {
-    return (real3) (m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2],
-                    m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2],
-                    m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2]);
+    return (real3) (m[0][0]*v.x + m[0][1]*v.y + m[0][2]*v.z,
+                    m[1][0]*v.x + m[1][1]*v.y + m[1][2]*v.z,
+                    m[2][0]*v.x + m[2][1]*v.y + m[2][2]*v.z);
 }
 
 real3 vectorMatrixProduct(real3 v, real (*m)[3]) {
-    return (real3) (m[0][0]*v[0] + m[1][0]*v[1] + m[2][0]*v[2],
-                    m[0][1]*v[0] + m[1][1]*v[1] + m[2][1]*v[2],
-                    m[0][2]*v[0] + m[1][2]*v[1] + m[2][2]*v[2]);
+    return (real3) (m[0][0]*v.x + m[1][0]*v.y + m[2][0]*v.z,
+                    m[0][1]*v.x + m[1][1]*v.y + m[2][1]*v.z,
+                    m[0][2]*v.x + m[1][2]*v.y + m[2][2]*v.z);
 }
 
 
@@ -306,31 +303,31 @@ void computeOneInteraction(AtomData* data1, AtomData* data2, real sigma, real ep
         real3 dchidq = cross(vectorMatrixProduct(iota, b), iota)*(-4*rInv2);
         real3 scale = (real3) (sig.y, sig.z, sig.w)*(-0.5f*eta/detG12);
         real d[3][3];
-        d[0][0] = scale[0]*(2*a[0][0]*(G12[1][1]*G12[2][2] - G12[1][2]*G12[2][1]) +
+        d[0][0] = scale.x*(2*a[0][0]*(G12[1][1]*G12[2][2] - G12[1][2]*G12[2][1]) +
                              a[0][2]*(G12[1][2]*G12[0][1] + G12[1][0]*G12[2][1] - G12[1][1]*(G12[0][2] + G12[2][0])) +
                              a[0][1]*(G12[0][2]*G12[2][1] + G12[2][0]*G12[1][2] - G12[2][2]*(G12[0][1] + G12[1][0])));
-        d[0][1] = scale[0]*(  a[0][0]*(G12[0][2]*G12[2][1] + G12[2][0]*G12[1][2] - G12[2][2]*(G12[0][1] + G12[1][0])) +
+        d[0][1] = scale.x*(  a[0][0]*(G12[0][2]*G12[2][1] + G12[2][0]*G12[1][2] - G12[2][2]*(G12[0][1] + G12[1][0])) +
                            2*a[0][1]*(G12[0][0]*G12[2][2] - G12[2][0]*G12[0][2]) +
                              a[0][2]*(G12[1][0]*G12[0][2] + G12[2][0]*G12[0][1] - G12[0][0]*(G12[1][2] + G12[2][1])));
-        d[0][2] = scale[0]*(  a[0][0]*(G12[0][1]*G12[1][2] + G12[1][0]*G12[2][1] - G12[1][1]*(G12[0][2] + G12[2][0])) +
+        d[0][2] = scale.x*(  a[0][0]*(G12[0][1]*G12[1][2] + G12[1][0]*G12[2][1] - G12[1][1]*(G12[0][2] + G12[2][0])) +
                              a[0][1]*(G12[1][0]*G12[0][2] + G12[2][0]*G12[0][1] - G12[0][0]*(G12[1][2] + G12[2][1])) +
                            2*a[0][2]*(G12[1][1]*G12[0][0] - G12[1][0]*G12[0][1]));
-        d[1][0] = scale[1]*(2*a[1][0]*(G12[1][1]*G12[2][2] - G12[1][2]*G12[2][1]) +
+        d[1][0] = scale.y*(2*a[1][0]*(G12[1][1]*G12[2][2] - G12[1][2]*G12[2][1]) +
                              a[1][1]*(G12[0][2]*G12[2][1] + G12[2][0]*G12[1][2] - G12[2][2]*(G12[0][1] + G12[1][0])) +
                              a[1][2]*(G12[1][2]*G12[0][1] + G12[1][0]*G12[2][1] - G12[1][1]*(G12[0][2] + G12[2][0])));
-        d[1][1] = scale[1]*(  a[1][0]*(G12[0][2]*G12[2][1] + G12[2][0]*G12[1][2] - G12[2][2]*(G12[0][1] + G12[1][0])) +
+        d[1][1] = scale.y*(  a[1][0]*(G12[0][2]*G12[2][1] + G12[2][0]*G12[1][2] - G12[2][2]*(G12[0][1] + G12[1][0])) +
                            2*a[1][1]*(G12[2][2]*G12[0][0] - G12[2][0]*G12[0][2]) +
                              a[1][2]*(G12[1][0]*G12[0][2] + G12[0][1]*G12[2][0] - G12[0][0]*(G12[1][2] + G12[2][1])));
-        d[1][2] = scale[1]*(  a[1][0]*(G12[0][1]*G12[1][2] + G12[1][0]*G12[2][1] - G12[1][1]*(G12[0][2] + G12[2][0])) +
+        d[1][2] = scale.y*(  a[1][0]*(G12[0][1]*G12[1][2] + G12[1][0]*G12[2][1] - G12[1][1]*(G12[0][2] + G12[2][0])) +
                              a[1][1]*(G12[1][0]*G12[0][2] + G12[0][1]*G12[2][0] - G12[0][0]*(G12[1][2] + G12[2][1])) +
                            2*a[1][2]*(G12[1][1]*G12[0][0] - G12[1][0]*G12[0][1]));
-        d[2][0] = scale[2]*(2*a[2][0]*(G12[1][1]*G12[2][2] - G12[2][1]*G12[1][2]) +
+        d[2][0] = scale.z*(2*a[2][0]*(G12[1][1]*G12[2][2] - G12[2][1]*G12[1][2]) +
                              a[2][1]*(G12[0][2]*G12[2][1] + G12[1][2]*G12[2][0] - G12[2][2]*(G12[0][1] + G12[1][0])) +
                              a[2][2]*(G12[0][1]*G12[1][2] + G12[2][1]*G12[1][0] - G12[1][1]*(G12[0][2] + G12[2][0])));
-        d[2][1] = scale[2]*(  a[2][0]*(G12[0][2]*G12[2][1] + G12[1][2]*G12[2][0] - G12[2][2]*(G12[0][1] + G12[1][0])) +
+        d[2][1] = scale.z*(  a[2][0]*(G12[0][2]*G12[2][1] + G12[1][2]*G12[2][0] - G12[2][2]*(G12[0][1] + G12[1][0])) +
                            2*a[2][1]*(G12[0][0]*G12[2][2] - G12[0][2]*G12[2][0]) +
                              a[2][2]*(G12[1][0]*G12[0][2] + G12[0][1]*G12[2][0] - G12[0][0]*(G12[1][2] + G12[2][1])));
-        d[2][2] = scale[2]*(  a[2][0]*(G12[0][1]*G12[1][2] + G12[2][1]*G12[1][0] - G12[1][1]*(G12[0][2] + G12[2][0])) +
+        d[2][2] = scale.z*(  a[2][0]*(G12[0][1]*G12[1][2] + G12[2][1]*G12[1][0] - G12[1][1]*(G12[0][2] + G12[2][0])) +
                              a[2][1]*(G12[1][0]*G12[0][2] + G12[2][0]*G12[0][1] - G12[0][0]*(G12[1][2] + G12[2][1])) +
                            2*a[2][2]*(G12[1][1]*G12[0][0] - G12[1][0]*G12[0][1]));
         real3 detadq = 0;
@@ -365,6 +362,8 @@ __kernel void computeForce(
     mixed energy = 0;
 #ifdef USE_CUTOFF
     const int numBlocks = *neighborBlockCount;
+    if (numBlocks > maxNeighborBlocks)
+        return; // There wasn't enough memory for the neighbor list.
     for (int block = get_global_id(0); block < numBlocks; block += get_global_size(0)) {
         // Load parameters for atom1.
         

tests/TestGayBerneForce.h

@@ -133,6 +133,20 @@ void testEnergyScales() {
     state = context.getState(State::Forces | State::Energy);
     ASSERT_EQUAL_TOL(expectedEnergy*expectedScale, state.getPotentialEnergy(), 1e-5);
     ASSERT_EQUAL_VEC(Vec3(0, expectedForce*expectedScale, 0), state.getForces()[3], 1e-5);
+    
+    // Modify their parameters and see if the result is still correct.
+    
+    double newSigma = 1.1*sigma;
+    gb->setParticleParameters(0, newSigma, 1.5*epsilon, 1, 2, newSigma, newSigma, newSigma, 1.2, 1.6, 1.9);
+    gb->setParticleParameters(3, newSigma, epsilon, 4, 5, newSigma, newSigma, newSigma, 1.3, 1.7, 1.8);
+    gb->updateParametersInContext(context);
+    double combinedEpsilon = sqrt(1.5)*epsilon;
+    expectedEnergy = 4*combinedEpsilon*(pow(newSigma, 12.0)-pow(newSigma, 6.0));
+    expectedForce = 4*combinedEpsilon*(12*pow(newSigma, 12.0)-6*pow(newSigma, 6.0));
+    expectedScale = pow(2.0/(1/sqrt(1.6) + 1/sqrt(1.8)), 2.0);
+    state = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(expectedEnergy*expectedScale, state.getPotentialEnergy(), 1e-5);
+    ASSERT_EQUAL_VEC(Vec3(0, expectedForce*expectedScale, 0), state.getForces()[3], 1e-5);
 }
 
 void testEnergyConservation() {
@@ -218,6 +232,16 @@ void testExceptions() {
     State state = context.getState(State::Forces | State::Energy);
     ASSERT_EQUAL_TOL(expectedEnergy*expectedScale, state.getPotentialEnergy(), 1e-5);
     ASSERT_EQUAL_VEC(Vec3(expectedForce*expectedScale, 0, 0), state.getForces()[3], 1e-5);
+    
+    // Modify the exception and see if the results are still correct.
+
+    gb->setExceptionParameters(0, 0, 3, sigma, 3.1*epsilon);
+    gb->updateParametersInContext(context);
+    expectedEnergy = 3.1*4*epsilon*(pow(sigma, 12.0)-pow(sigma, 6.0));
+    expectedForce = 3.1*4*epsilon*(12*pow(sigma, 12.0)-6*pow(sigma, 6.0));
+    state = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(expectedEnergy*expectedScale, state.getPotentialEnergy(), 1e-5);
+    ASSERT_EQUAL_VEC(Vec3(expectedForce*expectedScale, 0, 0), state.getForces()[3], 1e-5);
 }
 
 void runPlatformTests();