Merge pull request #604 from peastman/many

Created CustomManyParticleForce

platforms/cpu/tests/TestCpuCustomManyParticleForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2014 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the CPU implementation of CustomManyParticleForce.
+ */
+
+#ifdef WIN32
+  #define _USE_MATH_DEFINES // Needed to get M_PI
+#endif
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CpuPlatform.h"
+#include "openmm/CustomCompoundBondForce.h"
+#include "openmm/CustomManyParticleForce.h"
+#include "openmm/System.h"
+#include "openmm/TabulatedFunction.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void validateAxilrodTeller(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);
+    CpuPlatform platform;
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    State state1 = context.getState(State::Forces | State::Energy);
+    double c = context.getParameter("C");
+    
+    // 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);
+        double rprod = r12*r13*r23;
+        expectedEnergy += c*(1+3*ctheta1*ctheta2*ctheta3)/(rprod*rprod*rprod);
+    }
+    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 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);
+    CpuPlatform platform;
+    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;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    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[4][3] = {{0,1,2}, {1,2,3}, {2,3,0}, {3,0,1}};
+    vector<const int*> expectedSets(&sets[0], &sets[4]);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+}
+
+void testCutoff() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffNonPeriodic);
+    force->setCutoffDistance(1.55);
+    vector<double> params;
+    force->addParticle(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));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    int sets[7][3] = {{0,1,2}, {0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,2,4}, {1,3,4}};
+    vector<const int*> expectedSets(&sets[0], &sets[7]);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+}
+
+void testPeriodic() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
+    force->setCutoffDistance(1.05);
+    vector<double> params;
+    force->addParticle(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));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    double boxSize = 2.1;
+    int sets[5][3] = {{0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,3,4}};
+    vector<const int*> expectedSets(&sets[0], &sets[5]);
+    validateAxilrodTeller(force, positions, expectedSets, boxSize);
+}
+
+void testExclusions() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    vector<double> params;
+    force->addParticle(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));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    force->addExclusion(0, 2);
+    force->addExclusion(0, 3);
+    int sets[5][3] = {{0,1,4}, {1,2,3}, {1,2,4}, {1,3,4}, {2,3,4}};
+    vector<const int*> expectedSets(&sets[0], &sets[5]);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+}
+
+void testAllTerms() {
+    int numParticles = 4;
+    CpuPlatform platform;
+    
+    // Create a system with a CustomManyParticleForce.
+    
+    System system1;
+    CustomManyParticleForce* force1 = new CustomManyParticleForce(4,
+        "distance(p1,p2)+angle(p1,p4,p3)+dihedral(p1,p3,p2,p4)+x1+y4+z3");
+    system1.addForce(force1);
+    vector<double> params;
+    for (int i = 0; i < numParticles; i++) {
+        system1.addParticle(1.0);
+        force1->addParticle(params, i);
+    }
+    set<int> filter;
+    filter.insert(0);
+    force1->setTypeFilter(0, filter);
+    filter.clear();
+    filter.insert(1);
+    force1->setTypeFilter(1, filter);
+    filter.clear();
+    filter.insert(3);
+    force1->setTypeFilter(2, filter);
+    filter.clear();
+    filter.insert(2);
+    force1->setTypeFilter(3, filter);
+    
+    // Create a system that use a CustomCompoundBondForce to compute exactly the same interactions.
+    
+    System system2;
+    CustomCompoundBondForce* force2 = new CustomCompoundBondForce(4,
+        "distance(p1,p2)+angle(p1,p3,p4)+dihedral(p1,p4,p2,p3)+x1+y3+z4");
+    system2.addForce(force2);
+    vector<int> particles;
+    particles.push_back(0);
+    particles.push_back(1);
+    particles.push_back(2);
+    particles.push_back(3);
+    force2->addBond(particles, params);
+    for (int i = 0; i < numParticles; i++)
+        system2.addParticle(1.0);
+    
+    // Create contexts for both of them.
+
+    vector<Vec3> positions;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; i++)
+        positions.push_back(Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
+    VerletIntegrator integrator1(0.001);
+    VerletIntegrator integrator2(0.001);
+    Context context1(system1, integrator1, platform);
+    Context context2(system2, integrator2, platform);
+    context1.setPositions(positions);
+    context2.setPositions(positions);
+    
+    // See if they produce identical forces and energies.
+    
+    State state1 = context1.getState(State::Forces | State::Energy);
+    State state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state2.getPotentialEnergy(), state1.getPotentialEnergy(), 1e-4);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state2.getForces()[i], state1.getForces()[i], 1e-4);
+}
+
+void testParameters() {
+    // Create a system.
+    
+    int numParticles = 5;
+    System system;
+    CustomManyParticleForce* force = new CustomManyParticleForce(3, "C*scale1*scale2*scale3*(distance(p1,p2)+distance(p2,p3)+distance(p1,p3))");
+    force->addGlobalParameter("C", 2.0);
+    force->addPerParticleParameter("scale");
+    vector<double> params(1);
+    vector<Vec3> positions;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; i++) {
+        params[0] = i+1;
+        force->addParticle(params);
+        positions.push_back(Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
+        system.addParticle(1.0);
+    }
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    CpuPlatform platform;
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    
+    // See if the energy is correct.
+
+    State state = context.getState(State::Energy);
+    double expectedEnergy = 0;
+    for (int i = 0; i < numParticles; i++)
+        for (int j = i+1; j < numParticles; j++)
+            for (int k = j+1; k < numParticles; k++) {
+                Vec3 d12 = positions[j]-positions[i];
+                Vec3 d13 = positions[k]-positions[i];
+                Vec3 d23 = positions[k]-positions[j];
+                double r12 = sqrt(d12.dot(d12));
+                double r13 = sqrt(d13.dot(d13));
+                double r23 = sqrt(d23.dot(d23));
+                expectedEnergy += 2.0*(i+1)*(j+1)*(k+1)*(r12+r13+r23);
+            }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+    
+    // Modify the parameters.
+    
+    context.setParameter("C", 3.5);
+    for (int i = 0; i < numParticles; i++) {
+        params[0] = 0.5*i-0.1;
+        force->setParticleParameters(i, params, 0);
+    }
+    force->updateParametersInContext(context);
+    
+    // See if the energy is still correct.
+    
+    state = context.getState(State::Energy);
+    expectedEnergy = 0;
+    for (int i = 0; i < numParticles; i++)
+        for (int j = i+1; j < numParticles; j++)
+            for (int k = j+1; k < numParticles; k++) {
+                Vec3 d12 = positions[j]-positions[i];
+                Vec3 d13 = positions[k]-positions[i];
+                Vec3 d23 = positions[k]-positions[j];
+                double r12 = sqrt(d12.dot(d12));
+                double r13 = sqrt(d13.dot(d13));
+                double r23 = sqrt(d23.dot(d23));
+                expectedEnergy += 3.5*(0.5*i-0.1)*(0.5*j-0.1)*(0.5*k-0.1)*(r12+r13+r23);
+            }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+}
+
+void testTabulatedFunctions() {
+    int numParticles = 5;
+    
+    // Create two tabulated functions.
+    
+    vector<double> values;
+    values.push_back(0.0);
+    values.push_back(50.0);
+    Continuous1DFunction* f1 = new Continuous1DFunction(values, 0, 100);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<double> c(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        c[i] = genrand_real2(sfmt);
+    values.resize(numParticles*numParticles*numParticles);
+    for (int i = 0; i < numParticles; i++)
+        for (int j = 0; j < numParticles; j++)
+            for (int k = 0; k < numParticles; k++)
+                values[i+numParticles*j+numParticles*numParticles*k] = c[i]+c[j]+c[k];
+    Discrete3DFunction* f2 = new Discrete3DFunction(numParticles, numParticles, numParticles, values);
+    
+    // Create a system.
+    
+    System system;
+    CustomManyParticleForce* force = new CustomManyParticleForce(3, "f1(distance(p1,p2)+distance(p2,p3)+distance(p1,p3))*f2(atom1, atom2, atom3)");
+    force->addPerParticleParameter("atom");
+    force->addTabulatedFunction("f1", f1);
+    force->addTabulatedFunction("f2", f2);
+    vector<double> params(1);
+    vector<Vec3> positions;
+    for (int i = 0; i < numParticles; i++) {
+        params[0] = i;
+        force->addParticle(params);
+        positions.push_back(Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
+        system.addParticle(1.0);
+    }
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    CpuPlatform platform;
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    
+    // See if the energy is correct.
+
+    State state = context.getState(State::Energy);
+    double expectedEnergy = 0;
+    for (int i = 0; i < numParticles; i++)
+        for (int j = i+1; j < numParticles; j++)
+            for (int k = j+1; k < numParticles; k++) {
+                Vec3 d12 = positions[j]-positions[i];
+                Vec3 d13 = positions[k]-positions[i];
+                Vec3 d23 = positions[k]-positions[j];
+                double r12 = sqrt(d12.dot(d12));
+                double r13 = sqrt(d13.dot(d13));
+                double r23 = sqrt(d23.dot(d23));
+                expectedEnergy += 0.5*(r12+r13+r23)*(c[i]+c[j]+c[k]);
+            }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+}
+
+void testTypeFilters() {
+    // Create a system.
+    
+    System system;
+    for (int i = 0; i < 5; i++)
+        system.addParticle(1.0);
+    CustomManyParticleForce* force = new CustomManyParticleForce(3, "c1*(distance(p1,p2)+distance(p1,p3))");
+    force->addPerParticleParameter("c");
+    double c[] = {1.0, 2.0, 1.3, 1.5, -2.1};
+    int type[] = {0, 1, 0, 1, 5};
+    vector<double> params(1);
+    for (int i = 0; i < 5; i++) {
+        params[0] = c[i];
+        force->addParticle(params, type[i]);
+    }
+    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));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    set<int> f1, f2;
+    f1.insert(0);
+    f2.insert(1);
+    f2.insert(5);
+    force->setTypeFilter(0, f1);
+    force->setTypeFilter(1, f2);
+    force->setTypeFilter(2, f2);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    CpuPlatform platform;
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    
+    // See if the energy is correct.
+
+    State state = context.getState(State::Energy);
+    double expectedEnergy = 0;
+    int sets[6][3] = {{0,1,3}, {0,1,4}, {0,3,4}, {2,1,3}, {2,1,4}, {2,3,4}};
+    for (int i = 0; i < 6; i++) {
+        int p1 = sets[i][0];
+        int p2 = sets[i][1];
+        int p3 = sets[i][2];
+            Vec3 d12 = positions[p2]-positions[p1];
+            Vec3 d13 = positions[p3]-positions[p1];
+            double r12 = sqrt(d12.dot(d12));
+            double r13 = sqrt(d13.dot(d13));
+            expectedEnergy += c[p1]*(r12+r13);
+    }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+}
+
+void testLargeSystem() {
+    int gridSize = 8;
+    int numParticles = gridSize*gridSize*gridSize;
+    double boxSize = 3.0;
+    double spacing = boxSize/gridSize;
+    CpuPlatform platform;
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
+    force->setCutoffDistance(0.6);
+    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);
+}
+
+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(0.1, 0, 0));
+    positions.push_back(Vec3(0, 0.11, 0.03));
+    positions.push_back(Vec3(0.04, 0, -0.08));
+    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(0.155);
+    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(0.1, 0, 0));
+    positions.push_back(Vec3(0, 0.11, 0.03));
+    positions.push_back(Vec3(0.04, 0, -0.08));
+    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;
+    CpuPlatform platform;
+    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() {
+    try {
+        testNoCutoff();
+        testCutoff();
+        testPeriodic();
+        testExclusions();
+        testAllTerms();
+        testParameters();
+        testTabulatedFunctions();
+        testTypeFilters();
+        testLargeSystem();
+        testCentralParticleModeNoCutoff();
+        testCentralParticleModeCutoff();
+        testCentralParticleModeLargeSystem();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/cuda/include/CudaKernels.h

@@ -925,6 +925,69 @@ private:
     const System& system;
 };
 
+/**
+ * This kernel is invoked by CustomManyParticleForce to calculate the forces acting on the system.
+ */
+class CudaCalcCustomManyParticleForceKernel : public CalcCustomManyParticleForceKernel {
+public:
+    CudaCalcCustomManyParticleForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system) : CalcCustomManyParticleForceKernel(name, platform),
+            hasInitializedKernel(false), cu(cu), params(NULL), particleTypes(NULL), orderIndex(NULL), particleOrder(NULL), exclusions(NULL),
+            exclusionStartIndex(NULL), blockCenter(NULL), blockBoundingBox(NULL), neighborPairs(NULL), numNeighborPairs(NULL), neighborStartIndex(NULL),
+            numNeighborsForAtom(NULL), neighbors(NULL), system(system) {
+    }
+    ~CudaCalcCustomManyParticleForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomManyParticleForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomManyParticleForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomManyParticleForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force);
+
+private:
+    CudaContext& cu;
+    bool hasInitializedKernel;
+    NonbondedMethod nonbondedMethod;
+    int maxNeighborPairs, forceWorkgroupSize, findNeighborsWorkgroupSize;
+    CudaParameterSet* params;
+    CudaArray* particleTypes;
+    CudaArray* orderIndex;
+    CudaArray* particleOrder;
+    CudaArray* exclusions;
+    CudaArray* exclusionStartIndex;
+    CudaArray* blockCenter;
+    CudaArray* blockBoundingBox;
+    CudaArray* neighborPairs;
+    CudaArray* numNeighborPairs;
+    CudaArray* neighborStartIndex;
+    CudaArray* numNeighborsForAtom;
+    CudaArray* neighbors;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+    std::vector<CudaArray*> tabulatedFunctions;
+    std::vector<void*> forceArgs, blockBoundsArgs, neighborsArgs, startIndicesArgs, copyPairsArgs;
+    const System& system;
+    CUfunction forceKernel, blockBoundsKernel, neighborsKernel, startIndicesKernel, copyPairsKernel;
+    CUdeviceptr globalsPtr;
+    CUevent event;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

platforms/cuda/src/CudaKernelFactory.cpp

@@ -106,6 +106,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaCalcCustomHbondForceKernel(name, platform, cu, context.getSystem());
     if (name == CalcCustomCompoundBondForceKernel::Name())
         return new CudaCalcCustomCompoundBondForceKernel(name, platform, cu, context.getSystem());
+    if (name == CalcCustomManyParticleForceKernel::Name())
+        return new CudaCalcCustomManyParticleForceKernel(name, platform, cu, context.getSystem());
     if (name == IntegrateVerletStepKernel::Name())
         return new CudaIntegrateVerletStepKernel(name, platform, cu);
     if (name == IntegrateLangevinStepKernel::Name())

platforms/cuda/src/CudaKernels.cpp

@@ -33,6 +33,7 @@
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/CustomCompoundBondForceImpl.h"
 #include "openmm/internal/CustomHbondForceImpl.h"
+#include "openmm/internal/CustomManyParticleForceImpl.h"
 #include "openmm/internal/CustomNonbondedForceImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
 #include "CudaBondedUtilities.h"
@@ -55,6 +56,13 @@ using namespace std;
 using Lepton::ExpressionTreeNode;
 using Lepton::Operation;
 
