Implemented UniqueCentralParticle mode in CUDA version of CustomManyParticleForce

platforms/cuda/src/CudaKernels.cpp

@@ -4472,6 +4472,7 @@ void CudaCalcCustomManyParticleForceKernel::initialize(const System& system, con
     cu.setAsCurrent();
     int numParticles = force.getNumParticles();
     int particlesPerSet = force.getNumParticlesPerSet();
+    bool centralParticleMode = (force.getPermutationMode() == CustomManyParticleForce::UniqueCentralParticle);
     nonbondedMethod = CalcCustomManyParticleForceKernel::NonbondedMethod(force.getNonbondedMethod());
     forceWorkgroupSize = 128;
     findNeighborsWorkgroupSize = 128;
@@ -4790,35 +4791,62 @@ void CudaCalcCustomManyParticleForceKernel::initialize(const System& system, con
         for (int j = 0; j < (int) params->getBuffers().size(); j++)
             loadData<<params->getBuffers()[j].getType()<<" params"<<(j+1)<<(i+1)<<" = global_params"<<(j+1)<<"[atom"<<(i+1)<<"];\n";
     }
+    if (centralParticleMode) {
+        for (int i = 1; i < particlesPerSet; i++) {
+            if (i > 1)
+                isValidCombination<<" && p"<<(i+1)<<">p"<<i<<" && ";
+            isValidCombination<<"p"<<(i+1)<<"!=p1";
+        }
+    }
+    else {
         for (int i = 2; i < particlesPerSet; i++) {
             if (i > 2)
                 isValidCombination<<" && ";
             isValidCombination<<"a"<<(i+1)<<">a"<<i;
         }
+    }
     atomsForCombination<<"int tempIndex = index;\n";
     for (int i = 1; i < particlesPerSet; i++) {
         if (i > 1)
             numCombinations<<"*";
         numCombinations<<"numNeighbors";
+        if (centralParticleMode)
+            atomsForCombination<<"int a"<<(i+1)<<" = tempIndex%numNeighbors;\n";
+        else
             atomsForCombination<<"int a"<<(i+1)<<" = 1+tempIndex%numNeighbors;\n";
         if (i < particlesPerSet-1)
             atomsForCombination<<"tempIndex /= numNeighbors;\n";
     }
-    if (particlesPerSet > 2)
+    if (particlesPerSet > 2) {
+        if (centralParticleMode)
+            atomsForCombination<<"a2 = (a3%2 == 0 ? a2 : numNeighbors-a2-1);\n";
+        else
             atomsForCombination<<"a2 = (a3%2 == 0 ? a2 : numNeighbors-a2+1);\n";
+    }
     for (int i = 1; i < particlesPerSet; i++) {
-        if (nonbondedMethod == NoCutoff)
+        if (nonbondedMethod == NoCutoff) {
+            if (centralParticleMode)
+                atomsForCombination<<"int p"<<(i+1)<<" = a"<<(i+1)<<";\n";
+            else
                 atomsForCombination<<"int p"<<(i+1)<<" = p1+a"<<(i+1)<<";\n";
+        }
+        else {
+            if (centralParticleMode)
+                atomsForCombination<<"int p"<<(i+1)<<" = neighbors[firstNeighbor+a"<<(i+1)<<"];\n";
             else
                 atomsForCombination<<"int p"<<(i+1)<<" = neighbors[firstNeighbor-1+a"<<(i+1)<<"];\n";
         }
+    }
     if (nonbondedMethod != NoCutoff) {
         for (int i = 1; i < particlesPerSet; i++)
             verifyCutoff<<"real3 pos"<<(i+1)<<" = trim(posq[p"<<(i+1)<<"]);\n";
-        for (int i = 1; i < particlesPerSet; i++)
+        if (!centralParticleMode) {
+            for (int i = 1; i < particlesPerSet; i++) {
                 for (int j = i+1; j < particlesPerSet; j++)
                     verifyCutoff<<"includeInteraction &= (delta(pos"<<(i+1)<<", pos"<<(j+1)<<", periodicBoxSize, invPeriodicBoxSize).w < CUTOFF_SQUARED);\n";
             }
+        }
+    }
     if (force.getNumExclusions() > 0) {
         int startCheckFrom = (nonbondedMethod == NoCutoff ? 0 : 1);
         for (int i = startCheckFrom; i < particlesPerSet; i++)
@@ -4855,6 +4883,8 @@ void CudaCalcCustomManyParticleForceKernel::initialize(const System& system, con
         defines["USE_CUTOFF"] = "1";
     if (nonbondedMethod == CutoffPeriodic)
         defines["USE_PERIODIC"] = "1";
+    if (centralParticleMode)
+        defines["USE_CENTRAL_PARTICLE"] = "1";
     if (hasTypeFilters)
         defines["USE_FILTERS"] = "1";
     if (force.getNumExclusions() > 0)

platforms/cuda/src/kernels/customManyParticle.cu

@@ -97,12 +97,20 @@ extern "C" __global__ void computeInteraction(
     // Loop over particles to be the first one in the set.
     
     for (int p1 = blockIdx.x; p1 < NUM_ATOMS; p1 += gridDim.x) {
+#ifdef USE_CENTRAL_PARTICLE
+        const int a1 = p1;
+#else
         const int a1 = 0;
+#endif
 #ifdef USE_CUTOFF
         int firstNeighbor = neighborStartIndex[p1];
         int numNeighbors = neighborStartIndex[p1+1]-firstNeighbor;
 #else
+  #ifdef USE_CENTRAL_PARTICLE
+        int numNeighbors = NUM_ATOMS;
+  #else
         int numNeighbors = NUM_ATOMS-p1-1;
+  #endif
 #endif
         int numCombinations = NUM_CANDIDATE_COMBINATIONS;
         for (int index = threadIdx.x; index < numCombinations; index += blockDim.x) {
@@ -195,7 +203,12 @@ extern "C" __global__ void findNeighbors(real4 periodicBoxSize, real4 invPeriodi
         // Loop over atom blocks to search for neighbors.  The threads in a warp compare block1 against 32
         // other blocks in parallel.
 
-        for (int block2Base = block1; block2Base < NUM_BLOCKS; block2Base += 32) {
+#ifdef USE_CENTRAL_PARTICLE
+        int startBlock = 0;
+#else
+        int startBlock = block1;
+#endif
+        for (int block2Base = startBlock; block2Base < NUM_BLOCKS; block2Base += 32) {
             int block2 = block2Base+indexInWarp;
             bool includeBlock2 = (block2 < NUM_BLOCKS);
             if (includeBlock2) {
@@ -234,7 +247,11 @@ extern "C" __global__ void findNeighbors(real4 periodicBoxSize, real4 invPeriodi
                     // Decide whether to include this atom pair in the neighbor list.
 
                     real4 atomDelta = delta(pos1, pos2, periodicBoxSize, invPeriodicBoxSize);
+#ifdef USE_CENTRAL_PARTICLE
+                    bool includeAtom = (atom2 != atom1 && atom2 < NUM_ATOMS && atomDelta.w < CUTOFF_SQUARED);
+#else
                     bool includeAtom = (atom2 > atom1 && atom2 < NUM_ATOMS && atomDelta.w < CUTOFF_SQUARED);
+#endif
 #ifdef USE_EXCLUSIONS
                     if (includeAtom)
                         includeAtom &= !isInteractionExcluded(atom1, atom2, exclusions, exclusionStartIndex);

platforms/cuda/tests/TestCudaCustomManyParticleForce.cpp

@@ -30,7 +30,7 @@
  * -------------------------------------------------------------------------- */
 
 /**
- * This tests the CPU implementation of CustomManyParticleForce.
+ * This tests the CUDA implementation of CustomManyParticleForce.
  */
 
 #ifdef WIN32
@@ -122,6 +122,76 @@ void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& p
     ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4);
 }
 
+void validateStillingerWeber(CustomManyParticleForce* force, const vector<Vec3>& positions, const vector<const int*>& expectedSets, double boxSize) {
+    // Create a System and Context.
+    
+    int numParticles = force->getNumParticles();
+    CustomManyParticleForce::NonbondedMethod nonbondedMethod = force->getNonbondedMethod();
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    State state1 = context.getState(State::Forces | State::Energy);
+    double L = context.getParameter("L");
+    double eps = context.getParameter("eps");
+    double a = context.getParameter("a");
+    double gamma = context.getParameter("gamma");
+    double sigma = context.getParameter("sigma");
+    
+    // See if the energy matches the expected value.
+    
+    double expectedEnergy = 0;
+    for (int i = 0; i < (int) expectedSets.size(); i++) {
+        int p1 = expectedSets[i][0];
+        int p2 = expectedSets[i][1];
+        int p3 = expectedSets[i][2];
+        Vec3 d12 = positions[p2]-positions[p1];
+        Vec3 d13 = positions[p3]-positions[p1];
+        Vec3 d23 = positions[p3]-positions[p2];
+        if (nonbondedMethod == CustomManyParticleForce::CutoffPeriodic) {
+            for (int j = 0; j < 3; j++) {
+                d12[j] -= floor(d12[j]/boxSize+0.5f)*boxSize;
+                d13[j] -= floor(d13[j]/boxSize+0.5f)*boxSize;
+                d23[j] -= floor(d23[j]/boxSize+0.5f)*boxSize;
+            }
+        }
+        double r12 = sqrt(d12.dot(d12));
+        double r13 = sqrt(d13.dot(d13));
+        double r23 = sqrt(d23.dot(d23));
+        double ctheta1 = d12.dot(d13)/(r12*r13);
+        double ctheta2 = -d12.dot(d23)/(r12*r23);
+        double ctheta3 = d13.dot(d23)/(r13*r23);
+        expectedEnergy += L*eps*(ctheta1+1.0/3.0)*(ctheta1+1.0/3.0)*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));
+    }
+    ASSERT_EQUAL_TOL(expectedEnergy, state1.getPotentialEnergy(), 1e-5);
+
+    // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
+
+    const vector<Vec3>& forces = state1.getForces();
+    double norm = 0.0;
+    for (int i = 0; i < (int) forces.size(); ++i)
+        norm += forces[i].dot(forces[i]);
+    norm = std::sqrt(norm);
+    const double stepSize = 1e-3;
+    double step = 0.5*stepSize/norm;
+    vector<Vec3> positions2(numParticles), positions3(numParticles);
+    for (int i = 0; i < (int) positions.size(); ++i) {
+        Vec3 p = positions[i];
+        Vec3 f = forces[i];
+        positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+        positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
+    }
+    context.setPositions(positions2);
+    State state2 = context.getState(State::Energy);
+    context.setPositions(positions3);
+    State state3 = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4);
+}
+
 void testNoCutoff() {
     CustomManyParticleForce* force = new CustomManyParticleForce(3,
         "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
@@ -500,6 +570,101 @@ void testLargeSystem() {
         ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
 }
 
+void testCentralParticleModeNoCutoff() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "L*eps*(cos(theta1)+1/3)^2*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));"
+        "r12 = distance(p1,p2); r13 = distance(p1,p3); theta1 = angle(p3,p1,p2)");
+    force->setPermutationMode(CustomManyParticleForce::UniqueCentralParticle);
+    force->addGlobalParameter("L", 23.13);
+    force->addGlobalParameter("eps", 25.894776);
+    force->addGlobalParameter("a", 1.8);
+    force->addGlobalParameter("sigma", 0.23925);
+    force->addGlobalParameter("gamma", 1.2);
+    vector<double> params;
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    vector<Vec3> positions;
+    positions.push_back(Vec3(0, 0, 0));
+    positions.push_back(Vec3(1, 0, 0));
+    positions.push_back(Vec3(0, 1.1, 0.3));
+    positions.push_back(Vec3(0.4, 0, -0.8));
+    int sets[12][3] = {{0,1,2}, {0,1,3}, {0,2,3}, {1,0,2}, {1,0,3}, {1, 2, 3}, {2,0,1}, {2,0,3}, {2, 1, 3}, {3,0,1}, {3,0,2}, {3,1,2}};
+    vector<const int*> expectedSets(&sets[0], &sets[12]);
+    validateStillingerWeber(force, positions, expectedSets, 2.0);
+}
+
+void testCentralParticleModeCutoff() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "L*eps*(cos(theta1)+1/3)^2*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));"
+        "r12 = distance(p1,p2); r13 = distance(p1,p3); theta1 = angle(p3,p1,p2)");
+    force->setPermutationMode(CustomManyParticleForce::UniqueCentralParticle);
+    force->addGlobalParameter("L", 23.13);
+    force->addGlobalParameter("eps", 25.894776);
+    force->addGlobalParameter("a", 1.8);
+    force->addGlobalParameter("sigma", 0.23925);
+    force->addGlobalParameter("gamma", 1.2);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffNonPeriodic);
+    force->setCutoffDistance(1.55);
+    vector<double> params;
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    vector<Vec3> positions;
+    positions.push_back(Vec3(0, 0, 0));
+    positions.push_back(Vec3(1, 0, 0));
+    positions.push_back(Vec3(0, 1.1, 0.3));
+    positions.push_back(Vec3(0.4, 0, -0.8));
+    int sets[8][3] = {{0,1,2}, {0,1,3}, {0,2,3}, {1,0,2}, {1,0,3}, {1, 2, 3}, {2,0,1}, {3,0,1}};
+    vector<const int*> expectedSets(&sets[0], &sets[8]);
+    validateStillingerWeber(force, positions, expectedSets, 2.0);
+}
+
+void testCentralParticleModeLargeSystem() {
+    int gridSize = 8;
+    int numParticles = gridSize*gridSize*gridSize;
+    double boxSize = 2.0;
+    double spacing = boxSize/gridSize;
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "L*eps*(cos(theta1)+1/3)^2*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));"
+        "r12 = distance(p1,p2); r13 = distance(p1,p3); theta1 = angle(p3,p1,p2)");
+    force->setPermutationMode(CustomManyParticleForce::UniqueCentralParticle);
+    force->addGlobalParameter("L", 23.13);
+    force->addGlobalParameter("eps", 25.894776);
+    force->addGlobalParameter("a", 1.8);
+    force->addGlobalParameter("sigma", 0.23925);
+    force->addGlobalParameter("gamma", 1.2);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
+    force->setCutoffDistance(1.8*0.23925);
+    vector<double> params;
+    vector<Vec3> positions;
+    System system;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < gridSize; i++)
+        for (int j = 0; j < gridSize; j++)
+            for (int k = 0; k < gridSize; k++) {
+                force->addParticle(params);
+                positions.push_back(Vec3((i+0.4*genrand_real2(sfmt))*spacing, (j+0.4*genrand_real2(sfmt))*spacing, (k+0.4*genrand_real2(sfmt))*spacing));
+                system.addParticle(1.0);
+            }
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(force);
+    VerletIntegrator integrator1(0.001);
+    VerletIntegrator integrator2(0.001);
+    Context context1(system, integrator1, Platform::getPlatformByName("Reference"));
+    Context context2(system, integrator2, platform);
+    context1.setPositions(positions);
+    context2.setPositions(positions);
+    State state1 = context1.getState(State::Forces | State::Energy);
+    State state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
@@ -513,6 +678,9 @@ int main(int argc, char* argv[]) {
         testTabulatedFunctions();
         testTypeFilters();
         testLargeSystem();
+        testCentralParticleModeNoCutoff();
+        testCentralParticleModeCutoff();
+        testCentralParticleModeLargeSystem();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;