+#define CHECK_RESULT(result, prefix) \
+    if (result != CUDA_SUCCESS) { \
+        std::stringstream m; \
+        m<<prefix<<": "<<CudaContext::getErrorString(result)<<" ("<<result<<")"<<" at "<<__FILE__<<":"<<__LINE__; \
+        throw OpenMMException(m.str());\
+    }
+
 static bool isZeroExpression(const Lepton::ParsedExpression& expression) {
     const Lepton::Operation& op = expression.getRootNode().getOperation();
     if (op.getId() != Lepton::Operation::CONSTANT)
@@ -1962,7 +1970,7 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const
     map<string, Lepton::CustomFunction*> functions;
     vector<pair<string, string> > functionDefinitions;
     vector<const TabulatedFunction*> functionList;
-    for (int i = 0; i < force.getNumFunctions(); i++) {
+    for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
         functionList.push_back(&force.getTabulatedFunction(i));
         string name = force.getTabulatedFunctionName(i);
         string arrayName = prefix+"table"+cu.intToString(i);
@@ -3765,7 +3773,7 @@ void CudaCalcCustomHbondForceKernel::initialize(const System& system, const Cust
     vector<pair<string, string> > functionDefinitions;
     vector<const TabulatedFunction*> functionList;
     stringstream tableArgs;
-    for (int i = 0; i < force.getNumFunctions(); i++) {
+    for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
         functionList.push_back(&force.getTabulatedFunction(i));
         string name = force.getTabulatedFunctionName(i);
         string arrayName = "table"+cu.intToString(i);
@@ -4146,7 +4154,7 @@ void CudaCalcCustomCompoundBondForceKernel::initialize(const System& system, con
     vector<pair<string, string> > functionDefinitions;
     vector<const TabulatedFunction*> functionList;
     stringstream tableArgs;
-    for (int i = 0; i < force.getNumFunctions(); i++) {
+    for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
         functionList.push_back(&force.getTabulatedFunction(i));
         string name = force.getTabulatedFunctionName(i);
         functions[name] = cu.getExpressionUtilities().getFunctionPlaceholder(force.getTabulatedFunction(i));
@@ -4394,6 +4402,660 @@ void CudaCalcCustomCompoundBondForceKernel::copyParametersToContext(ContextImpl&
     cu.invalidateMolecules();
 }
 
+class CudaCustomManyParticleForceInfo : public CudaForceInfo {
+public:
+    CudaCustomManyParticleForceInfo(const CustomManyParticleForce& force) : force(force) {
+    }
+    bool areParticlesIdentical(int particle1, int particle2) {
+        vector<double> params1, params2;
+        int type1, type2;
+        force.getParticleParameters(particle1, params1, type1);
+        force.getParticleParameters(particle2, params2, type2);
+        if (type1 != type2)
+            return false;
+        for (int i = 0; i < (int) params1.size(); i++)
+            if (params1[i] != params2[i])
+                return false;
+        return true;
+    }
+    int getNumParticleGroups() {
+        return force.getNumExclusions();
+    }
+    void getParticlesInGroup(int index, vector<int>& particles) {
+        int particle1, particle2;
+        force.getExclusionParticles(index, particle1, particle2);
+        particles.resize(2);
+        particles[0] = particle1;
+        particles[1] = particle2;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        return true;
+    }
+private:
+    const CustomManyParticleForce& force;
+};
+
+CudaCalcCustomManyParticleForceKernel::~CudaCalcCustomManyParticleForceKernel() {
+    cu.setAsCurrent();
+    if (params != NULL)
+        delete params;
+    if (orderIndex != NULL)
+        delete orderIndex;
+    if (particleOrder != NULL)
+        delete particleOrder;
+    if (particleTypes != NULL)
+        delete particleTypes;
+    if (exclusions != NULL)
+        delete exclusions;
+    if (exclusionStartIndex != NULL)
+        delete exclusionStartIndex;
+    if (blockCenter != NULL)
+        delete blockCenter;
+    if (blockBoundingBox != NULL)
+        delete blockBoundingBox;
+    if (neighborPairs != NULL)
+        delete neighborPairs;
+    if (numNeighborPairs != NULL)
+        delete numNeighborPairs;
+    if (neighborStartIndex != NULL)
+        delete neighborStartIndex;
+    if (neighbors != NULL)
+        delete neighbors;
+    if (numNeighborsForAtom != NULL)
+        delete numNeighborsForAtom;
+    for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
+        delete tabulatedFunctions[i];
+}
+
+void CudaCalcCustomManyParticleForceKernel::initialize(const System& system, const CustomManyParticleForce& force) {
+    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;
+    
+    // Record parameter values.
+    
+    params = new CudaParameterSet(cu, force.getNumPerParticleParameters(), numParticles, "customManyParticleParameters");
+    vector<vector<float> > paramVector(numParticles);
+    for (int i = 0; i < numParticles; i++) {
+        vector<double> parameters;
+        int type;
+        force.getParticleParameters(i, parameters, type);
+        paramVector[i].resize(parameters.size());
+        for (int j = 0; j < (int) parameters.size(); j++)
+            paramVector[i][j] = (float) parameters[j];
+    }
+    params->setParameterValues(paramVector);
+    cu.addForce(new CudaCustomManyParticleForceInfo(force));
+
+    // Record the tabulated functions.
+
+    map<string, Lepton::CustomFunction*> functions;
+    vector<pair<string, string> > functionDefinitions;
+    vector<const TabulatedFunction*> functionList;
+    stringstream tableArgs;
+    for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
+        functionList.push_back(&force.getTabulatedFunction(i));
+        string name = force.getTabulatedFunctionName(i);
+        string arrayName = "table"+cu.intToString(i);
+        functionDefinitions.push_back(make_pair(name, arrayName));
+        functions[name] = cu.getExpressionUtilities().getFunctionPlaceholder(force.getTabulatedFunction(i));
+        int width;
+        vector<float> f = cu.getExpressionUtilities().computeFunctionCoefficients(force.getTabulatedFunction(i), width);
+        tabulatedFunctions.push_back(CudaArray::create<float>(cu, f.size(), "TabulatedFunction"));
+        tabulatedFunctions[tabulatedFunctions.size()-1]->upload(f);
+        tableArgs << ", const float";
+        if (width > 1)
+            tableArgs << width;
+        tableArgs << "* __restrict__ " << arrayName;
+    }
+    
+    // Record information about parameters.
+
+    globalParamNames.resize(force.getNumGlobalParameters());
+    globalParamValues.resize(force.getNumGlobalParameters());
+    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+        globalParamNames[i] = force.getGlobalParameterName(i);
+        globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
+    }
+    vector<pair<ExpressionTreeNode, string> > variables;
+    for (int i = 0; i < particlesPerSet; i++) {
+        string index = cu.intToString(i+1);
+        variables.push_back(makeVariable("x"+index, "pos"+index+".x"));
+        variables.push_back(makeVariable("y"+index, "pos"+index+".y"));
+        variables.push_back(makeVariable("z"+index, "pos"+index+".z"));
+    }
+    for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
+        const string& name = force.getPerParticleParameterName(i);
+        for (int j = 0; j < particlesPerSet; j++) {
+            string index = cu.intToString(j+1);
+            variables.push_back(makeVariable(name+index, "params"+params->getParameterSuffix(i, index)));
+        }
+    }
+    if (force.getNumGlobalParameters() > 0) {
+        for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+            const string& name = force.getGlobalParameterName(i);
+            string value = "globals["+cu.intToString(i)+"]";
+            variables.push_back(makeVariable(name, value));
+        }
+    }
+    
+    // Build data structures for type filters.
+    
+    vector<int> particleTypesVec;
+    vector<int> orderIndexVec;
+    vector<std::vector<int> > particleOrderVec;
+    int numTypes;
+    CustomManyParticleForceImpl::buildFilterArrays(force, numTypes, particleTypesVec, orderIndexVec, particleOrderVec);
+    bool hasTypeFilters = (particleOrderVec.size() > 1);
+    if (hasTypeFilters) {
+        particleTypes = CudaArray::create<int>(cu, particleTypesVec.size(), "customManyParticleTypes");
+        orderIndex = CudaArray::create<int>(cu, orderIndexVec.size(), "customManyParticleOrderIndex");
+        particleOrder = CudaArray::create<int>(cu, particleOrderVec.size()*particlesPerSet, "customManyParticleOrder");
+        particleTypes->upload(particleTypesVec);
+        orderIndex->upload(orderIndexVec);
+        vector<int> flattenedOrder(particleOrder->getSize());
+        for (int i = 0; i < (int) particleOrderVec.size(); i++)
+            for (int j = 0; j < particlesPerSet; j++)
+                flattenedOrder[i*particlesPerSet+j] = particleOrderVec[i][j];
+        particleOrder->upload(flattenedOrder);
+    }
+    
+    // Build data structures for exclusions.
+    
+    if (force.getNumExclusions() > 0) {
+        vector<vector<int> > particleExclusions(numParticles);
+        for (int i = 0; i < force.getNumExclusions(); i++) {
+            int p1, p2;
+            force.getExclusionParticles(i, p1, p2);
+            particleExclusions[p1].push_back(p2);
+            particleExclusions[p2].push_back(p1);
+        }
+        vector<int> exclusionsVec;
+        vector<int> exclusionStartIndexVec(numParticles+1);
+        exclusionStartIndexVec[0] = 0;
+        for (int i = 0; i < numParticles; i++) {
+            sort(particleExclusions[i].begin(), particleExclusions[i].end());
+            exclusionsVec.insert(exclusionsVec.end(), particleExclusions[i].begin(), particleExclusions[i].end());
+            exclusionStartIndexVec[i+1] = exclusionsVec.size();
+        }
+        exclusions = CudaArray::create<int>(cu, exclusionsVec.size(), "customManyParticleExclusions");
+        exclusionStartIndex = CudaArray::create<int>(cu, exclusionStartIndexVec.size(), "customManyParticleExclusionStart");
+        exclusions->upload(exclusionsVec);
+        exclusionStartIndex->upload(exclusionStartIndexVec);
+    }
+    
+    // Build data structures for the neighbor list.
+    
+    if (nonbondedMethod != NoCutoff) {
+        int numAtomBlocks = cu.getNumAtomBlocks();
+        int elementSize = (cu.getUseDoublePrecision() ? sizeof(double) : sizeof(float));
+        blockCenter = new CudaArray(cu, numAtomBlocks, 4*elementSize, "blockCenter");
+        blockBoundingBox = new CudaArray(cu, numAtomBlocks, 4*elementSize, "blockBoundingBox");
+        numNeighborPairs = CudaArray::create<int>(cu, 1, "customManyParticleNumNeighborPairs");
+        neighborStartIndex = CudaArray::create<int>(cu, numParticles+1, "customManyParticleNeighborStartIndex");
+        numNeighborsForAtom = CudaArray::create<int>(cu, numParticles, "customManyParticleNumNeighborsForAtom");
+        CHECK_RESULT(cuEventCreate(&event, CU_EVENT_DISABLE_TIMING), "Error creating event for CustomManyParticleForce");
+
+        // Select a size for the array that holds the neighbor list.  We have to make a fairly
+        // arbitrary guess, but if this turns out to be too small we'll increase it later.
+
+        maxNeighborPairs = 150*numParticles;
+        neighborPairs = CudaArray::create<int2>(cu, maxNeighborPairs, "customManyParticleNeighborPairs");
+        neighbors = CudaArray::create<int>(cu, maxNeighborPairs, "customManyParticleNeighbors");
+    }
+
+    // Now to generate the kernel.  First, it needs to calculate all distances, angles,
+    // and dihedrals the expression depends on.
+
+    map<string, vector<int> > distances;
+    map<string, vector<int> > angles;
+    map<string, vector<int> > dihedrals;
+    Lepton::ParsedExpression energyExpression = CustomManyParticleForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
+    map<string, Lepton::ParsedExpression> forceExpressions;
+    set<string> computedDeltas;
+    vector<string> atomNames, posNames;
+    for (int i = 0; i < particlesPerSet; i++) {
+        string index = cu.intToString(i+1);
+        atomNames.push_back("P"+index);
+        posNames.push_back("pos"+index);
+    }
+    stringstream compute;
+    int index = 0;
+    for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName = atomNames[atoms[0]]+atomNames[atoms[1]];
+        if (computedDeltas.count(deltaName) == 0) {
+            compute<<"real4 delta"<<deltaName<<" = delta("<<posNames[atoms[0]]<<", "<<posNames[atoms[1]]<<", periodicBoxSize, invPeriodicBoxSize);\n";
+            computedDeltas.insert(deltaName);
+        }
+        compute<<"real r_"<<deltaName<<" = sqrt(delta"<<deltaName<<".w);\n";
+        variables.push_back(makeVariable(iter->first, "r_"+deltaName));
+        forceExpressions["real dEdDistance"+cu.intToString(index)+" = "] = energyExpression.differentiate(iter->first).optimize();
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName1 = atomNames[atoms[1]]+atomNames[atoms[0]];
+        string deltaName2 = atomNames[atoms[1]]+atomNames[atoms[2]];
+        string angleName = "angle_"+atomNames[atoms[0]]+atomNames[atoms[1]]+atomNames[atoms[2]];
+        if (computedDeltas.count(deltaName1) == 0) {
+            compute<<"real4 delta"<<deltaName1<<" = delta("<<posNames[atoms[1]]<<", "<<posNames[atoms[0]]<<", periodicBoxSize, invPeriodicBoxSize);\n";
+            computedDeltas.insert(deltaName1);
+        }
+        if (computedDeltas.count(deltaName2) == 0) {
+            compute<<"real4 delta"<<deltaName2<<" = delta("<<posNames[atoms[1]]<<", "<<posNames[atoms[2]]<<", periodicBoxSize, invPeriodicBoxSize);\n";
+            computedDeltas.insert(deltaName2);
+        }
+        compute<<"real "<<angleName<<" = computeAngle(delta"<<deltaName1<<", delta"<<deltaName2<<");\n";
+        variables.push_back(makeVariable(iter->first, angleName));
+        forceExpressions["real dEdAngle"+cu.intToString(index)+" = "] = energyExpression.differentiate(iter->first).optimize();
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName1 = atomNames[atoms[0]]+atomNames[atoms[1]];
+        string deltaName2 = atomNames[atoms[2]]+atomNames[atoms[1]];
+        string deltaName3 = atomNames[atoms[2]]+atomNames[atoms[3]];
+        string crossName1 = "cross_"+deltaName1+"_"+deltaName2;
+        string crossName2 = "cross_"+deltaName2+"_"+deltaName3;
+        string dihedralName = "dihedral_"+atomNames[atoms[0]]+atomNames[atoms[1]]+atomNames[atoms[2]]+atomNames[atoms[3]];
+        if (computedDeltas.count(deltaName1) == 0) {
+            compute<<"real4 delta"<<deltaName1<<" = delta("<<posNames[atoms[0]]<<", "<<posNames[atoms[1]]<<", periodicBoxSize, invPeriodicBoxSize);\n";
+            computedDeltas.insert(deltaName1);
+        }
+        if (computedDeltas.count(deltaName2) == 0) {
+            compute<<"real4 delta"<<deltaName2<<" = delta("<<posNames[atoms[2]]<<", "<<posNames[atoms[1]]<<", periodicBoxSize, invPeriodicBoxSize);\n";
+            computedDeltas.insert(deltaName2);
+        }
+        if (computedDeltas.count(deltaName3) == 0) {
+            compute<<"real4 delta"<<deltaName3<<" = delta("<<posNames[atoms[2]]<<", "<<posNames[atoms[3]]<<", periodicBoxSize, invPeriodicBoxSize);\n";
+            computedDeltas.insert(deltaName3);
+        }
+        compute<<"real4 "<<crossName1<<" = computeCross(delta"<<deltaName1<<", delta"<<deltaName2<<");\n";
+        compute<<"real4 "<<crossName2<<" = computeCross(delta"<<deltaName2<<", delta"<<deltaName3<<");\n";
+        compute<<"real "<<dihedralName<<" = computeAngle("<<crossName1<<", "<<crossName2<<");\n";
+        compute<<dihedralName<<" *= (delta"<<deltaName1<<".x*"<<crossName2<<".x + delta"<<deltaName1<<".y*"<<crossName2<<".y + delta"<<deltaName1<<".z*"<<crossName2<<".z < 0 ? -1 : 1);\n";
+        variables.push_back(makeVariable(iter->first, dihedralName));
+        forceExpressions["real dEdDihedral"+cu.intToString(index)+" = "] = energyExpression.differentiate(iter->first).optimize();
+    }
+
+    // Now evaluate the expressions.
+
+    for (int i = 0; i < (int) params->getBuffers().size(); i++) {
+        CudaNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
+        compute<<buffer.getType()<<" params"<<(i+1)<<" = global_params"<<(i+1)<<"[index];\n";
+    }
+    forceExpressions["energy += "] = energyExpression;
+    compute << cu.getExpressionUtilities().createExpressions(forceExpressions, variables, functionList, functionDefinitions, "temp");
+
+    // Apply forces to atoms.
+
+    vector<string> forceNames;
+    for (int i = 0; i < particlesPerSet; i++) {
+        string istr = cu.intToString(i+1);
+        string forceName = "force"+istr;
+        forceNames.push_back(forceName);
+        compute<<"real3 "<<forceName<<" = make_real3(0);\n";
+        compute<<"{\n";
+        Lepton::ParsedExpression forceExpressionX = energyExpression.differentiate("x"+istr).optimize();
+        Lepton::ParsedExpression forceExpressionY = energyExpression.differentiate("y"+istr).optimize();
+        Lepton::ParsedExpression forceExpressionZ = energyExpression.differentiate("z"+istr).optimize();
+        map<string, Lepton::ParsedExpression> expressions;
+        if (!isZeroExpression(forceExpressionX))
+            expressions[forceName+".x -= "] = forceExpressionX;
+        if (!isZeroExpression(forceExpressionY))
+            expressions[forceName+".y -= "] = forceExpressionY;
+        if (!isZeroExpression(forceExpressionZ))
+            expressions[forceName+".z -= "] = forceExpressionZ;
+        if (expressions.size() > 0)
+            compute<<cu.getExpressionUtilities().createExpressions(expressions, variables, functionList, functionDefinitions, "coordtemp");
+        compute<<"}\n";
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName = atomNames[atoms[0]]+atomNames[atoms[1]];
+        string value = "(dEdDistance"+cu.intToString(index)+"/r_"+deltaName+")*trim(delta"+deltaName+")";
+        compute<<forceNames[atoms[0]]<<" += "<<"-"<<value<<";\n";
+        compute<<forceNames[atoms[1]]<<" += "<<value<<";\n";
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName1 = atomNames[atoms[1]]+atomNames[atoms[0]];
+        string deltaName2 = atomNames[atoms[1]]+atomNames[atoms[2]];
+        compute<<"{\n";
+        compute<<"real3 crossProd = cross(delta"<<deltaName2<<", delta"<<deltaName1<<");\n";
+        compute<<"real lengthCross = max(SQRT(dot(crossProd, crossProd)), 1e-6f);\n";
+        compute<<"real3 deltaCross0 = -cross(trim(delta"<<deltaName1<<"), crossProd)*dEdAngle"<<cu.intToString(index)<<"/(delta"<<deltaName1<<".w*lengthCross);\n";
+        compute<<"real3 deltaCross2 = cross(trim(delta"<<deltaName2<<"), crossProd)*dEdAngle"<<cu.intToString(index)<<"/(delta"<<deltaName2<<".w*lengthCross);\n";
+        compute<<"real3 deltaCross1 = -(deltaCross0+deltaCross2);\n";
+        compute<<forceNames[atoms[0]]<<" += deltaCross0;\n";
+        compute<<forceNames[atoms[1]]<<" += deltaCross1;\n";
+        compute<<forceNames[atoms[2]]<<" += deltaCross2;\n";
+        compute<<"}\n";
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName1 = atomNames[atoms[0]]+atomNames[atoms[1]];
+        string deltaName2 = atomNames[atoms[2]]+atomNames[atoms[1]];
+        string deltaName3 = atomNames[atoms[2]]+atomNames[atoms[3]];
+        string crossName1 = "cross_"+deltaName1+"_"+deltaName2;
+        string crossName2 = "cross_"+deltaName2+"_"+deltaName3;
+        compute<<"{\n";
+        compute<<"real r = sqrt(delta"<<deltaName2<<".w);\n";
+        compute<<"real4 ff;\n";
+        compute<<"ff.x = (-dEdDihedral"<<cu.intToString(index)<<"*r)/"<<crossName1<<".w;\n";
+        compute<<"ff.y = (delta"<<deltaName1<<".x*delta"<<deltaName2<<".x + delta"<<deltaName1<<".y*delta"<<deltaName2<<".y + delta"<<deltaName1<<".z*delta"<<deltaName2<<".z)/delta"<<deltaName2<<".w;\n";
+        compute<<"ff.z = (delta"<<deltaName3<<".x*delta"<<deltaName2<<".x + delta"<<deltaName3<<".y*delta"<<deltaName2<<".y + delta"<<deltaName3<<".z*delta"<<deltaName2<<".z)/delta"<<deltaName2<<".w;\n";
+        compute<<"ff.w = (dEdDihedral"<<cu.intToString(index)<<"*r)/"<<crossName2<<".w;\n";
+        compute<<"real3 internalF0 = ff.x*trim("<<crossName1<<");\n";
+        compute<<"real3 internalF3 = ff.w*trim("<<crossName2<<");\n";
+        compute<<"real3 s = ff.y*internalF0 - ff.z*internalF3;\n";
+        compute<<forceNames[atoms[0]]<<" += internalF0;\n";
+        compute<<forceNames[atoms[1]]<<" += s-internalF0;\n";
+        compute<<forceNames[atoms[2]]<<" += -s-internalF3;\n";
+        compute<<forceNames[atoms[3]]<<" += internalF3;\n";
+        compute<<"}\n";
+    }
+    
+    // Store forces to global memory.
+    
+    for (int i = 0; i < particlesPerSet; i++)
+        compute<<"storeForce(atom"<<(i+1)<<", "<<forceNames[i]<<", forceBuffers);\n";
+    
+    // Create other replacements that depend on the number of particles per set.
+    
+    stringstream numCombinations, atomsForCombination, isValidCombination, permute, loadData, verifyCutoff, verifyExclusions;
+    if (hasTypeFilters) {
+        permute<<"int particleSet[] = {";
+        for (int i = 0; i < particlesPerSet; i++) {
+            permute<<"p"<<(i+1);
+            if (i < particlesPerSet-1)
+                permute<<", ";
+        }
+        permute<<"};\n";
+    }
+    for (int i = 0; i < particlesPerSet; i++) {
+        if (hasTypeFilters)
+            permute<<"int atom"<<(i+1)<<" = particleSet[particleOrder["<<particlesPerSet<<"*order+"<<i<<"]];\n";
+        else
+            permute<<"int atom"<<(i+1)<<" = p"<<(i+1)<<";\n";
+        loadData<<"real3 pos"<<(i+1)<<" = trim(posq[atom"<<(i+1)<<"]);\n";
+        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 (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 (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";
+        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++)
+            for (int j = i+1; j < particlesPerSet; j++)
+                verifyExclusions<<"includeInteraction &= !isInteractionExcluded(p"<<(i+1)<<", p"<<(j+1)<<", exclusions, exclusionStartIndex);\n";
+    }
+    string computeTypeIndex = "particleTypes[p"+cu.intToString(particlesPerSet)+"]";
+    for (int i = particlesPerSet-2; i >= 0; i--)
+        computeTypeIndex = "particleTypes[p"+cu.intToString(i+1)+"]+"+cu.intToString(numTypes)+"*("+computeTypeIndex+")";
+    
+    // Create replacements for extra arguments.
+    
+    stringstream extraArgs;
+    for (int i = 0; i < (int) params->getBuffers().size(); i++) {
+        CudaNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
+        extraArgs<<", const "<<buffer.getType()<<"* __restrict__ global_params"<<(i+1);
+    }
+
+    // Create the kernels.
+
+    map<string, string> replacements;
+    replacements["COMPUTE_INTERACTION"] = compute.str();
+    replacements["NUM_CANDIDATE_COMBINATIONS"] = numCombinations.str();
+    replacements["FIND_ATOMS_FOR_COMBINATION_INDEX"] = atomsForCombination.str();
+    replacements["IS_VALID_COMBINATION"] = isValidCombination.str();
+    replacements["VERIFY_CUTOFF"] = verifyCutoff.str();
+    replacements["VERIFY_EXCLUSIONS"] = verifyExclusions.str();
+    replacements["PERMUTE_ATOMS"] = permute.str();
+    replacements["LOAD_PARTICLE_DATA"] = loadData.str();
+    replacements["COMPUTE_TYPE_INDEX"] = computeTypeIndex;
+    replacements["PARAMETER_ARGUMENTS"] = extraArgs.str()+tableArgs.str();
+    map<string, string> defines;
+    if (nonbondedMethod != NoCutoff)
+        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)
+        defines["USE_EXCLUSIONS"] = "1";
+    defines["NUM_ATOMS"] = cu.intToString(cu.getNumAtoms());
+    defines["PADDED_NUM_ATOMS"] = cu.intToString(cu.getPaddedNumAtoms());
+    defines["M_PI"] = cu.doubleToString(M_PI);
+    defines["CUTOFF_SQUARED"] = cu.doubleToString(force.getCutoffDistance()*force.getCutoffDistance());
+    defines["TILE_SIZE"] = cu.intToString(CudaContext::TileSize);
+    defines["NUM_BLOCKS"] = cu.intToString(cu.getNumAtomBlocks());
+    defines["NUM_GLOBALS"] = cu.intToString(max(1, force.getNumGlobalParameters()));
+    defines["FIND_NEIGHBORS_WORKGROUP_SIZE"] = cu.intToString(findNeighborsWorkgroupSize);
+    CUmodule module = cu.createModule(cu.replaceStrings(CudaKernelSources::vectorOps+CudaKernelSources::customManyParticle, replacements), defines);
+    forceKernel = cu.getKernel(module, "computeInteraction");
+    blockBoundsKernel = cu.getKernel(module, "findBlockBounds");
+    neighborsKernel = cu.getKernel(module, "findNeighbors");
+    startIndicesKernel = cu.getKernel(module, "computeNeighborStartIndices");
+    copyPairsKernel = cu.getKernel(module, "copyPairsToNeighborList");
+    cuFuncSetCacheConfig(forceKernel, CU_FUNC_CACHE_PREFER_L1);
+    cuFuncSetCacheConfig(neighborsKernel, CU_FUNC_CACHE_PREFER_L1);
+    size_t bytes;
+    CHECK_RESULT(cuModuleGetGlobal(&globalsPtr, &bytes, module, "globals"), "Error getting address for constant memory")
+    cuMemcpyHtoD(globalsPtr, &globalParamValues[0], globalParamValues.size()*sizeof(float));
+}
+
+double CudaCalcCustomManyParticleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    if (!hasInitializedKernel) {
+        hasInitializedKernel = true;
+        
+        // Set arguments for the force kernel.
+        
+        forceArgs.push_back(&cu.getForce().getDevicePointer());
+        forceArgs.push_back(&cu.getEnergyBuffer().getDevicePointer());
+        forceArgs.push_back(&cu.getPosq().getDevicePointer());
+        forceArgs.push_back(cu.getPeriodicBoxSizePointer());
+        forceArgs.push_back(cu.getInvPeriodicBoxSizePointer());
+        if (nonbondedMethod != NoCutoff) {
+            forceArgs.push_back(&neighbors->getDevicePointer());
+            forceArgs.push_back(&neighborStartIndex->getDevicePointer());
+        }
+        if (particleTypes != NULL) {
+            forceArgs.push_back(&particleTypes->getDevicePointer());
+            forceArgs.push_back(&orderIndex->getDevicePointer());
+            forceArgs.push_back(&particleOrder->getDevicePointer());
+        }
+        if (exclusions != NULL) {
+            forceArgs.push_back(&exclusions->getDevicePointer());
+            forceArgs.push_back(&exclusionStartIndex->getDevicePointer());
+        }
+        for (int i = 0; i < (int) params->getBuffers().size(); i++) {
+            CudaNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
+            forceArgs.push_back(&buffer.getMemory());
+        }
+        for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
+            forceArgs.push_back(&tabulatedFunctions[i]->getDevicePointer());
+        
+        if (nonbondedMethod != NoCutoff) {
+            // Set arguments for the block bounds kernel.
+
+            blockBoundsArgs.push_back(cu.getPeriodicBoxSizePointer());
+            blockBoundsArgs.push_back(cu.getInvPeriodicBoxSizePointer());
+            blockBoundsArgs.push_back(&cu.getPosq().getDevicePointer());
+            blockBoundsArgs.push_back(&blockCenter->getDevicePointer());
+            blockBoundsArgs.push_back(&blockBoundingBox->getDevicePointer());
+            blockBoundsArgs.push_back(&numNeighborPairs->getDevicePointer());
+
+            // Set arguments for the neighbor list kernel.
+
+            neighborsArgs.push_back(cu.getPeriodicBoxSizePointer());
+            neighborsArgs.push_back(cu.getInvPeriodicBoxSizePointer());
+            neighborsArgs.push_back(&cu.getPosq().getDevicePointer());
+            neighborsArgs.push_back(&blockCenter->getDevicePointer());
+            neighborsArgs.push_back(&blockBoundingBox->getDevicePointer());
+            neighborsArgs.push_back(&neighborPairs->getDevicePointer());
+            neighborsArgs.push_back(&numNeighborPairs->getDevicePointer());
+            neighborsArgs.push_back(&numNeighborsForAtom->getDevicePointer());
+            neighborsArgs.push_back(&maxNeighborPairs);
+            if (exclusions != NULL) {
+                neighborsArgs.push_back(&exclusions->getDevicePointer());
+                neighborsArgs.push_back(&exclusionStartIndex->getDevicePointer());
+            }
+            
+            // Set arguments for the kernel to find neighbor list start indices.
+            
+            startIndicesArgs.push_back(&numNeighborsForAtom->getDevicePointer());
+            startIndicesArgs.push_back(&neighborStartIndex->getDevicePointer());
+            startIndicesArgs.push_back(&numNeighborPairs->getDevicePointer());
+            startIndicesArgs.push_back(&maxNeighborPairs);
+
+            // Set arguments for the kernel to assemble the final neighbor list.
+            
+            copyPairsArgs.push_back(&neighborPairs->getDevicePointer());
+            copyPairsArgs.push_back(&neighbors->getDevicePointer());
+            copyPairsArgs.push_back(&numNeighborPairs->getDevicePointer());
+            copyPairsArgs.push_back(&maxNeighborPairs);
+            copyPairsArgs.push_back(&numNeighborsForAtom->getDevicePointer());
+            copyPairsArgs.push_back(&neighborStartIndex->getDevicePointer());
+       }
+    }
+    if (globalParamValues.size() > 0) {
+        bool changed = false;
+        for (int i = 0; i < (int) globalParamNames.size(); i++) {
+            float value = (float) context.getParameter(globalParamNames[i]);
+            if (value != globalParamValues[i])
+                changed = true;
+            globalParamValues[i] = value;
+        }
+        if (changed)
+            cuMemcpyHtoD(globalsPtr, &globalParamValues[0], globalParamValues.size()*sizeof(float));
+    }
+    while (true) {
+        int* numPairs = (int*) cu.getPinnedBuffer();
+        if (nonbondedMethod != NoCutoff) {
+            cu.executeKernel(blockBoundsKernel, &blockBoundsArgs[0], cu.getNumAtomBlocks());
+            cu.executeKernel(neighborsKernel, &neighborsArgs[0], cu.getNumAtoms(), findNeighborsWorkgroupSize);
+
+            // We need to make sure there was enough memory for the neighbor list.  Download the
+            // information asynchronously so kernels can be running at the same time.
+
+            numNeighborPairs->download(numPairs, false);
+            CHECK_RESULT(cuEventRecord(event, 0), "Error recording event for CustomManyParticleForce");
+            cu.executeKernel(startIndicesKernel, &startIndicesArgs[0], 256, 256, 256*sizeof(int));
+            cu.executeKernel(copyPairsKernel, &copyPairsArgs[0], maxNeighborPairs);
+        }
+        int maxThreads = min(cu.getNumAtoms()*forceWorkgroupSize, cu.getEnergyBuffer().getSize());
+        cu.executeKernel(forceKernel, &forceArgs[0], maxThreads, forceWorkgroupSize);
+        if (nonbondedMethod != NoCutoff) {
+            // Make sure there was enough memory for the neighbor list.
+
+            CHECK_RESULT(cuEventSynchronize(event), "Error synchronizing on event for CustomManyParticleForce");
+            if (*numPairs > maxNeighborPairs) {
+                // Resize the arrays and run the calculation again.
+
+                delete neighborPairs;
+                neighborPairs = NULL;
+                delete neighbors;
+                neighbors = NULL;
+                maxNeighborPairs = (int) (1.1*(*numPairs));
+                neighborPairs = CudaArray::create<int2>(cu, maxNeighborPairs, "customManyParticleNeighborPairs");
+                neighbors = CudaArray::create<int>(cu, maxNeighborPairs, "customManyParticleNeighbors");
+                continue;
+            }
+        }
+        break;
+    }
+    return 0.0;
+}
+
+void CudaCalcCustomManyParticleForceKernel::copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force) {
+    cu.setAsCurrent();
+    int numParticles = force.getNumParticles();
+    if (numParticles != cu.getNumAtoms())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+    
+    // Record the per-particle parameters.
+    
+    vector<vector<float> > paramVector(numParticles);
+    vector<double> parameters;
+    int type;
+    for (int i = 0; i < numParticles; i++) {
+        force.getParticleParameters(i, parameters, type);
+        paramVector[i].resize(parameters.size());
+        for (int j = 0; j < (int) parameters.size(); j++)
+            paramVector[i][j] = (float) parameters[j];
+    }
+    params->setParameterValues(paramVector);
+    
+    // Mark that the current reordering may be invalid.
+    
+    cu.invalidateMolecules();
+}
+
 CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
 }
 

platforms/cuda/src/CudaPlatform.cpp

@@ -81,6 +81,7 @@ CudaPlatform::CudaPlatform() {
     registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);

platforms/cuda/src/kernels/customManyParticle.cu

+/**
+ * Record the force on an atom to global memory.
+ */
+inline __device__ void storeForce(int atom, real3 force, unsigned long long* __restrict__ forceBuffers) {
+    atomicAdd(&forceBuffers[atom], static_cast<unsigned long long>((long long) (force.x*0x100000000)));
+    atomicAdd(&forceBuffers[atom+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (force.y*0x100000000)));
+    atomicAdd(&forceBuffers[atom+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (force.z*0x100000000)));
+}
+
+/**
+ * Convert a real4 to a real3 by removing its last element.
+ */
+inline __device__ real3 trim(real4 v) {
+    return make_real3(v.x, v.y, v.z);
+}
+
+/**
+ * Compute the difference between two vectors, taking periodic boundary conditions into account
+ * and setting the fourth component to the squared magnitude.
+ */
+inline __device__ real4 delta(real3 vec1, real3 vec2, real4 periodicBoxSize, real4 invPeriodicBoxSize) {
+    real4 result = make_real4(vec1.x-vec2.x, vec1.y-vec2.y, vec1.z-vec2.z, 0.0f);
+#ifdef USE_PERIODIC
+    result.x -= floor(result.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+    result.y -= floor(result.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+    result.z -= floor(result.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
+#endif
+    result.w = result.x*result.x + result.y*result.y + result.z*result.z;
+    return result;
+}
+
+/**
+ * Compute the angle between two vectors.  The w component of each vector should contain the squared magnitude.
+ */
+__device__ real computeAngle(real4 vec1, real4 vec2) {
+    real dotProduct = vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
+    real cosine = dotProduct*RSQRT(vec1.w*vec2.w);
+    real angle;
+    if (cosine > 0.99f || cosine < -0.99f) {
+        // We're close to the singularity in acos(), so take the cross product and use asin() instead.
+
+        real3 crossProduct = cross(vec1, vec2);
+        real scale = vec1.w*vec2.w;
+        angle = ASIN(SQRT(dot(crossProduct, crossProduct)/scale));
+        if (cosine < 0.0f)
+            angle = M_PI-angle;
+    }
+    else
+       angle = ACOS(cosine);
+    return angle;
+}
+
+/**
+ * Compute the cross product of two vectors, setting the fourth component to the squared magnitude.
+ */
+inline __device__ real4 computeCross(real4 vec1, real4 vec2) {
+    real3 cp = cross(vec1, vec2);
+    return make_real4(cp.x, cp.y, cp.z, cp.x*cp.x+cp.y*cp.y+cp.z*cp.z);
+}
+
+/**
+ * Determine whether a particular interaction is in the list of exclusions.
+ */
+inline __device__ bool isInteractionExcluded(int atom1, int atom2, int* __restrict__ exclusions, int* __restrict__ exclusionStartIndex) {
+    int first = exclusionStartIndex[atom1];
+    int last = exclusionStartIndex[atom1+1];
+    for (int i = last-1; i >= first; i--) {
+        int excluded = exclusions[i];
+        if (excluded == atom2)
+            return true;
+        if (excluded <= atom1)
+            return false;
+    }
+    return false;
+}
+
+__constant__ float globals[NUM_GLOBALS];
+
+/**
+ * Compute the interaction.
+ */
+extern "C" __global__ void computeInteraction(
+        unsigned long long* __restrict__ forceBuffers, real* __restrict__ energyBuffer, const real4* __restrict__ posq,
+        real4 periodicBoxSize, real4 invPeriodicBoxSize
+#ifdef USE_CUTOFF
+        , const int* __restrict__ neighbors, const int* __restrict__ neighborStartIndex
+#endif
+#ifdef USE_FILTERS
+        , int* __restrict__ particleTypes, int* __restrict__ orderIndex, int* __restrict__ particleOrder
+#endif
+#ifdef USE_EXCLUSIONS
+        , int* __restrict__ exclusions, int* __restrict__ exclusionStartIndex
+#endif
+        PARAMETER_ARGUMENTS) {
+    real energy = 0.0f;
+    
+    // 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) {
+            FIND_ATOMS_FOR_COMBINATION_INDEX;
+            bool includeInteraction = IS_VALID_COMBINATION;
+#ifdef USE_CUTOFF
+            if (includeInteraction) {
+                VERIFY_CUTOFF;
+            }
+#endif
+#ifdef USE_FILTERS
+            int order = orderIndex[COMPUTE_TYPE_INDEX];
+            if (order == -1)
+                includeInteraction = false;
+#endif
+#ifdef USE_EXCLUSIONS
+            if (includeInteraction) {
+                VERIFY_EXCLUSIONS;
+            }
+#endif
+            if (includeInteraction) {
+                PERMUTE_ATOMS;
+                LOAD_PARTICLE_DATA;
+                COMPUTE_INTERACTION;
+            }
+        }
+    }
+    energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += energy;
+}
+
+/**
+ * Find a bounding box for the atoms in each block.
+ */
+extern "C" __global__ void findBlockBounds(real4 periodicBoxSize, real4 invPeriodicBoxSize, const real4* __restrict__ posq,
+        real4* __restrict__ blockCenter, real4* __restrict__ blockBoundingBox, int* __restrict__ numNeighborPairs) {
+    int index = blockIdx.x*blockDim.x+threadIdx.x;
+    int base = index*TILE_SIZE;
+    while (base < NUM_ATOMS) {
+        real4 pos = posq[base];
+#ifdef USE_PERIODIC
+        pos.x -= floor(pos.x*invPeriodicBoxSize.x)*periodicBoxSize.x;
+        pos.y -= floor(pos.y*invPeriodicBoxSize.y)*periodicBoxSize.y;
+        pos.z -= floor(pos.z*invPeriodicBoxSize.z)*periodicBoxSize.z;
+#endif
+        real4 minPos = pos;
+        real4 maxPos = pos;
+        int last = min(base+TILE_SIZE, NUM_ATOMS);
+        for (int i = base+1; i < last; i++) {
+            pos = posq[i];
+#ifdef USE_PERIODIC
+            real4 center = 0.5f*(maxPos+minPos);
+            pos.x -= floor((pos.x-center.x)*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+            pos.y -= floor((pos.y-center.y)*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+            pos.z -= floor((pos.z-center.z)*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
+#endif
+            minPos = make_real4(min(minPos.x,pos.x), min(minPos.y,pos.y), min(minPos.z,pos.z), 0);
+            maxPos = make_real4(max(maxPos.x,pos.x), max(maxPos.y,pos.y), max(maxPos.z,pos.z), 0);
+        }
+        real4 blockSize = 0.5f*(maxPos-minPos);
+        blockBoundingBox[index] = blockSize;
+        blockCenter[index] = 0.5f*(maxPos+minPos);
+        index += blockDim.x*gridDim.x;
+        base = index*TILE_SIZE;
+    }
+    if (blockIdx.x == 0 && threadIdx.x == 0)
+        *numNeighborPairs = 0;
+}
+
+/**
+ * Find a list of neighbors for each atom.
+ */
+extern "C" __global__ void findNeighbors(real4 periodicBoxSize, real4 invPeriodicBoxSize, const real4* __restrict__ posq,
+        const real4* __restrict__ blockCenter, const real4* __restrict__ blockBoundingBox, int2* __restrict__ neighborPairs,
+        int* __restrict__ numNeighborPairs, int* __restrict__ numNeighborsForAtom, int maxNeighborPairs
+#ifdef USE_EXCLUSIONS
+        , int* __restrict__ exclusions, int* __restrict__ exclusionStartIndex
+#endif
+        ) {
+    __shared__ real3 positionCache[FIND_NEIGHBORS_WORKGROUP_SIZE];
+    int indexInWarp = threadIdx.x%32;
+    for (int atom1 = blockIdx.x*blockDim.x+threadIdx.x; atom1 < PADDED_NUM_ATOMS; atom1 += blockDim.x*gridDim.x) {
+        // Load data for this atom.  Note that all threads in a warp are processing atoms from the same block.
+        
+        real3 pos1 = trim(posq[atom1]);
+        int block1 = atom1/TILE_SIZE;
+        real4 blockCenter1 = blockCenter[block1];
+        real4 blockSize1 = blockBoundingBox[block1];
+        int totalNeighborsForAtom1 = 0;
+        
+        // Loop over atom blocks to search for neighbors.  The threads in a warp compare block1 against 32
+        // other blocks in parallel.
+
+#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) {
+                real4 blockCenter2 = blockCenter[block2];
+                real4 blockSize2 = blockBoundingBox[block2];
+                real4 blockDelta = blockCenter1-blockCenter2;
+#ifdef USE_PERIODIC
+                blockDelta.x -= floor(blockDelta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                blockDelta.y -= floor(blockDelta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                blockDelta.z -= floor(blockDelta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
+#endif
+                blockDelta.x = max(0.0f, fabs(blockDelta.x)-blockSize1.x-blockSize2.x);
+                blockDelta.y = max(0.0f, fabs(blockDelta.y)-blockSize1.y-blockSize2.y);
+                blockDelta.z = max(0.0f, fabs(blockDelta.z)-blockSize1.z-blockSize2.z);
+                includeBlock2 &= (blockDelta.x*blockDelta.x+blockDelta.y*blockDelta.y+blockDelta.z*blockDelta.z < CUTOFF_SQUARED);
+            }
+            
+            // Loop over any blocks we identified as potentially containing neighbors.
+            
+            int includeBlockFlags = __ballot(includeBlock2);
+            while (includeBlockFlags != 0) {
+                int i = __ffs(includeBlockFlags)-1;
+                includeBlockFlags &= includeBlockFlags-1;
+                int block2 = block2Base+i;
+
+                // Loop over atoms in this block.
+
+                int start = block2*TILE_SIZE;
+                int included[TILE_SIZE];
+                int numIncluded = 0;
+                positionCache[threadIdx.x] = trim(posq[start+indexInWarp]);
+                if (atom1 < NUM_ATOMS) {
+                    for (int j = 0; j < 32; j++) {
+                        int atom2 = start+j;
+                        real3 pos2 = positionCache[threadIdx.x-indexInWarp+j];
+
+                        // 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);
+#endif
+                        if (includeAtom)
+                            included[numIncluded++] = atom2;
+                    }
+                }
+
+                // If we found any neighbors, store them to the neighbor list.
+
+                if (numIncluded > 0) {
+                    int baseIndex = atomicAdd(numNeighborPairs, numIncluded);
+                    if (baseIndex+numIncluded <= maxNeighborPairs)
+                        for (int j = 0; j < numIncluded; j++)
+                            neighborPairs[baseIndex+j] = make_int2(atom1, included[j]);
+                    totalNeighborsForAtom1 += numIncluded;
+                }
+            }
+        }
+        numNeighborsForAtom[atom1] = totalNeighborsForAtom1;
+    }
+}
+
+/**
+ * Sum the neighbor counts to compute the start position of each atom.  This kernel
+ * is executed as a single work group.
+ */
+extern "C" __global__ void computeNeighborStartIndices(int* __restrict__ numNeighborsForAtom, int* __restrict__ neighborStartIndex,
+            int* __restrict__ numNeighborPairs, int maxNeighborPairs) {
+    extern __shared__ unsigned int posBuffer[];
+    if (*numNeighborPairs > maxNeighborPairs) {
+        // There wasn't enough memory for the neighbor list, so we'll need to rebuild it.  Set the neighbor start
+        // indices to indicate no neighbors for any atom.
+        
+        for (int i = threadIdx.x; i <= NUM_ATOMS; i += blockDim.x)
+            neighborStartIndex[i] = 0;
+        return;
+    }
+    unsigned int globalOffset = 0;
+    for (unsigned int startAtom = 0; startAtom < NUM_ATOMS; startAtom += blockDim.x) {
+        // Load the neighbor counts into local memory.
+
+        unsigned int globalIndex = startAtom+threadIdx.x;
+        posBuffer[threadIdx.x] = (globalIndex < NUM_ATOMS ? numNeighborsForAtom[globalIndex] : 0);
+        __syncthreads();
+
+        // Perform a parallel prefix sum.
+
+        for (unsigned int step = 1; step < blockDim.x; step *= 2) {
+            unsigned int add = (threadIdx.x >= step ? posBuffer[threadIdx.x-step] : 0);
+            __syncthreads();
+            posBuffer[threadIdx.x] += add;
+            __syncthreads();
+        }
+
+        // Write the results back to global memory.
+
+        if (globalIndex < NUM_ATOMS) {
+            neighborStartIndex[globalIndex+1] = posBuffer[threadIdx.x]+globalOffset;
+            numNeighborsForAtom[globalIndex] = 0; // Clear this so the next kernel can use it as a counter
+        }
+        globalOffset += posBuffer[blockDim.x-1];
+    }
+    if (threadIdx.x == 0)
+        neighborStartIndex[0] = 0;
+}
+
+/**
+ * Assemble the final neighbor list.
+ */
+extern "C" __global__ void copyPairsToNeighborList(const int2* __restrict__ neighborPairs, int* __restrict__ neighbors, int* __restrict__ numNeighborPairs,
+            int maxNeighborPairs, int* __restrict__ numNeighborsForAtom, const int* __restrict__ neighborStartIndex) {
+    int actualPairs = *numNeighborPairs;
+    if (actualPairs > maxNeighborPairs)
+        return; // There wasn't enough memory for the neighbor list, so we'll need to rebuild it.
+    for (unsigned int index = blockDim.x*blockIdx.x+threadIdx.x; index < actualPairs; index += blockDim.x*gridDim.x) {
+        int2 pair = neighborPairs[index];
+        int startIndex = neighborStartIndex[pair.x];
+        int offset = atomicAdd(numNeighborsForAtom+pair.x, 1);
+        neighbors[startIndex+offset] = pair.y;
+    }
+}

platforms/cuda/tests/TestCudaCustomManyParticleForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2014 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the CUDA implementation of CustomManyParticleForce.
+ */
+
+#ifdef WIN32
+  #define _USE_MATH_DEFINES // Needed to get M_PI
+#endif
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "openmm/CustomCompoundBondForce.h"
+#include "openmm/CustomManyParticleForce.h"
+#include "openmm/System.h"
+#include "openmm/TabulatedFunction.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+CudaPlatform platform;
+
+void validateAxilrodTeller(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 c = context.getParameter("C");
+    
+    // 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);
+        double rprod = r12*r13*r23;
+        expectedEnergy += c*(1+3*ctheta1*ctheta2*ctheta3)/(rprod*rprod*rprod);
+    }
+    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 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;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    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[4][3] = {{0,1,2}, {1,2,3}, {2,3,0}, {3,0,1}};
+    vector<const int*> expectedSets(&sets[0], &sets[4]);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+}
+
+void testCutoff() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffNonPeriodic);
+    force->setCutoffDistance(1.55);
+    vector<double> params;
+    force->addParticle(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));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    int sets[7][3] = {{0,1,2}, {0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,2,4}, {1,3,4}};
+    vector<const int*> expectedSets(&sets[0], &sets[7]);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+}
+
+void testPeriodic() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
+    force->setCutoffDistance(1.05);
+    vector<double> params;
+    force->addParticle(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));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    double boxSize = 2.1;
+    int sets[5][3] = {{0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,3,4}};
+    vector<const int*> expectedSets(&sets[0], &sets[5]);
+    validateAxilrodTeller(force, positions, expectedSets, boxSize);
+}
+
+void testExclusions() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    vector<double> params;
+    force->addParticle(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));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    force->addExclusion(0, 2);
+    force->addExclusion(0, 3);
+    int sets[5][3] = {{0,1,4}, {1,2,3}, {1,2,4}, {1,3,4}, {2,3,4}};
+    vector<const int*> expectedSets(&sets[0], &sets[5]);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+}
+
+void testAllTerms() {
+    int numParticles = 4;
+    
+    // Create a system with a CustomManyParticleForce.
+    
+    System system1;
+    CustomManyParticleForce* force1 = new CustomManyParticleForce(4,
+        "distance(p1,p2)+angle(p1,p4,p3)+dihedral(p1,p3,p2,p4)+x1+y4+z3");
+    system1.addForce(force1);
+    vector<double> params;
+    for (int i = 0; i < numParticles; i++) {
+        system1.addParticle(1.0);
+        force1->addParticle(params, i);
+    }
+    set<int> filter;
+    filter.insert(0);
+    force1->setTypeFilter(0, filter);
+    filter.clear();
+    filter.insert(1);
+    force1->setTypeFilter(1, filter);
+    filter.clear();
+    filter.insert(3);
+    force1->setTypeFilter(2, filter);
+    filter.clear();
+    filter.insert(2);
+    force1->setTypeFilter(3, filter);
+    
+    // Create a system that use a CustomCompoundBondForce to compute exactly the same interactions.
+    
+    System system2;
+    CustomCompoundBondForce* force2 = new CustomCompoundBondForce(4,
+        "distance(p1,p2)+angle(p1,p3,p4)+dihedral(p1,p4,p2,p3)+x1+y3+z4");
+    system2.addForce(force2);
+    vector<int> particles;
+    particles.push_back(0);
+    particles.push_back(1);
+    particles.push_back(2);
+    particles.push_back(3);
+    force2->addBond(particles, params);
+    for (int i = 0; i < numParticles; i++)
+        system2.addParticle(1.0);
+    
+    // Create contexts for both of them.
+
+    vector<Vec3> positions;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; i++)
+        positions.push_back(Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
+    VerletIntegrator integrator1(0.001);
+    VerletIntegrator integrator2(0.001);
+    Context context1(system1, integrator1, platform);
+    Context context2(system2, integrator2, platform);
+    context1.setPositions(positions);
+    context2.setPositions(positions);
+    
+    // See if they produce identical forces and energies.
+    
+    State state1 = context1.getState(State::Forces | State::Energy);
+    State state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state2.getPotentialEnergy(), state1.getPotentialEnergy(), 1e-4);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state2.getForces()[i], state1.getForces()[i], 1e-4);
+}
+
+void testParameters() {
+    // Create a system.
+    
+    int numParticles = 5;
+    System system;
+    CustomManyParticleForce* force = new CustomManyParticleForce(3, "C*scale1*scale2*scale3*(distance(p1,p2)+distance(p2,p3)+distance(p1,p3))");
+    force->addGlobalParameter("C", 2.0);
+    force->addPerParticleParameter("scale");
+    vector<double> params(1);
+    vector<Vec3> positions;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; i++) {
+        params[0] = i+1;
+        force->addParticle(params);
+        positions.push_back(Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
+        system.addParticle(1.0);
+    }
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    
+    // See if the energy is correct.
+
+    State state = context.getState(State::Energy);
+    double expectedEnergy = 0;
+    for (int i = 0; i < numParticles; i++)
+        for (int j = i+1; j < numParticles; j++)
+            for (int k = j+1; k < numParticles; k++) {
+                Vec3 d12 = positions[j]-positions[i];
+                Vec3 d13 = positions[k]-positions[i];
+                Vec3 d23 = positions[k]-positions[j];
+                double r12 = sqrt(d12.dot(d12));
+                double r13 = sqrt(d13.dot(d13));
+                double r23 = sqrt(d23.dot(d23));
+                expectedEnergy += 2.0*(i+1)*(j+1)*(k+1)*(r12+r13+r23);
+            }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+    
+    // Modify the parameters.
+    
+    context.setParameter("C", 3.5);
+    for (int i = 0; i < numParticles; i++) {
+        params[0] = 0.5*i-0.1;
+        force->setParticleParameters(i, params, 0);
+    }
+    force->updateParametersInContext(context);
+    
+    // See if the energy is still correct.
+    
+    state = context.getState(State::Energy);
+    expectedEnergy = 0;
+    for (int i = 0; i < numParticles; i++)
+        for (int j = i+1; j < numParticles; j++)
+            for (int k = j+1; k < numParticles; k++) {
+                Vec3 d12 = positions[j]-positions[i];
+                Vec3 d13 = positions[k]-positions[i];
+                Vec3 d23 = positions[k]-positions[j];
+                double r12 = sqrt(d12.dot(d12));
+                double r13 = sqrt(d13.dot(d13));
+                double r23 = sqrt(d23.dot(d23));
+                expectedEnergy += 3.5*(0.5*i-0.1)*(0.5*j-0.1)*(0.5*k-0.1)*(r12+r13+r23);
+            }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+}
+
+void testTabulatedFunctions() {
+    int numParticles = 5;
+    
+    // Create two tabulated functions.
+    
+    vector<double> values;
+    values.push_back(0.0);
+    values.push_back(50.0);
+    Continuous1DFunction* f1 = new Continuous1DFunction(values, 0, 100);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<double> c(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        c[i] = genrand_real2(sfmt);
+    values.resize(numParticles*numParticles*numParticles);
+    for (int i = 0; i < numParticles; i++)
+        for (int j = 0; j < numParticles; j++)
+            for (int k = 0; k < numParticles; k++)
+                values[i+numParticles*j+numParticles*numParticles*k] = c[i]+c[j]+c[k];
+    Discrete3DFunction* f2 = new Discrete3DFunction(numParticles, numParticles, numParticles, values);
+    
+    // Create a system.
+    
+    System system;
+    CustomManyParticleForce* force = new CustomManyParticleForce(3, "f1(distance(p1,p2)+distance(p2,p3)+distance(p1,p3))*f2(atom1, atom2, atom3)");
+    force->addPerParticleParameter("atom");
+    force->addTabulatedFunction("f1", f1);
+    force->addTabulatedFunction("f2", f2);
+    vector<double> params(1);
+    vector<Vec3> positions;
+    for (int i = 0; i < numParticles; i++) {
+        params[0] = i;
+        force->addParticle(params);
+        positions.push_back(Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
+        system.addParticle(1.0);
+    }
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    
+    // See if the energy is correct.
+
+    State state = context.getState(State::Energy);
+    double expectedEnergy = 0;
+    for (int i = 0; i < numParticles; i++)
+        for (int j = i+1; j < numParticles; j++)
+            for (int k = j+1; k < numParticles; k++) {
+                Vec3 d12 = positions[j]-positions[i];
+                Vec3 d13 = positions[k]-positions[i];
+                Vec3 d23 = positions[k]-positions[j];
+                double r12 = sqrt(d12.dot(d12));
+                double r13 = sqrt(d13.dot(d13));
+                double r23 = sqrt(d23.dot(d23));
+                expectedEnergy += 0.5*(r12+r13+r23)*(c[i]+c[j]+c[k]);
+            }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+}
+
+void testTypeFilters() {
+    // Create a system.
+    
+    System system;
+    for (int i = 0; i < 5; i++)
+        system.addParticle(1.0);
+    CustomManyParticleForce* force = new CustomManyParticleForce(3, "c1*(distance(p1,p2)+distance(p1,p3))");
+    force->addPerParticleParameter("c");
+    double c[] = {1.0, 2.0, 1.3, 1.5, -2.1};
+    int type[] = {0, 1, 0, 1, 5};
+    vector<double> params(1);
+    for (int i = 0; i < 5; i++) {
+        params[0] = c[i];
+        force->addParticle(params, type[i]);
+    }
+    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));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    set<int> f1, f2;
+    f1.insert(0);
+    f2.insert(1);
+    f2.insert(5);
+    force->setTypeFilter(0, f1);
+    force->setTypeFilter(1, f2);
+    force->setTypeFilter(2, f2);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    
+    // See if the energy is correct.
+
+    State state = context.getState(State::Energy);
+    double expectedEnergy = 0;
+    int sets[6][3] = {{0,1,3}, {0,1,4}, {0,3,4}, {2,1,3}, {2,1,4}, {2,3,4}};
+    for (int i = 0; i < 6; i++) {
+        int p1 = sets[i][0];
+        int p2 = sets[i][1];
+        int p3 = sets[i][2];
+            Vec3 d12 = positions[p2]-positions[p1];
+            Vec3 d13 = positions[p3]-positions[p1];
+            double r12 = sqrt(d12.dot(d12));
+            double r13 = sqrt(d13.dot(d13));
+            expectedEnergy += c[p1]*(r12+r13);
+    }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+}
+
+void testLargeSystem() {
+    int gridSize = 8;
+    int numParticles = gridSize*gridSize*gridSize;
+    double boxSize = 3.0;
+    double spacing = boxSize/gridSize;
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
+    force->setCutoffDistance(0.6);
+    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);
+}
+
+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(0.1, 0, 0));
+    positions.push_back(Vec3(0, 0.11, 0.03));
+    positions.push_back(Vec3(0.04, 0, -0.08));
+    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(0.155);
+    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(0.1, 0, 0));
+    positions.push_back(Vec3(0, 0.11, 0.03));
+    positions.push_back(Vec3(0.04, 0, -0.08));
+    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)
+            platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
+        testNoCutoff();
+        testCutoff();
+        testPeriodic();
+        testExclusions();
+        testAllTerms();
+        testParameters();
+        testTabulatedFunctions();
+        testTypeFilters();
+        testLargeSystem();
+        testCentralParticleModeNoCutoff();
+        testCentralParticleModeCutoff();
+        testCentralParticleModeLargeSystem();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/opencl/include/OpenCLContext.h

@@ -32,6 +32,7 @@
 #include <string>
 #include <pthread.h>
 #define __CL_ENABLE_EXCEPTIONS
+#define CL_USE_DEPRECATED_OPENCL_1_1_APIS
 #ifdef _MSC_VER
     // Prevent Windows from defining macros that interfere with other code.
     #define NOMINMAX

platforms/opencl/include/OpenCLKernels.h

@@ -929,6 +929,67 @@ private:
     const System& system;
 };
 
+/**
+ * This kernel is invoked by CustomManyParticleForce to calculate the forces acting on the system.
+ */
+class OpenCLCalcCustomManyParticleForceKernel : public CalcCustomManyParticleForceKernel {
+public:
+    OpenCLCalcCustomManyParticleForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomManyParticleForceKernel(name, platform),
+            hasInitializedKernel(false), cl(cl), params(NULL), globals(NULL), particleTypes(NULL), orderIndex(NULL), particleOrder(NULL), exclusions(NULL),
+            exclusionStartIndex(NULL), blockCenter(NULL), blockBoundingBox(NULL), neighborPairs(NULL), numNeighborPairs(NULL), neighborStartIndex(NULL),
+            numNeighborsForAtom(NULL), neighbors(NULL), system(system) {
+    }
+    ~OpenCLCalcCustomManyParticleForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomManyParticleForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomManyParticleForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomManyParticleForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force);
+
+private:
+    OpenCLContext& cl;
+    bool hasInitializedKernel;
+    NonbondedMethod nonbondedMethod;
+    int maxNeighborPairs, forceWorkgroupSize, findNeighborsWorkgroupSize;
+    OpenCLParameterSet* params;
+    OpenCLArray* globals;
+    OpenCLArray* particleTypes;
+    OpenCLArray* orderIndex;
+    OpenCLArray* particleOrder;
+    OpenCLArray* exclusions;
+    OpenCLArray* exclusionStartIndex;
+    OpenCLArray* blockCenter;
+    OpenCLArray* blockBoundingBox;
+    OpenCLArray* neighborPairs;
+    OpenCLArray* numNeighborPairs;
+    OpenCLArray* neighborStartIndex;
+    OpenCLArray* numNeighborsForAtom;
+    OpenCLArray* neighbors;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+    std::vector<OpenCLArray*> tabulatedFunctions;
+    const System& system;
+    cl::Kernel forceKernel, blockBoundsKernel, neighborsKernel, startIndicesKernel, copyPairsKernel;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -104,6 +104,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
         return new OpenCLCalcCustomHbondForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcCustomCompoundBondForceKernel::Name())
         return new OpenCLCalcCustomCompoundBondForceKernel(name, platform, cl, context.getSystem());
+    if (name == CalcCustomManyParticleForceKernel::Name())
+        return new OpenCLCalcCustomManyParticleForceKernel(name, platform, cl, context.getSystem());
     if (name == IntegrateVerletStepKernel::Name())
         return new OpenCLIntegrateVerletStepKernel(name, platform, cl);
     if (name == IntegrateLangevinStepKernel::Name())

platforms/opencl/src/OpenCLKernels.cpp

@@ -33,6 +33,7 @@
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/CustomCompoundBondForceImpl.h"
 #include "openmm/internal/CustomHbondForceImpl.h"
+#include "openmm/internal/CustomManyParticleForceImpl.h"
 #include "openmm/internal/CustomNonbondedForceImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
 #include "OpenCLBondedUtilities.h"
@@ -4558,6 +4559,668 @@ void OpenCLCalcCustomCompoundBondForceKernel::copyParametersToContext(ContextImp
     cl.invalidateMolecules();
 }
 
+class OpenCLCustomManyParticleForceInfo : public OpenCLForceInfo {
+public:
+    OpenCLCustomManyParticleForceInfo(const CustomManyParticleForce& force) : OpenCLForceInfo(0), force(force) {
+    }
+    bool areParticlesIdentical(int particle1, int particle2) {
+        vector<double> params1, params2;
+        int type1, type2;
+        force.getParticleParameters(particle1, params1, type1);
+        force.getParticleParameters(particle2, params2, type2);
+        if (type1 != type2)
+            return false;
+        for (int i = 0; i < (int) params1.size(); i++)
+            if (params1[i] != params2[i])
+                return false;
+        return true;
+    }
+    int getNumParticleGroups() {
+        return force.getNumExclusions();
+    }
+    void getParticlesInGroup(int index, vector<int>& particles) {
+        int particle1, particle2;
+        force.getExclusionParticles(index, particle1, particle2);
+        particles.resize(2);
+        particles[0] = particle1;
+        particles[1] = particle2;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        return true;
+    }
+private:
+    const CustomManyParticleForce& force;
+};
+
+OpenCLCalcCustomManyParticleForceKernel::~OpenCLCalcCustomManyParticleForceKernel() {
+    if (params != NULL)
+        delete params;
+    if (globals != NULL)
+        delete globals;
+    if (orderIndex != NULL)
+        delete orderIndex;
+    if (particleOrder != NULL)
+        delete particleOrder;
+    if (particleTypes != NULL)
+        delete particleTypes;
+    if (exclusions != NULL)
+        delete exclusions;
+    if (exclusionStartIndex != NULL)
+        delete exclusionStartIndex;
+    if (blockCenter != NULL)
+        delete blockCenter;
+    if (blockBoundingBox != NULL)
+        delete blockBoundingBox;
+    if (neighborPairs != NULL)
+        delete neighborPairs;
+    if (numNeighborPairs != NULL)
+        delete numNeighborPairs;
+    if (neighborStartIndex != NULL)
+        delete neighborStartIndex;
+    if (neighbors != NULL)
+        delete neighbors;
+    if (numNeighborsForAtom != NULL)
+        delete numNeighborsForAtom;
+    for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
+        delete tabulatedFunctions[i];
+}
+
+void OpenCLCalcCustomManyParticleForceKernel::initialize(const System& system, const CustomManyParticleForce& force) {
+    if (!cl.getSupports64BitGlobalAtomics())
+        throw OpenMMException("CustomManyParticleForce requires a device that supports 64 bit atomic operations");
+    int numParticles = force.getNumParticles();
+    int particlesPerSet = force.getNumParticlesPerSet();
+    bool centralParticleMode = (force.getPermutationMode() == CustomManyParticleForce::UniqueCentralParticle);
+    nonbondedMethod = CalcCustomManyParticleForceKernel::NonbondedMethod(force.getNonbondedMethod());
+    forceWorkgroupSize = 128;
+    findNeighborsWorkgroupSize = (cl.getSIMDWidth() >= 32 ? 128 : 32);
+    
+    // Record parameter values.
+    
+    params = new OpenCLParameterSet(cl, force.getNumPerParticleParameters(), numParticles, "customManyParticleParameters");
+    vector<vector<float> > paramVector(numParticles);
+    for (int i = 0; i < numParticles; i++) {
+        vector<double> parameters;
+        int type;
+        force.getParticleParameters(i, parameters, type);
+        paramVector[i].resize(parameters.size());
+        for (int j = 0; j < (int) parameters.size(); j++)
+            paramVector[i][j] = (float) parameters[j];
+    }
+    params->setParameterValues(paramVector);
+    cl.addForce(new OpenCLCustomManyParticleForceInfo(force));
+
+    // Record the tabulated functions.
+
+    map<string, Lepton::CustomFunction*> functions;
+    vector<pair<string, string> > functionDefinitions;
+    vector<const TabulatedFunction*> functionList;
+    stringstream tableArgs;
+    for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
+        functionList.push_back(&force.getTabulatedFunction(i));
+        string name = force.getTabulatedFunctionName(i);
+        string arrayName = "table"+cl.intToString(i);
+        functionDefinitions.push_back(make_pair(name, arrayName));
+        functions[name] = cl.getExpressionUtilities().getFunctionPlaceholder(force.getTabulatedFunction(i));
+        int width;
+        vector<float> f = cl.getExpressionUtilities().computeFunctionCoefficients(force.getTabulatedFunction(i), width);
+        tabulatedFunctions.push_back(OpenCLArray::create<float>(cl, f.size(), "TabulatedFunction"));
+        tabulatedFunctions[tabulatedFunctions.size()-1]->upload(f);
+        tableArgs << ", __global const float";
+        if (width > 1)
+            tableArgs << width;
+        tableArgs << "* restrict " << arrayName;
+    }
+    
+    // Record information about parameters.
+
+    globalParamNames.resize(force.getNumGlobalParameters());
+    globalParamValues.resize(force.getNumGlobalParameters());
+    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+        globalParamNames[i] = force.getGlobalParameterName(i);
+        globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
+    }
+    vector<pair<ExpressionTreeNode, string> > variables;
+    for (int i = 0; i < particlesPerSet; i++) {
+        string index = cl.intToString(i+1);
+        variables.push_back(makeVariable("x"+index, "pos"+index+".x"));
+        variables.push_back(makeVariable("y"+index, "pos"+index+".y"));
+        variables.push_back(makeVariable("z"+index, "pos"+index+".z"));
+    }
+    for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
+        const string& name = force.getPerParticleParameterName(i);
+        for (int j = 0; j < particlesPerSet; j++) {
+            string index = cl.intToString(j+1);
+            variables.push_back(makeVariable(name+index, "params"+params->getParameterSuffix(i, index)));
+        }
+    }
+    if (force.getNumGlobalParameters() > 0) {
+        globals = OpenCLArray::create<cl_float>(cl, force.getNumGlobalParameters(), "customManyParticleGlobals", CL_MEM_READ_ONLY);
+        globals->upload(globalParamValues);
+        for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+            const string& name = force.getGlobalParameterName(i);
+            string value = "globals["+cl.intToString(i)+"]";
+            variables.push_back(makeVariable(name, value));
+        }
+    }
+    
+    // Build data structures for type filters.
+    
+    vector<int> particleTypesVec;
+    vector<int> orderIndexVec;
+    vector<std::vector<int> > particleOrderVec;
+    int numTypes;
+    CustomManyParticleForceImpl::buildFilterArrays(force, numTypes, particleTypesVec, orderIndexVec, particleOrderVec);
+    bool hasTypeFilters = (particleOrderVec.size() > 1);
+    if (hasTypeFilters) {
+        particleTypes = OpenCLArray::create<int>(cl, particleTypesVec.size(), "customManyParticleTypes");
+        orderIndex = OpenCLArray::create<int>(cl, orderIndexVec.size(), "customManyParticleOrderIndex");
+        particleOrder = OpenCLArray::create<int>(cl, particleOrderVec.size()*particlesPerSet, "customManyParticleOrder");
+        particleTypes->upload(particleTypesVec);
+        orderIndex->upload(orderIndexVec);
+        vector<int> flattenedOrder(particleOrder->getSize());
+        for (int i = 0; i < (int) particleOrderVec.size(); i++)
+            for (int j = 0; j < particlesPerSet; j++)
+                flattenedOrder[i*particlesPerSet+j] = particleOrderVec[i][j];
+        particleOrder->upload(flattenedOrder);
+    }
+    
+    // Build data structures for exclusions.
+    
+    if (force.getNumExclusions() > 0) {
+        vector<vector<int> > particleExclusions(numParticles);
+        for (int i = 0; i < force.getNumExclusions(); i++) {
+            int p1, p2;
+            force.getExclusionParticles(i, p1, p2);
+            particleExclusions[p1].push_back(p2);
+            particleExclusions[p2].push_back(p1);
+        }
+        vector<int> exclusionsVec;
+        vector<int> exclusionStartIndexVec(numParticles+1);
+        exclusionStartIndexVec[0] = 0;
+        for (int i = 0; i < numParticles; i++) {
+            sort(particleExclusions[i].begin(), particleExclusions[i].end());
+            exclusionsVec.insert(exclusionsVec.end(), particleExclusions[i].begin(), particleExclusions[i].end());
+            exclusionStartIndexVec[i+1] = exclusionsVec.size();
+        }
+        exclusions = OpenCLArray::create<int>(cl, exclusionsVec.size(), "customManyParticleExclusions");
+        exclusionStartIndex = OpenCLArray::create<int>(cl, exclusionStartIndexVec.size(), "customManyParticleExclusionStart");
+        exclusions->upload(exclusionsVec);
+        exclusionStartIndex->upload(exclusionStartIndexVec);
+    }
+    
+    // Build data structures for the neighbor list.
+    
+    if (nonbondedMethod != NoCutoff) {
+        int numAtomBlocks = cl.getNumAtomBlocks();
+        int elementSize = (cl.getUseDoublePrecision() ? sizeof(double) : sizeof(float));
+        blockCenter = new OpenCLArray(cl, numAtomBlocks, 4*elementSize, "blockCenter");
+        blockBoundingBox = new OpenCLArray(cl, numAtomBlocks, 4*elementSize, "blockBoundingBox");
+        numNeighborPairs = OpenCLArray::create<int>(cl, 1, "customManyParticleNumNeighborPairs");
+        neighborStartIndex = OpenCLArray::create<int>(cl, numParticles+1, "customManyParticleNeighborStartIndex");
+        numNeighborsForAtom = OpenCLArray::create<int>(cl, numParticles, "customManyParticleNumNeighborsForAtom");
+
+        // Select a size for the array that holds the neighbor list.  We have to make a fairly
+        // arbitrary guess, but if this turns out to be too small we'll increase it later.
+
+        maxNeighborPairs = 150*numParticles;
+        neighborPairs = OpenCLArray::create<mm_int2>(cl, maxNeighborPairs, "customManyParticleNeighborPairs");
+        neighbors = OpenCLArray::create<int>(cl, maxNeighborPairs, "customManyParticleNeighbors");
+    }
+
+    // Now to generate the kernel.  First, it needs to calculate all distances, angles,
+    // and dihedrals the expression depends on.
+
+    map<string, vector<int> > distances;
+    map<string, vector<int> > angles;
+    map<string, vector<int> > dihedrals;
+    Lepton::ParsedExpression energyExpression = CustomManyParticleForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
+    map<string, Lepton::ParsedExpression> forceExpressions;
+    set<string> computedDeltas;
+    vector<string> atomNames, posNames;
+    for (int i = 0; i < particlesPerSet; i++) {
+        string index = cl.intToString(i+1);
+        atomNames.push_back("P"+index);
+        posNames.push_back("pos"+index);
+    }
+    stringstream compute;
+    int index = 0;
+    for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName = atomNames[atoms[0]]+atomNames[atoms[1]];
+        if (computedDeltas.count(deltaName) == 0) {
+            compute<<"real4 delta"<<deltaName<<" = delta("<<posNames[atoms[0]]<<", "<<posNames[atoms[1]]<<", periodicBoxSize, invPeriodicBoxSize);\n";
+            computedDeltas.insert(deltaName);
+        }
+        compute<<"real r_"<<deltaName<<" = sqrt(delta"<<deltaName<<".w);\n";
+        variables.push_back(makeVariable(iter->first, "r_"+deltaName));
+        forceExpressions["real dEdDistance"+cl.intToString(index)+" = "] = energyExpression.differentiate(iter->first).optimize();
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName1 = atomNames[atoms[1]]+atomNames[atoms[0]];
+        string deltaName2 = atomNames[atoms[1]]+atomNames[atoms[2]];
+        string angleName = "angle_"+atomNames[atoms[0]]+atomNames[atoms[1]]+atomNames[atoms[2]];
+        if (computedDeltas.count(deltaName1) == 0) {
+            compute<<"real4 delta"<<deltaName1<<" = delta("<<posNames[atoms[1]]<<", "<<posNames[atoms[0]]<<", periodicBoxSize, invPeriodicBoxSize);\n";
+            computedDeltas.insert(deltaName1);
+        }
+        if (computedDeltas.count(deltaName2) == 0) {
+            compute<<"real4 delta"<<deltaName2<<" = delta("<<posNames[atoms[1]]<<", "<<posNames[atoms[2]]<<", periodicBoxSize, invPeriodicBoxSize);\n";
+            computedDeltas.insert(deltaName2);
+        }
+        compute<<"real "<<angleName<<" = computeAngle(delta"<<deltaName1<<", delta"<<deltaName2<<");\n";
+        variables.push_back(makeVariable(iter->first, angleName));
+        forceExpressions["real dEdAngle"+cl.intToString(index)+" = "] = energyExpression.differentiate(iter->first).optimize();
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName1 = atomNames[atoms[0]]+atomNames[atoms[1]];
+        string deltaName2 = atomNames[atoms[2]]+atomNames[atoms[1]];
+        string deltaName3 = atomNames[atoms[2]]+atomNames[atoms[3]];
+        string crossName1 = "cross_"+deltaName1+"_"+deltaName2;
+        string crossName2 = "cross_"+deltaName2+"_"+deltaName3;
+        string dihedralName = "dihedral_"+atomNames[atoms[0]]+atomNames[atoms[1]]+atomNames[atoms[2]]+atomNames[atoms[3]];
+        if (computedDeltas.count(deltaName1) == 0) {
+            compute<<"real4 delta"<<deltaName1<<" = delta("<<posNames[atoms[0]]<<", "<<posNames[atoms[1]]<<", periodicBoxSize, invPeriodicBoxSize);\n";
+            computedDeltas.insert(deltaName1);
+        }
+        if (computedDeltas.count(deltaName2) == 0) {
+            compute<<"real4 delta"<<deltaName2<<" = delta("<<posNames[atoms[2]]<<", "<<posNames[atoms[1]]<<", periodicBoxSize, invPeriodicBoxSize);\n";
+            computedDeltas.insert(deltaName2);
+        }
+        if (computedDeltas.count(deltaName3) == 0) {
+            compute<<"real4 delta"<<deltaName3<<" = delta("<<posNames[atoms[2]]<<", "<<posNames[atoms[3]]<<", periodicBoxSize, invPeriodicBoxSize);\n";
+            computedDeltas.insert(deltaName3);
+        }
+        compute<<"real4 "<<crossName1<<" = computeCross(delta"<<deltaName1<<", delta"<<deltaName2<<");\n";
+        compute<<"real4 "<<crossName2<<" = computeCross(delta"<<deltaName2<<", delta"<<deltaName3<<");\n";
+        compute<<"real "<<dihedralName<<" = computeAngle("<<crossName1<<", "<<crossName2<<");\n";
+        compute<<dihedralName<<" *= (delta"<<deltaName1<<".x*"<<crossName2<<".x + delta"<<deltaName1<<".y*"<<crossName2<<".y + delta"<<deltaName1<<".z*"<<crossName2<<".z < 0 ? -1 : 1);\n";
+        variables.push_back(makeVariable(iter->first, dihedralName));
+        forceExpressions["real dEdDihedral"+cl.intToString(index)+" = "] = energyExpression.differentiate(iter->first).optimize();
+    }
+
+    // Now evaluate the expressions.
+
+    for (int i = 0; i < (int) params->getBuffers().size(); i++) {
+        OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
+        compute<<buffer.getType()<<" params"<<(i+1)<<" = global_params"<<(i+1)<<"[index];\n";
+    }
+    forceExpressions["energy += "] = energyExpression;
+    compute << cl.getExpressionUtilities().createExpressions(forceExpressions, variables, functionList, functionDefinitions, "temp");
+
+    // Apply forces to atoms.
+
+    vector<string> forceNames;
+    for (int i = 0; i < particlesPerSet; i++) {
+        string istr = cl.intToString(i+1);
+        string forceName = "force"+istr;
+        forceNames.push_back(forceName);
+        compute<<"real4 "<<forceName<<" = (real4) 0;\n";
+        compute<<"{\n";
+        Lepton::ParsedExpression forceExpressionX = energyExpression.differentiate("x"+istr).optimize();
+        Lepton::ParsedExpression forceExpressionY = energyExpression.differentiate("y"+istr).optimize();
+        Lepton::ParsedExpression forceExpressionZ = energyExpression.differentiate("z"+istr).optimize();
+        map<string, Lepton::ParsedExpression> expressions;
+        if (!isZeroExpression(forceExpressionX))
+            expressions[forceName+".x -= "] = forceExpressionX;
+        if (!isZeroExpression(forceExpressionY))
+            expressions[forceName+".y -= "] = forceExpressionY;
+        if (!isZeroExpression(forceExpressionZ))
+            expressions[forceName+".z -= "] = forceExpressionZ;
+        if (expressions.size() > 0)
+            compute<<cl.getExpressionUtilities().createExpressions(expressions, variables, functionList, functionDefinitions, "coordtemp");
+        compute<<"}\n";
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName = atomNames[atoms[0]]+atomNames[atoms[1]];
+        string value = "(dEdDistance"+cl.intToString(index)+"/r_"+deltaName+")*delta"+deltaName+".xyz";
+        compute<<forceNames[atoms[0]]<<".xyz += "<<"-"<<value<<";\n";
+        compute<<forceNames[atoms[1]]<<".xyz += "<<value<<";\n";
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName1 = atomNames[atoms[1]]+atomNames[atoms[0]];
+        string deltaName2 = atomNames[atoms[1]]+atomNames[atoms[2]];
+        compute<<"{\n";
+        compute<<"real4 crossProd = cross(delta"<<deltaName2<<", delta"<<deltaName1<<");\n";
+        compute<<"real lengthCross = max(SQRT(dot(crossProd, crossProd)), (real) 1e-6f);\n";
+        compute<<"real4 deltaCross0 = -cross(delta"<<deltaName1<<", crossProd)*dEdAngle"<<cl.intToString(index)<<"/(delta"<<deltaName1<<".w*lengthCross);\n";
+        compute<<"real4 deltaCross2 = cross(delta"<<deltaName2<<", crossProd)*dEdAngle"<<cl.intToString(index)<<"/(delta"<<deltaName2<<".w*lengthCross);\n";
+        compute<<"real4 deltaCross1 = -(deltaCross0+deltaCross2);\n";
+        compute<<forceNames[atoms[0]]<<".xyz += deltaCross0.xyz;\n";
+        compute<<forceNames[atoms[1]]<<".xyz += deltaCross1.xyz;\n";
+        compute<<forceNames[atoms[2]]<<".xyz += deltaCross2.xyz;\n";
+        compute<<"}\n";
+    }
+    index = 0;
+    for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter, ++index) {
+        const vector<int>& atoms = iter->second;
+        string deltaName1 = atomNames[atoms[0]]+atomNames[atoms[1]];
+        string deltaName2 = atomNames[atoms[2]]+atomNames[atoms[1]];
+        string deltaName3 = atomNames[atoms[2]]+atomNames[atoms[3]];
+        string crossName1 = "cross_"+deltaName1+"_"+deltaName2;
+        string crossName2 = "cross_"+deltaName2+"_"+deltaName3;
+        compute<<"{\n";
+        compute<<"real r = sqrt(delta"<<deltaName2<<".w);\n";
+        compute<<"real4 ff;\n";
+        compute<<"ff.x = (-dEdDihedral"<<cl.intToString(index)<<"*r)/"<<crossName1<<".w;\n";
+        compute<<"ff.y = (delta"<<deltaName1<<".x*delta"<<deltaName2<<".x + delta"<<deltaName1<<".y*delta"<<deltaName2<<".y + delta"<<deltaName1<<".z*delta"<<deltaName2<<".z)/delta"<<deltaName2<<".w;\n";
+        compute<<"ff.z = (delta"<<deltaName3<<".x*delta"<<deltaName2<<".x + delta"<<deltaName3<<".y*delta"<<deltaName2<<".y + delta"<<deltaName3<<".z*delta"<<deltaName2<<".z)/delta"<<deltaName2<<".w;\n";
+        compute<<"ff.w = (dEdDihedral"<<cl.intToString(index)<<"*r)/"<<crossName2<<".w;\n";
+        compute<<"real4 internalF0 = ff.x*"<<crossName1<<";\n";
+        compute<<"real4 internalF3 = ff.w*"<<crossName2<<";\n";
+        compute<<"real4 s = ff.y*internalF0 - ff.z*internalF3;\n";
+        compute<<forceNames[atoms[0]]<<".xyz += internalF0.xyz;\n";
+        compute<<forceNames[atoms[1]]<<".xyz += s.xyz-internalF0.xyz;\n";
+        compute<<forceNames[atoms[2]]<<".xyz += -s.xyz-internalF3.xyz;\n";
+        compute<<forceNames[atoms[3]]<<".xyz += internalF3.xyz;\n";
+        compute<<"}\n";
+    }
+    
+    // Store forces to global memory.
+    
+    for (int i = 0; i < particlesPerSet; i++)
+        compute<<"storeForce(atom"<<(i+1)<<", "<<forceNames[i]<<", forceBuffers);\n";
+    
+    // Create other replacements that depend on the number of particles per set.
+    
+    stringstream numCombinations, atomsForCombination, isValidCombination, permute, loadData, verifyCutoff, verifyExclusions;
+    if (hasTypeFilters) {
+        permute<<"int particleSet[] = {";
+        for (int i = 0; i < particlesPerSet; i++) {
+            permute<<"p"<<(i+1);
+            if (i < particlesPerSet-1)
+                permute<<", ";
+        }
+        permute<<"};\n";
+    }
+    for (int i = 0; i < particlesPerSet; i++) {
+        if (hasTypeFilters)
+            permute<<"int atom"<<(i+1)<<" = particleSet[particleOrder["<<particlesPerSet<<"*order+"<<i<<"]];\n";
+        else
+            permute<<"int atom"<<(i+1)<<" = p"<<(i+1)<<";\n";
+        loadData<<"real4 pos"<<(i+1)<<" = posq[atom"<<(i+1)<<"];\n";
+        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 (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 (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<<"real4 pos"<<(i+1)<<" = posq[p"<<(i+1)<<"];\n";
+        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++)
+            for (int j = i+1; j < particlesPerSet; j++)
+                verifyExclusions<<"includeInteraction &= !isInteractionExcluded(p"<<(i+1)<<", p"<<(j+1)<<", exclusions, exclusionStartIndex);\n";
+    }
+    string computeTypeIndex = "particleTypes[p"+cl.intToString(particlesPerSet)+"]";
+    for (int i = particlesPerSet-2; i >= 0; i--)
+        computeTypeIndex = "particleTypes[p"+cl.intToString(i+1)+"]+"+cl.intToString(numTypes)+"*("+computeTypeIndex+")";
+    
+    // Create replacements for extra arguments.
+    
+    stringstream extraArgs;
+    if (force.getNumGlobalParameters() > 0)
+        extraArgs << ", __global const float* globals";
+    for (int i = 0; i < (int) params->getBuffers().size(); i++) {
+        OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
+        extraArgs<<", __global const "<<buffer.getType()<<"* restrict global_params"<<(i+1);
+    }
+
+    // Create the kernels.
+
+    map<string, string> replacements;
+    replacements["COMPUTE_INTERACTION"] = compute.str();
+    replacements["NUM_CANDIDATE_COMBINATIONS"] = numCombinations.str();
+    replacements["FIND_ATOMS_FOR_COMBINATION_INDEX"] = atomsForCombination.str();
+    replacements["IS_VALID_COMBINATION"] = isValidCombination.str();
+    replacements["VERIFY_CUTOFF"] = verifyCutoff.str();
+    replacements["VERIFY_EXCLUSIONS"] = verifyExclusions.str();
+    replacements["PERMUTE_ATOMS"] = permute.str();
+    replacements["LOAD_PARTICLE_DATA"] = loadData.str();
+    replacements["COMPUTE_TYPE_INDEX"] = computeTypeIndex;
+    replacements["PARAMETER_ARGUMENTS"] = extraArgs.str()+tableArgs.str();
+    map<string, string> defines;
+    if (nonbondedMethod != NoCutoff)
+        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)
+        defines["USE_EXCLUSIONS"] = "1";
+    defines["NUM_ATOMS"] = cl.intToString(cl.getNumAtoms());
+    defines["PADDED_NUM_ATOMS"] = cl.intToString(cl.getPaddedNumAtoms());
+    defines["M_PI"] = cl.doubleToString(M_PI);
+    defines["CUTOFF_SQUARED"] = cl.doubleToString(force.getCutoffDistance()*force.getCutoffDistance());
+    defines["TILE_SIZE"] = cl.intToString(OpenCLContext::TileSize);
+    defines["NUM_BLOCKS"] = cl.intToString(cl.getNumAtomBlocks());
+    defines["FIND_NEIGHBORS_WORKGROUP_SIZE"] = cl.intToString(findNeighborsWorkgroupSize);
+    cl::Program program = cl.createProgram(cl.replaceStrings(OpenCLKernelSources::customManyParticle, replacements), defines);
+    forceKernel = cl::Kernel(program, "computeInteraction");
+    blockBoundsKernel = cl::Kernel(program, "findBlockBounds");
+    neighborsKernel = cl::Kernel(program, "findNeighbors");
+    startIndicesKernel = cl::Kernel(program, "computeNeighborStartIndices");
+    copyPairsKernel = cl::Kernel(program, "copyPairsToNeighborList");
+}
+
+double OpenCLCalcCustomManyParticleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    if (!hasInitializedKernel) {
+        hasInitializedKernel = true;
+        
+        // Set arguments for the force kernel.
+        
+        int index = 0;
+        forceKernel.setArg<cl::Buffer>(index++, cl.getLongForceBuffer().getDeviceBuffer());
+        forceKernel.setArg<cl::Buffer>(index++, cl.getEnergyBuffer().getDeviceBuffer());
+        forceKernel.setArg<cl::Buffer>(index++, cl.getPosq().getDeviceBuffer());
+        setPeriodicBoxSizeArg(cl, forceKernel, index++);
+        setInvPeriodicBoxSizeArg(cl, forceKernel, index++);
+        if (nonbondedMethod != NoCutoff) {
+            forceKernel.setArg<cl::Buffer>(index++, neighbors->getDeviceBuffer());
+            forceKernel.setArg<cl::Buffer>(index++, neighborStartIndex->getDeviceBuffer());
+        }
+        if (particleTypes != NULL) {
+            forceKernel.setArg<cl::Buffer>(index++, particleTypes->getDeviceBuffer());
+            forceKernel.setArg<cl::Buffer>(index++, orderIndex->getDeviceBuffer());
+            forceKernel.setArg<cl::Buffer>(index++, particleOrder->getDeviceBuffer());
+        }
+        if (exclusions != NULL) {
+            forceKernel.setArg<cl::Buffer>(index++, exclusions->getDeviceBuffer());
+            forceKernel.setArg<cl::Buffer>(index++, exclusionStartIndex->getDeviceBuffer());
+        }
+        if (globals != NULL)
+            forceKernel.setArg<cl::Buffer>(index++, globals->getDeviceBuffer());
+        for (int i = 0; i < (int) params->getBuffers().size(); i++) {
+            OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
+            forceKernel.setArg<cl::Memory>(index++, buffer.getMemory());
+        }
+        for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
+            forceKernel.setArg<cl::Buffer>(index++, tabulatedFunctions[i]->getDeviceBuffer());
+        
+        if (nonbondedMethod != NoCutoff) {
+            // Set arguments for the block bounds kernel.
+
+            index = 0;
+            setPeriodicBoxSizeArg(cl, blockBoundsKernel, index++);
+            setInvPeriodicBoxSizeArg(cl, blockBoundsKernel, index++);
+            blockBoundsKernel.setArg<cl::Buffer>(index++, cl.getPosq().getDeviceBuffer());
+            blockBoundsKernel.setArg<cl::Buffer>(index++, blockCenter->getDeviceBuffer());
+            blockBoundsKernel.setArg<cl::Buffer>(index++, blockBoundingBox->getDeviceBuffer());
+            blockBoundsKernel.setArg<cl::Buffer>(index++, numNeighborPairs->getDeviceBuffer());
+
+            // Set arguments for the neighbor list kernel.
+
+            index = 0;
+            setPeriodicBoxSizeArg(cl, neighborsKernel, index++);
+            setInvPeriodicBoxSizeArg(cl, neighborsKernel, index++);
+            neighborsKernel.setArg<cl::Buffer>(index++, cl.getPosq().getDeviceBuffer());
+            neighborsKernel.setArg<cl::Buffer>(index++, blockCenter->getDeviceBuffer());
+            neighborsKernel.setArg<cl::Buffer>(index++, blockBoundingBox->getDeviceBuffer());
+            neighborsKernel.setArg<cl::Buffer>(index++, neighborPairs->getDeviceBuffer());
+            neighborsKernel.setArg<cl::Buffer>(index++, numNeighborPairs->getDeviceBuffer());
+            neighborsKernel.setArg<cl::Buffer>(index++, numNeighborsForAtom->getDeviceBuffer());
+            index++;
+            if (exclusions != NULL) {
+                neighborsKernel.setArg<cl::Buffer>(index++, exclusions->getDeviceBuffer());
+                neighborsKernel.setArg<cl::Buffer>(index++, exclusionStartIndex->getDeviceBuffer());
+            }
+            
+            // Set arguments for the kernel to find neighbor list start indices.
+            
+            index = 0;
+            startIndicesKernel.setArg<cl::Buffer>(index++, numNeighborsForAtom->getDeviceBuffer());
+            startIndicesKernel.setArg<cl::Buffer>(index++, neighborStartIndex->getDeviceBuffer());
+            startIndicesKernel.setArg<cl::Buffer>(index++, numNeighborPairs->getDeviceBuffer());
+
+            // Set arguments for the kernel to assemble the final neighbor list.
+            
+            index = 0;
+            copyPairsKernel.setArg<cl::Buffer>(index++, neighborPairs->getDeviceBuffer());
+            copyPairsKernel.setArg<cl::Buffer>(index++, neighbors->getDeviceBuffer());
+            copyPairsKernel.setArg<cl::Buffer>(index++, numNeighborPairs->getDeviceBuffer());
+            index++;
+            copyPairsKernel.setArg<cl::Buffer>(index++, numNeighborsForAtom->getDeviceBuffer());
+            copyPairsKernel.setArg<cl::Buffer>(index++, neighborStartIndex->getDeviceBuffer());
+       }
+    }
+    if (globals != NULL) {
+        bool changed = false;
+        for (int i = 0; i < (int) globalParamNames.size(); i++) {
+            cl_float value = (cl_float) context.getParameter(globalParamNames[i]);
+            if (value != globalParamValues[i])
+                changed = true;
+            globalParamValues[i] = value;
+        }
+        if (changed)
+            globals->upload(globalParamValues);
+    }
+    while (true) {
+        int* numPairs = (int*) cl.getPinnedBuffer();
+        cl::Event event;
+        if (nonbondedMethod != NoCutoff) {
+            neighborsKernel.setArg<int>(8, maxNeighborPairs);
+            startIndicesKernel.setArg<int>(3, maxNeighborPairs);
+            copyPairsKernel.setArg<int>(3, maxNeighborPairs);
+            cl.executeKernel(blockBoundsKernel, cl.getNumAtomBlocks());
+            cl.executeKernel(neighborsKernel, cl.getNumAtoms(), findNeighborsWorkgroupSize);
+
+            // We need to make sure there was enough memory for the neighbor list.  Download the
+            // information asynchronously so kernels can be running at the same time.
+
+            numNeighborPairs->download(numPairs, false);
+            cl.getQueue().enqueueMarker(&event);
+            cl.executeKernel(startIndicesKernel, 256, 256);
+            cl.executeKernel(copyPairsKernel, maxNeighborPairs);
+        }
+        int maxThreads = min(cl.getNumAtoms()*forceWorkgroupSize, cl.getEnergyBuffer().getSize());
+        cl.executeKernel(forceKernel, maxThreads, forceWorkgroupSize);
+        if (nonbondedMethod != NoCutoff) {
+            // Make sure there was enough memory for the neighbor list.
+
+            event.wait();
+            if (*numPairs > maxNeighborPairs) {
+                // Resize the arrays and run the calculation again.
+
+                delete neighborPairs;
+                neighborPairs = NULL;
+                delete neighbors;
+                neighbors = NULL;
+                maxNeighborPairs = (int) (1.1*(*numPairs));
+                neighborPairs = OpenCLArray::create<mm_int2>(cl, maxNeighborPairs, "customManyParticleNeighborPairs");
+                neighbors = OpenCLArray::create<int>(cl, maxNeighborPairs, "customManyParticleNeighbors");
+                continue;
+            }
+        }
+        break;
+    }
+    return 0.0;
+}
+
+void OpenCLCalcCustomManyParticleForceKernel::copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force) {
+    int numParticles = force.getNumParticles();
+    if (numParticles != cl.getNumAtoms())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+    
+    // Record the per-particle parameters.
+    
+    vector<vector<float> > paramVector(numParticles);
+    vector<double> parameters;
+    int type;
+    for (int i = 0; i < numParticles; i++) {
+        force.getParticleParameters(i, parameters, type);
+        paramVector[i].resize(parameters.size());
+        for (int j = 0; j < (int) parameters.size(); j++)
+            paramVector[i][j] = (float) parameters[j];
+    }
+    params->setParameterValues(paramVector);
+    
+    // Mark that the current reordering may be invalid.
+    
+    cl.invalidateMolecules();
+}
+
 OpenCLIntegrateVerletStepKernel::~OpenCLIntegrateVerletStepKernel() {
 }
 

platforms/opencl/src/OpenCLPlatform.cpp

@@ -75,6 +75,7 @@ OpenCLPlatform::OpenCLPlatform() {
     registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);

platforms/opencl/src/kernels/customManyParticle.cl

+#pragma OPENCL EXTENSION cl_khr_global_int32_base_atomics : enable
+#pragma OPENCL EXTENSION cl_khr_int64_base_atomics : enable
+
+/**
+ * Record the force on an atom to global memory.
+ */
+inline void storeForce(int atom, real4 force, __global long* restrict forceBuffers) {
+    atom_add(&forceBuffers[atom], (long) (force.x*0x100000000));
+    atom_add(&forceBuffers[atom+PADDED_NUM_ATOMS], (long) (force.y*0x100000000));
+    atom_add(&forceBuffers[atom+2*PADDED_NUM_ATOMS], (long) (force.z*0x100000000));
+}
+
+/**
+ * Compute the difference between two vectors, taking periodic boundary conditions into account
+ * and setting the fourth component to the squared magnitude.
+ */
+inline real4 delta(real4 vec1, real4 vec2, real4 periodicBoxSize, real4 invPeriodicBoxSize) {
+    real4 result = (real4) (vec1.x-vec2.x, vec1.y-vec2.y, vec1.z-vec2.z, 0.0f);
+#ifdef USE_PERIODIC
+    result.x -= floor(result.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+    result.y -= floor(result.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+    result.z -= floor(result.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
+#endif
+    result.w = result.x*result.x + result.y*result.y + result.z*result.z;
+    return result;
+}
+
+/**
+ * Compute the angle between two vectors.  The w component of each vector should contain the squared magnitude.
+ */
+real computeAngle(real4 vec1, real4 vec2) {
+    real dotProduct = vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
+    real cosine = dotProduct*RSQRT(vec1.w*vec2.w);
+    real angle;
+    if (cosine > 0.99f || cosine < -0.99f) {
+        // We're close to the singularity in acos(), so take the cross product and use asin() instead.
+
+        real4 crossProduct = cross(vec1, vec2);
+        real scale = vec1.w*vec2.w;
+        angle = asin(SQRT(dot(crossProduct, crossProduct)/scale));
+        if (cosine < 0.0f)
+            angle = M_PI-angle;
+    }
+    else
+       angle = acos(cosine);
+    return angle;
+}
+
+/**
+ * Compute the cross product of two vectors, setting the fourth component to the squared magnitude.
+ */
+inline real4 computeCross(real4 vec1, real4 vec2) {
+    real4 cp = cross(vec1, vec2);
+    return (real4) (cp.x, cp.y, cp.z, cp.x*cp.x+cp.y*cp.y+cp.z*cp.z);
+}
+
+/**
+ * Determine whether a particular interaction is in the list of exclusions.
+ */
+inline bool isInteractionExcluded(int atom1, int atom2, __global int* restrict exclusions, __global int* restrict exclusionStartIndex) {
+    int first = exclusionStartIndex[atom1];
+    int last = exclusionStartIndex[atom1+1];
+    for (int i = last-1; i >= first; i--) {
+        int excluded = exclusions[i];
+        if (excluded == atom2)
+            return true;
+        if (excluded <= atom1)
+            return false;
+    }
+    return false;
+}
+
+/**
+ * Compute the interaction.
+ */
+__kernel void computeInteraction(
+        __global long* restrict forceBuffers, __global real* restrict energyBuffer, __global const real4* restrict posq,
+        real4 periodicBoxSize, real4 invPeriodicBoxSize
+#ifdef USE_CUTOFF
+        , __global const int* restrict neighbors, __global const int* restrict neighborStartIndex
+#endif
+#ifdef USE_FILTERS
+        , __global int* restrict particleTypes, __global int* restrict orderIndex, __global int* restrict particleOrder
+#endif
+#ifdef USE_EXCLUSIONS
+        , __global int* restrict exclusions, __global int* restrict exclusionStartIndex
+#endif
+        PARAMETER_ARGUMENTS) {
+    real energy = 0.0f;
+    
+    // Loop over particles to be the first one in the set.
+    
+    for (int p1 = get_group_id(0); p1 < NUM_ATOMS; p1 += get_num_groups(0)) {
+#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 = get_local_id(0); index < numCombinations; index += get_local_size(0)) {
+            FIND_ATOMS_FOR_COMBINATION_INDEX;
+            bool includeInteraction = IS_VALID_COMBINATION;
+#ifdef USE_CUTOFF
+            if (includeInteraction) {
+                VERIFY_CUTOFF;
+            }
+#endif
+#ifdef USE_FILTERS
+            int order = orderIndex[COMPUTE_TYPE_INDEX];
+            if (order == -1)
+                includeInteraction = false;
+#endif
+#ifdef USE_EXCLUSIONS
+            if (includeInteraction) {
+                VERIFY_EXCLUSIONS;
+            }
+#endif
+            if (includeInteraction) {
+                PERMUTE_ATOMS;
+                LOAD_PARTICLE_DATA;
+                COMPUTE_INTERACTION;
+            }
+        }
+    }
+    energyBuffer[get_global_id(0)] += energy;
+}
+
+/**
+ * Find a bounding box for the atoms in each block.
+ */
+__kernel void findBlockBounds(real4 periodicBoxSize, real4 invPeriodicBoxSize, __global const real4* restrict posq,
+        __global real4* restrict blockCenter, __global real4* restrict blockBoundingBox, __global int* restrict numNeighborPairs) {
+    int index = get_global_id(0);
+    int base = index*TILE_SIZE;
+    while (base < NUM_ATOMS) {
+        real4 pos = posq[base];
+#ifdef USE_PERIODIC
+        pos.x -= floor(pos.x*invPeriodicBoxSize.x)*periodicBoxSize.x;
+        pos.y -= floor(pos.y*invPeriodicBoxSize.y)*periodicBoxSize.y;
+        pos.z -= floor(pos.z*invPeriodicBoxSize.z)*periodicBoxSize.z;
+#endif
+        real4 minPos = pos;
+        real4 maxPos = pos;
+        int last = min(base+TILE_SIZE, NUM_ATOMS);
+        for (int i = base+1; i < last; i++) {
+            pos = posq[i];
+#ifdef USE_PERIODIC
+            real4 center = 0.5f*(maxPos+minPos);
+            pos.x -= floor((pos.x-center.x)*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+            pos.y -= floor((pos.y-center.y)*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+            pos.z -= floor((pos.z-center.z)*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
+#endif
+            minPos = (real4) (min(minPos.x,pos.x), min(minPos.y,pos.y), min(minPos.z,pos.z), 0);
+            maxPos = (real4) (max(maxPos.x,pos.x), max(maxPos.y,pos.y), max(maxPos.z,pos.z), 0);
+        }
+        real4 blockSize = 0.5f*(maxPos-minPos);
+        blockBoundingBox[index] = blockSize;
+        blockCenter[index] = 0.5f*(maxPos+minPos);
+        index += get_global_size(0);
+        base = index*TILE_SIZE;
+    }
+    if (get_group_id(0) == 0 && get_local_id(0) == 0)
+        *numNeighborPairs = 0;
+}
+
+/**
+ * Find a list of neighbors for each atom.
+ */
+__kernel void findNeighbors(real4 periodicBoxSize, real4 invPeriodicBoxSize, __global const real4* restrict posq,
+        __global const real4* restrict blockCenter, __global const real4* restrict blockBoundingBox, __global int2* restrict neighborPairs,
+        __global int* restrict numNeighborPairs, __global int* restrict numNeighborsForAtom, int maxNeighborPairs
+#ifdef USE_EXCLUSIONS
+        , __global int* restrict exclusions, __global int* restrict exclusionStartIndex
+#endif
+        ) {
+    __local real4 positionCache[FIND_NEIGHBORS_WORKGROUP_SIZE];
+    __local bool includeBlockFlags[FIND_NEIGHBORS_WORKGROUP_SIZE];
+    int indexInWarp = get_local_id(0)%32;
+    int warpStart = get_local_id(0)-indexInWarp;
+    for (int atom1 = get_global_id(0); atom1 < PADDED_NUM_ATOMS; atom1 += get_global_size(0)) {
+        // Load data for this atom.  Note that all threads in a warp are processing atoms from the same block.
+        
+        real4 pos1 = posq[atom1];
+        int block1 = atom1/TILE_SIZE;
+        real4 blockCenter1 = blockCenter[block1];
+        real4 blockSize1 = blockBoundingBox[block1];
+        int totalNeighborsForAtom1 = 0;
+        
+        // Loop over atom blocks to search for neighbors.  The threads in a warp compare block1 against 32
+        // other blocks in parallel.
+
+#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) {
+                real4 blockCenter2 = blockCenter[block2];
+                real4 blockSize2 = blockBoundingBox[block2];
+                real4 blockDelta = blockCenter1-blockCenter2;
+#ifdef USE_PERIODIC
+                blockDelta.x -= floor(blockDelta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                blockDelta.y -= floor(blockDelta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                blockDelta.z -= floor(blockDelta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
+#endif
+                blockDelta.x = max((real) 0, fabs(blockDelta.x)-blockSize1.x-blockSize2.x);
+                blockDelta.y = max((real) 0, fabs(blockDelta.y)-blockSize1.y-blockSize2.y);
+                blockDelta.z = max((real) 0, fabs(blockDelta.z)-blockSize1.z-blockSize2.z);
+                includeBlock2 &= (blockDelta.x*blockDelta.x+blockDelta.y*blockDelta.y+blockDelta.z*blockDelta.z < CUTOFF_SQUARED);
+            }
+            
+            // Loop over any blocks we identified as potentially containing neighbors.
+            
+            includeBlockFlags[get_local_id(0)] = includeBlock2;
+            SYNC_WARPS;
+            for (int i = 0; i < TILE_SIZE; i++) {
+                if (includeBlockFlags[warpStart+i]) {
+                    int block2 = block2Base+i;
+
+                    // Loop over atoms in this block.
+
+                    int start = block2*TILE_SIZE;
+                    int included[TILE_SIZE];
+                    int numIncluded = 0;
+                    positionCache[get_local_id(0)] = posq[start+indexInWarp];
+                    if (atom1 < NUM_ATOMS) {
+                        for (int j = 0; j < 32; j++) {
+                            int atom2 = start+j;
+                            real4 pos2 = positionCache[get_local_id(0)-indexInWarp+j];
+
+                            // 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);
+#endif
+                            if (includeAtom)
+                                included[numIncluded++] = atom2;
+                        }
+                    }
+
+                    // If we found any neighbors, store them to the neighbor list.
+
+                    if (numIncluded > 0) {
+                        int baseIndex = atom_add(numNeighborPairs, numIncluded);
+                        if (baseIndex+numIncluded <= maxNeighborPairs)
+                            for (int j = 0; j < numIncluded; j++)
+                                neighborPairs[baseIndex+j] = (int2) (atom1, included[j]);
+                        totalNeighborsForAtom1 += numIncluded;
+                    }
+                }
+            }
+        }
+        numNeighborsForAtom[atom1] = totalNeighborsForAtom1;
+    }
+}
+
+/**
+ * Sum the neighbor counts to compute the start position of each atom.  This kernel
+ * is executed as a single work group.
+ */
+__kernel void computeNeighborStartIndices(__global int* restrict numNeighborsForAtom, __global int* restrict neighborStartIndex,
+            __global int* restrict numNeighborPairs, int maxNeighborPairs) {
+    __local unsigned int posBuffer[256];
+    if (*numNeighborPairs > maxNeighborPairs) {
+        // There wasn't enough memory for the neighbor list, so we'll need to rebuild it.  Set the neighbor start
+        // indices to indicate no neighbors for any atom.
+        
+        for (int i = get_local_id(0); i <= NUM_ATOMS; i += get_local_size(0))
+            neighborStartIndex[i] = 0;
+        return;
+    }
+    unsigned int globalOffset = 0;
+    for (unsigned int startAtom = 0; startAtom < NUM_ATOMS; startAtom += get_local_size(0)) {
+        // Load the neighbor counts into local memory.
+
+        unsigned int globalIndex = startAtom+get_local_id(0);
+        posBuffer[get_local_id(0)] = (globalIndex < NUM_ATOMS ? numNeighborsForAtom[globalIndex] : 0);
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+        // Perform a parallel prefix sum.
+
+        for (unsigned int step = 1; step < get_local_size(0); step *= 2) {
+            unsigned int add = (get_local_id(0) >= step ? posBuffer[get_local_id(0)-step] : 0);
+            barrier(CLK_LOCAL_MEM_FENCE);
+            posBuffer[get_local_id(0)] += add;
+            barrier(CLK_LOCAL_MEM_FENCE);
+        }
+
+        // Write the results back to global memory.
+
+        if (globalIndex < NUM_ATOMS) {
+            neighborStartIndex[globalIndex+1] = posBuffer[get_local_id(0)]+globalOffset;
+            numNeighborsForAtom[globalIndex] = 0; // Clear this so the next kernel can use it as a counter
+        }
+        globalOffset += posBuffer[get_local_size(0)-1];
+    }
+    if (get_local_id(0) == 0)
+        neighborStartIndex[0] = 0;
+}
+
+/**
+ * Assemble the final neighbor list.
+ */
+__kernel void copyPairsToNeighborList(__global const int2* restrict neighborPairs, __global int* restrict neighbors, __global int* restrict numNeighborPairs,
+            int maxNeighborPairs, __global int* restrict numNeighborsForAtom, __global const int* restrict neighborStartIndex) {
+    int actualPairs = *numNeighborPairs;
+    if (actualPairs > maxNeighborPairs)
+        return; // There wasn't enough memory for the neighbor list, so we'll need to rebuild it.
+    for (unsigned int index = get_global_id(0); index < actualPairs; index += get_global_size(0)) {
+        int2 pair = neighborPairs[index];
+        int startIndex = neighborStartIndex[pair.x];
+        int offset = atom_add(numNeighborsForAtom+pair.x, 1);
+        neighbors[startIndex+offset] = pair.y;
+    }
+}

platforms/opencl/tests/TestOpenCLCustomManyParticleForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2014 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the OpenCL implementation of CustomManyParticleForce.
+ */
+
+#ifdef WIN32
+  #define _USE_MATH_DEFINES // Needed to get M_PI
+#endif
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "OpenCLPlatform.h"
+#include "openmm/CustomCompoundBondForce.h"
+#include "openmm/CustomManyParticleForce.h"
+#include "openmm/System.h"
+#include "openmm/TabulatedFunction.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+OpenCLPlatform platform;
+
+void validateAxilrodTeller(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 c = context.getParameter("C");
+    
+    // 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);
+        double rprod = r12*r13*r23;
+        expectedEnergy += c*(1+3*ctheta1*ctheta2*ctheta3)/(rprod*rprod*rprod);
+    }
+    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 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;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    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[4][3] = {{0,1,2}, {1,2,3}, {2,3,0}, {3,0,1}};
+    vector<const int*> expectedSets(&sets[0], &sets[4]);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+}
+
+void testCutoff() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffNonPeriodic);
+    force->setCutoffDistance(1.55);
+    vector<double> params;
+    force->addParticle(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));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    int sets[7][3] = {{0,1,2}, {0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,2,4}, {1,3,4}};
+    vector<const int*> expectedSets(&sets[0], &sets[7]);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+}
+
+void testPeriodic() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
+    force->setCutoffDistance(1.05);
+    vector<double> params;
+    force->addParticle(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));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    double boxSize = 2.1;
+    int sets[5][3] = {{0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,3,4}};
+    vector<const int*> expectedSets(&sets[0], &sets[5]);
+    validateAxilrodTeller(force, positions, expectedSets, boxSize);
+}
+
+void testExclusions() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    vector<double> params;
+    force->addParticle(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));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    force->addExclusion(0, 2);
+    force->addExclusion(0, 3);
+    int sets[5][3] = {{0,1,4}, {1,2,3}, {1,2,4}, {1,3,4}, {2,3,4}};
+    vector<const int*> expectedSets(&sets[0], &sets[5]);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+}
+
+void testAllTerms() {
+    int numParticles = 4;
+    
+    // Create a system with a CustomManyParticleForce.
+    
+    System system1;
+    CustomManyParticleForce* force1 = new CustomManyParticleForce(4,
+        "distance(p1,p2)+angle(p1,p4,p3)+dihedral(p1,p3,p2,p4)+x1+y4+z3");
+    system1.addForce(force1);
+    vector<double> params;
+    for (int i = 0; i < numParticles; i++) {
+        system1.addParticle(1.0);
+        force1->addParticle(params, i);
+    }
+    set<int> filter;
+    filter.insert(0);
+    force1->setTypeFilter(0, filter);
+    filter.clear();
+    filter.insert(1);
+    force1->setTypeFilter(1, filter);
+    filter.clear();
+    filter.insert(3);
+    force1->setTypeFilter(2, filter);
+    filter.clear();
+    filter.insert(2);
+    force1->setTypeFilter(3, filter);
+    
+    // Create a system that use a CustomCompoundBondForce to compute exactly the same interactions.
+    
+    System system2;
+    CustomCompoundBondForce* force2 = new CustomCompoundBondForce(4,
+        "distance(p1,p2)+angle(p1,p3,p4)+dihedral(p1,p4,p2,p3)+x1+y3+z4");
+    system2.addForce(force2);
+    vector<int> particles;
+    particles.push_back(0);
+    particles.push_back(1);
+    particles.push_back(2);
+    particles.push_back(3);
+    force2->addBond(particles, params);
+    for (int i = 0; i < numParticles; i++)
+        system2.addParticle(1.0);
+    
+    // Create contexts for both of them.
+
+    vector<Vec3> positions;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; i++)
+        positions.push_back(Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
+    VerletIntegrator integrator1(0.001);
+    VerletIntegrator integrator2(0.001);
+    Context context1(system1, integrator1, platform);
+    Context context2(system2, integrator2, platform);
+    context1.setPositions(positions);
+    context2.setPositions(positions);
+    
+    // See if they produce identical forces and energies.
+    
+    State state1 = context1.getState(State::Forces | State::Energy);
+    State state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state2.getPotentialEnergy(), state1.getPotentialEnergy(), 1e-4);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state2.getForces()[i], state1.getForces()[i], 1e-4);
+}
+
+void testParameters() {
+    // Create a system.
+    
+    int numParticles = 5;
+    System system;
+    CustomManyParticleForce* force = new CustomManyParticleForce(3, "C*scale1*scale2*scale3*(distance(p1,p2)+distance(p2,p3)+distance(p1,p3))");
+    force->addGlobalParameter("C", 2.0);
+    force->addPerParticleParameter("scale");
+    vector<double> params(1);
+    vector<Vec3> positions;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; i++) {
+        params[0] = i+1;
+        force->addParticle(params);
+        positions.push_back(Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
+        system.addParticle(1.0);
+    }
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    
+    // See if the energy is correct.
+
+    State state = context.getState(State::Energy);
+    double expectedEnergy = 0;
+    for (int i = 0; i < numParticles; i++)
+        for (int j = i+1; j < numParticles; j++)
+            for (int k = j+1; k < numParticles; k++) {
+                Vec3 d12 = positions[j]-positions[i];
+                Vec3 d13 = positions[k]-positions[i];
+                Vec3 d23 = positions[k]-positions[j];
+                double r12 = sqrt(d12.dot(d12));
+                double r13 = sqrt(d13.dot(d13));
+                double r23 = sqrt(d23.dot(d23));
+                expectedEnergy += 2.0*(i+1)*(j+1)*(k+1)*(r12+r13+r23);
+            }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+    
+    // Modify the parameters.
+    
+    context.setParameter("C", 3.5);
+    for (int i = 0; i < numParticles; i++) {
+        params[0] = 0.5*i-0.1;
+        force->setParticleParameters(i, params, 0);
+    }
+    force->updateParametersInContext(context);
+    
+    // See if the energy is still correct.
+    
+    state = context.getState(State::Energy);
+    expectedEnergy = 0;
+    for (int i = 0; i < numParticles; i++)
+        for (int j = i+1; j < numParticles; j++)
+            for (int k = j+1; k < numParticles; k++) {
+                Vec3 d12 = positions[j]-positions[i];
+                Vec3 d13 = positions[k]-positions[i];
+                Vec3 d23 = positions[k]-positions[j];
+                double r12 = sqrt(d12.dot(d12));
+                double r13 = sqrt(d13.dot(d13));
+                double r23 = sqrt(d23.dot(d23));
+                expectedEnergy += 3.5*(0.5*i-0.1)*(0.5*j-0.1)*(0.5*k-0.1)*(r12+r13+r23);
+            }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+}
+
+void testTabulatedFunctions() {
+    int numParticles = 5;
+    
+    // Create two tabulated functions.
+    
+    vector<double> values;
+    values.push_back(0.0);
+    values.push_back(50.0);
+    Continuous1DFunction* f1 = new Continuous1DFunction(values, 0, 100);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<double> c(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        c[i] = genrand_real2(sfmt);
+    values.resize(numParticles*numParticles*numParticles);
+    for (int i = 0; i < numParticles; i++)
+        for (int j = 0; j < numParticles; j++)
+            for (int k = 0; k < numParticles; k++)
+                values[i+numParticles*j+numParticles*numParticles*k] = c[i]+c[j]+c[k];
+    Discrete3DFunction* f2 = new Discrete3DFunction(numParticles, numParticles, numParticles, values);
+    
+    // Create a system.
+    
+    System system;
+    CustomManyParticleForce* force = new CustomManyParticleForce(3, "f1(distance(p1,p2)+distance(p2,p3)+distance(p1,p3))*f2(atom1, atom2, atom3)");
+    force->addPerParticleParameter("atom");
+    force->addTabulatedFunction("f1", f1);
+    force->addTabulatedFunction("f2", f2);
+    vector<double> params(1);
+    vector<Vec3> positions;
+    for (int i = 0; i < numParticles; i++) {
+        params[0] = i;
+        force->addParticle(params);
+        positions.push_back(Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
+        system.addParticle(1.0);
+    }
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    
+    // See if the energy is correct.
+
+    State state = context.getState(State::Energy);
+    double expectedEnergy = 0;
+    for (int i = 0; i < numParticles; i++)
+        for (int j = i+1; j < numParticles; j++)
+            for (int k = j+1; k < numParticles; k++) {
+                Vec3 d12 = positions[j]-positions[i];
+                Vec3 d13 = positions[k]-positions[i];
+                Vec3 d23 = positions[k]-positions[j];
+                double r12 = sqrt(d12.dot(d12));
+                double r13 = sqrt(d13.dot(d13));
+                double r23 = sqrt(d23.dot(d23));
+                expectedEnergy += 0.5*(r12+r13+r23)*(c[i]+c[j]+c[k]);
+            }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+}
+
+void testTypeFilters() {
+    // Create a system.
+    
+    System system;
+    for (int i = 0; i < 5; i++)
+        system.addParticle(1.0);
+    CustomManyParticleForce* force = new CustomManyParticleForce(3, "c1*(distance(p1,p2)+distance(p1,p3))");
+    force->addPerParticleParameter("c");
+    double c[] = {1.0, 2.0, 1.3, 1.5, -2.1};
+    int type[] = {0, 1, 0, 1, 5};
+    vector<double> params(1);
+    for (int i = 0; i < 5; i++) {
+        params[0] = c[i];
+        force->addParticle(params, type[i]);
+    }
+    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));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    set<int> f1, f2;
+    f1.insert(0);
+    f2.insert(1);
+    f2.insert(5);
+    force->setTypeFilter(0, f1);
+    force->setTypeFilter(1, f2);
+    force->setTypeFilter(2, f2);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    
+    // See if the energy is correct.
+
+    State state = context.getState(State::Energy);
+    double expectedEnergy = 0;
+    int sets[6][3] = {{0,1,3}, {0,1,4}, {0,3,4}, {2,1,3}, {2,1,4}, {2,3,4}};
+    for (int i = 0; i < 6; i++) {
+        int p1 = sets[i][0];
+        int p2 = sets[i][1];
+        int p3 = sets[i][2];
+            Vec3 d12 = positions[p2]-positions[p1];
+            Vec3 d13 = positions[p3]-positions[p1];
+            double r12 = sqrt(d12.dot(d12));
+            double r13 = sqrt(d13.dot(d13));
+            expectedEnergy += c[p1]*(r12+r13);
+    }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+}
+
+void testLargeSystem() {
+    int gridSize = 8;
+    int numParticles = gridSize*gridSize*gridSize;
+    double boxSize = 3.0;
+    double spacing = boxSize/gridSize;
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
+    force->setCutoffDistance(0.6);
+    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);
+}
+
+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(0.1, 0, 0));
+    positions.push_back(Vec3(0, 0.11, 0.03));
+    positions.push_back(Vec3(0.04, 0, -0.08));
+    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(0.155);
+    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(0.1, 0, 0));
+    positions.push_back(Vec3(0, 0.11, 0.03));
+    positions.push_back(Vec3(0.04, 0, -0.08));
+    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)
+            platform.setPropertyDefaultValue("OpenCLPrecision", string(argv[1]));
+        testNoCutoff();
+        testCutoff();
+        testPeriodic();
+        testExclusions();
+        testAllTerms();
+        testParameters();
+        testTabulatedFunctions();
+        testTypeFilters();
+        testLargeSystem();
+        testCentralParticleModeNoCutoff();
+        testCentralParticleModeCutoff();
+        testCentralParticleModeLargeSystem();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/reference/include/ReferenceBondIxn.h

@@ -82,7 +82,7 @@ class OPENMM_EXPORT ReferenceBondIxn {
       
          --------------------------------------------------------------------------------------- */
       
-      RealOpenMM getNormedDotProduct( RealOpenMM* vector1, RealOpenMM* vector2, int hasREntry ) const;
+      static RealOpenMM getNormedDotProduct( RealOpenMM* vector1, RealOpenMM* vector2, int hasREntry );
       
       /**---------------------------------------------------------------------------------------
       
@@ -98,8 +98,8 @@ class OPENMM_EXPORT ReferenceBondIxn {
       
          --------------------------------------------------------------------------------------- */
       
-      RealOpenMM getAngleBetweenTwoVectors( RealOpenMM* vector1, RealOpenMM* vector2, 
-                                            RealOpenMM* outputDotProduct, int hasREntry ) const;
+      static RealOpenMM getAngleBetweenTwoVectors( RealOpenMM* vector1, RealOpenMM* vector2, 
+                                                   RealOpenMM* outputDotProduct, int hasREntry );
       
       /**---------------------------------------------------------------------------------------
       
@@ -119,10 +119,10 @@ class OPENMM_EXPORT ReferenceBondIxn {
       
          --------------------------------------------------------------------------------------- */
       
-      RealOpenMM getDihedralAngleBetweenThreeVectors( RealOpenMM* vector1, RealOpenMM* vector2, 
+      static RealOpenMM getDihedralAngleBetweenThreeVectors( RealOpenMM* vector1, RealOpenMM* vector2, 
                                                              RealOpenMM* vector3, RealOpenMM** outputCrossProduct, 
                                                              RealOpenMM* cosineOfAngle, RealOpenMM* signVector, 
-                                                      RealOpenMM* signOfAngle, int hasREntry ) const;
+                                                             RealOpenMM* signOfAngle, int hasREntry );
       
 };
 

platforms/reference/include/ReferenceCustomManyParticleIxn.h

+
+/* Portions copyright (c) 2009-2014 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef __ReferenceCustomManyParticleIxn_H__
+#define __ReferenceCustomManyParticleIxn_H__
+
+#include "ReferenceBondIxn.h"
+#include "openmm/CustomManyParticleForce.h"
+#include "lepton/ExpressionProgram.h"
+#include "lepton/ParsedExpression.h"
+#include <map>
+#include <set>
+#include <vector>
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceCustomManyParticleIxn {
+
+   private:
+
+      class ParticleTermInfo;
+      class DistanceTermInfo;
+      class AngleTermInfo;
+      class DihedralTermInfo;
+      int numParticlesPerSet, numPerParticleParameters, numTypes;
+      bool useCutoff, usePeriodic, centralParticleMode;
+      RealOpenMM cutoffDistance;
+      RealOpenMM periodicBoxSize[3];
+      Lepton::ExpressionProgram energyExpression;
+      std::vector<std::vector<std::string> > particleParamNames;
+      std::vector<std::set<int> > exclusions;
+      std::vector<int> particleTypes;
+      std::vector<int> orderIndex;
+      std::vector<std::vector<int> > particleOrder;
+      std::vector<ParticleTermInfo> particleTerms;
+      std::vector<DistanceTermInfo> distanceTerms;
+      std::vector<AngleTermInfo> angleTerms;
+      std::vector<DihedralTermInfo> dihedralTerms;
+
+      void loopOverInteractions(std::vector<int>& particles, int loopIndex, std::vector<OpenMM::RealVec>& atomCoordinates,
+                                RealOpenMM** particleParameters, std::map<std::string, double>& variables, std::vector<OpenMM::RealVec>& forces,
+                                RealOpenMM* totalEnergy) const;
+
+      /**---------------------------------------------------------------------------------------
+
+         Calculate custom interaction for one set of particles
+
+         @param particles          the indices of the particles
+         @param atomCoordinates    atom coordinates
+         @param particleParameters particle parameter values (particleParameters[particleIndex][parameterIndex])
+         @param variables          the values of variables that may appear in expressions
+         @param forces             force array (forces added)
+         @param totalEnergy        total energy
+
+         --------------------------------------------------------------------------------------- */
+
+      void calculateOneIxn(const std::vector<int>& particles, std::vector<OpenMM::RealVec>& atomCoordinates,
+                           RealOpenMM** particleParameters, std::map<std::string, double>& variables, std::vector<OpenMM::RealVec>& forces,
+                           RealOpenMM* totalEnergy) const;
+
+      void computeDelta(int atom1, int atom2, RealOpenMM* delta, std::vector<OpenMM::RealVec>& atomCoordinates) const;
+
+      static RealOpenMM computeAngle(RealOpenMM* vec1, RealOpenMM* vec2);
+
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+
+         Constructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceCustomManyParticleIxn(const OpenMM::CustomManyParticleForce& force);
+
+      /**---------------------------------------------------------------------------------------
+
+         Destructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceCustomManyParticleIxn();
+
+      /**---------------------------------------------------------------------------------------
+
+         Set the force to use a cutoff.
+
+         @param distance            the cutoff distance
+
+         --------------------------------------------------------------------------------------- */
+
+      void setUseCutoff(RealOpenMM distance);
+
+      /**---------------------------------------------------------------------------------------
+
+         Set the force to use periodic boundary conditions.  This requires that a cutoff has
+         already been set, and the smallest side of the periodic box is at least twice the cutoff
+         distance.
+
+         @param boxSize             the X, Y, and Z widths of the periodic box
+
+         --------------------------------------------------------------------------------------- */
+
+      void setPeriodic(OpenMM::RealVec& boxSize);
+
+      /**---------------------------------------------------------------------------------------
+
+         Calculate the interaction
+
+         @param atomCoordinates    atom coordinates
+         @param particleParameters particle parameter values (particleParameters[particleIndex][parameterIndex])
+         @param globalParameters   the values of global parameters
+         @param forces             force array (forces added)
+         @param totalEnergy        total energy
+
+         --------------------------------------------------------------------------------------- */
+
+      void calculateIxn(std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** particleParameters,
+                        const std::map<std::string, double>& globalParameters,
+                        std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy) const;
+
+// ---------------------------------------------------------------------------------------
+
+};
+
+class ReferenceCustomManyParticleIxn::ParticleTermInfo {
+public:
+    std::string name;
+    int atom, component;
+    Lepton::ExpressionProgram forceExpression;
+    ParticleTermInfo(const std::string& name, int atom, int component, const Lepton::ExpressionProgram& forceExpression) :
+            name(name), atom(atom), component(component), forceExpression(forceExpression) {
+    }
+};
+
+class ReferenceCustomManyParticleIxn::DistanceTermInfo {
+public:
+    std::string name;
+    int p1, p2;
+    Lepton::ExpressionProgram forceExpression;
+    mutable RealOpenMM delta[ReferenceForce::LastDeltaRIndex];
+    DistanceTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
+            name(name), p1(atoms[0]), p2(atoms[1]), forceExpression(forceExpression) {
+    }
+};
+
+class ReferenceCustomManyParticleIxn::AngleTermInfo {
+public:
+    std::string name;
+    int p1, p2, p3;
+    Lepton::ExpressionProgram forceExpression;
+    mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
+    mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
+    AngleTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
+            name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), forceExpression(forceExpression) {
+    }
+};
+
+class ReferenceCustomManyParticleIxn::DihedralTermInfo {
+public:
+    std::string name;
+    int p1, p2, p3, p4;
+    Lepton::ExpressionProgram forceExpression;
+    mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
+    mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
+    mutable RealOpenMM delta3[ReferenceForce::LastDeltaRIndex];
+    mutable RealOpenMM cross1[3];
+    mutable RealOpenMM cross2[3];
+    DihedralTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
+            name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), p4(atoms[3]), forceExpression(forceExpression) {
+    }
+};
+
+#endif // __ReferenceCustomManyParticleIxn_H__

platforms/reference/include/ReferenceKernels.h

@@ -44,6 +44,7 @@ class CpuGBVI;
 class ReferenceAndersenThermostat;
 class ReferenceCustomCompoundBondIxn;
 class ReferenceCustomHbondIxn;
+class ReferenceCustomManyParticleIxn;
 class ReferenceBrownianDynamics;
 class ReferenceStochasticDynamics;
 class ReferenceConstraintAlgorithm;
@@ -861,6 +862,46 @@ private:
     std::vector<std::string> globalParameterNames;
 };
 
+/**
+ * This kernel is invoked by CustomManyParticleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcCustomManyParticleForceKernel : public CalcCustomManyParticleForceKernel {
+public:
+    ReferenceCalcCustomManyParticleForceKernel(std::string name, const Platform& platform) : CalcCustomManyParticleForceKernel(name, platform), ixn(NULL) {
+    }
+    ~ReferenceCalcCustomManyParticleForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomManyParticleForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomManyParticleForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomManyParticleForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force);
+private:
+    int numParticles;
+    RealOpenMM cutoffDistance;
+    RealOpenMM **particleParamArray;
+    ReferenceCustomManyParticleIxn* ixn;
+    std::vector<std::string> globalParameterNames;
+    NonbondedMethod nonbondedMethod;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -80,6 +80,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcCustomHbondForceKernel(name, platform);
     if (name == CalcCustomCompoundBondForceKernel::Name())
         return new ReferenceCalcCustomCompoundBondForceKernel(name, platform);
+    if (name == CalcCustomManyParticleForceKernel::Name())
+        return new ReferenceCalcCustomManyParticleForceKernel(name, platform);
     if (name == IntegrateVerletStepKernel::Name())
         return new ReferenceIntegrateVerletStepKernel(name, platform, data);
     if (name == IntegrateLangevinStepKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -47,6 +47,7 @@
 #include "ReferenceCustomGBIxn.h"
 #include "ReferenceCustomHbondIxn.h"
 #include "ReferenceCustomNonbondedIxn.h"
+#include "ReferenceCustomManyParticleIxn.h"
 #include "ReferenceCustomTorsionIxn.h"
 #include "ReferenceHarmonicBondIxn.h"
 #include "ReferenceLJCoulomb14.h"
@@ -1615,6 +1616,68 @@ void ReferenceCalcCustomCompoundBondForceKernel::copyParametersToContext(Context
     }
 }
 
+ReferenceCalcCustomManyParticleForceKernel::~ReferenceCalcCustomManyParticleForceKernel() {
+    disposeRealArray(particleParamArray, numParticles);
+    if (ixn != NULL)
+        delete ixn;
+}
+
+void ReferenceCalcCustomManyParticleForceKernel::initialize(const System& system, const CustomManyParticleForce& force) {
+
+    // Build the arrays.
+
+    numParticles = system.getNumParticles();
+    int numParticleParameters = force.getNumPerParticleParameters();
+    particleParamArray = allocateRealArray(numParticles, numParticleParameters);
+    for (int i = 0; i < numParticles; ++i) {
+        vector<double> parameters;
+        int type;
+        force.getParticleParameters(i, parameters, type);
+        for (int j = 0; j < numParticleParameters; j++)
+            particleParamArray[i][j] = parameters[j];
+    }
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        globalParameterNames.push_back(force.getGlobalParameterName(i));
+    ixn = new ReferenceCustomManyParticleIxn(force);
+    nonbondedMethod = CalcCustomManyParticleForceKernel::NonbondedMethod(force.getNonbondedMethod());
+    cutoffDistance = force.getCutoffDistance();
+}
+
+double ReferenceCalcCustomManyParticleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    vector<RealVec>& posData = extractPositions(context);
+    vector<RealVec>& forceData = extractForces(context);
+    RealOpenMM energy = 0;
+    map<string, double> globalParameters;
+    for (int i = 0; i < (int) globalParameterNames.size(); i++)
+        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    if (nonbondedMethod == CutoffPeriodic) {
+        RealVec& box = extractBoxSize(context);
+        double minAllowedSize = 2*cutoffDistance;
+        if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize)
+            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
+        ixn->setPeriodic(box);
+    }
+    ixn->calculateIxn(posData, particleParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
+    return energy;
+}
+
+void ReferenceCalcCustomManyParticleForceKernel::copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force) {
+    if (numParticles != force.getNumParticles())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+
+    // Record the values.
+
+    int numParameters = force.getNumPerParticleParameters();
+    vector<double> params;
+    for (int i = 0; i < numParticles; ++i) {
+        vector<double> parameters;
+        int type;
+        force.getParticleParameters(i, parameters, type);
+        for (int j = 0; j < numParameters; j++)
+            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
+    }
+}
+
 ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
     if (dynamics)
         delete dynamics;

platforms/reference/src/ReferencePlatform.cpp

@@ -63,6 +63,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);