CPU version of CustomManyParticleForce is multithreaded

platforms/cpu/include/CpuCustomManyParticleForce.h

@@ -36,6 +36,7 @@
 #include "lepton/ParsedExpression.h"
 #include <map>
 #include <set>
+#include <utility>
 #include <vector>
 
 namespace OpenMM {
@@ -47,26 +48,36 @@ private:
     class DistanceTermInfo;
     class AngleTermInfo;
     class DihedralTermInfo;
+    class ComputeForceTask;
+    class ThreadData;
     int numParticlesPerSet, numPerParticleParameters, numTypes;
     bool useCutoff, usePeriodic;
     RealOpenMM cutoffDistance;
     RealOpenMM periodicBoxSize[3];
     CpuNeighborList* neighborList;
     ThreadPool& threads;
-    CompiledExpressionSet expressionSet;
-    Lepton::CompiledExpression energyExpression;
-    std::vector<std::vector<int> > particleParamIndices;
     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;
+    std::vector<ThreadData*> threadData;
+    // The following variables are used to make information accessible to the individual threads.
+    int numParticles;
+    float* posq;
+    RealVec const* atomCoordinates;
+    RealOpenMM** particleParameters;        
+    const std::map<std::string, double>* globalParameters;
+    std::vector<AlignedArray<float> >* threadForce;
+    bool includeForces, includeEnergy;
+    void* atomicCounter;
+
+    /**
+     * This routine contains the code executed by each thread.
+     */
+    void threadComputeForce(ThreadPool& threads, int threadIndex);
 
-    void loopOverInteractions(std::vector<int>& availableParticles, std::vector<int>& particleSet, int loopIndex, int startIndex, float* posq, std::vector<OpenMM::RealVec>& atomCoordinates,
-                              RealOpenMM** particleParameters, std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize);
+    void loopOverInteractions(std::vector<int>& availableParticles, std::vector<int>& particleSet, int loopIndex, int startIndex,
+                              RealOpenMM** particleParameters, float* forces, ThreadData& data, const fvec4& boxSize, const fvec4& invBoxSize);
 
     /**---------------------------------------------------------------------------------------
 
@@ -81,8 +92,8 @@ private:
 
        --------------------------------------------------------------------------------------- */
 
-    void calculateOneIxn(std::vector<int>& particleSet, float* posq, std::vector<OpenMM::RealVec>& atomCoordinates,
-                         RealOpenMM** particleParameters, std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize);
+    void calculateOneIxn(std::vector<int>& particleSet,
+                         RealOpenMM** particleParameters, float* forces, ThreadData& data, const fvec4& boxSize, const fvec4& invBoxSize);
 
     /**
      * Compute the displacement and squared distance between two points, optionally using
@@ -90,9 +101,9 @@ private:
      */
     void getDeltaR(const fvec4& posI, const fvec4& posJ, fvec4& deltaR, float& r2, const fvec4& boxSize, const fvec4& invBoxSize) const;
 
-    void computeDelta(int atom1, int atom2, RealOpenMM* delta, std::vector<OpenMM::RealVec>& atomCoordinates) const;
+    void computeDelta(int atom1, int atom2, RealOpenMM* delta, const OpenMM::RealVec* atomCoordinates) const;
 
-    static RealOpenMM computeAngle(RealOpenMM* vec1, RealOpenMM* vec2);
+    static RealOpenMM computeAngle(RealOpenMM* vec1, RealOpenMM* vec2, float sign);
 
 
 public:
@@ -150,10 +161,7 @@ public:
 
     void calculateIxn(AlignedArray<float>& posq, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** particleParameters,
                       const std::map<std::string, double>& globalParameters,
-                      std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy);
-
-// ---------------------------------------------------------------------------------------
-
+                      std::vector<AlignedArray<float> >& threadForce, bool includeForces, bool includeEnergy, double& energy);
 };
 
 class CpuCustomManyParticleForce::ParticleTermInfo {
@@ -161,10 +169,7 @@ public:
     std::string name;
     int atom, component, variableIndex;
     Lepton::CompiledExpression forceExpression;
-    ParticleTermInfo(const std::string& name, int atom, int component, const Lepton::CompiledExpression& forceExpression, CompiledExpressionSet& set) :
-            name(name), atom(atom), component(component), forceExpression(forceExpression) {
-        variableIndex = set.getVariableIndex(name);
-    }
+    ParticleTermInfo(const std::string& name, int atom, int component, const Lepton::CompiledExpression& forceExpression, ThreadData& data);
 };
 
 class CpuCustomManyParticleForce::DistanceTermInfo {
@@ -172,11 +177,9 @@ public:
     std::string name;
     int p1, p2, variableIndex;
     Lepton::CompiledExpression forceExpression;
-    mutable RealOpenMM delta[ReferenceForce::LastDeltaRIndex];
-    DistanceTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression, CompiledExpressionSet& set) :
-            name(name), p1(atoms[0]), p2(atoms[1]), forceExpression(forceExpression) {
-        variableIndex = set.getVariableIndex(name);
-    }
+    int delta;
+    float deltaSign;
+    DistanceTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression, ThreadData& data);
 };
 
 class CpuCustomManyParticleForce::AngleTermInfo {
@@ -184,12 +187,9 @@ public:
     std::string name;
     int p1, p2, p3, variableIndex;
     Lepton::CompiledExpression forceExpression;
-    mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
-    mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
-    AngleTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression, CompiledExpressionSet& set) :
-            name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), forceExpression(forceExpression) {
-        variableIndex = set.getVariableIndex(name);
-    }
+    int delta1, delta2;
+    float delta1Sign, delta2Sign;
+    AngleTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression, ThreadData& data);
 };
 
 class CpuCustomManyParticleForce::DihedralTermInfo {
@@ -197,15 +197,29 @@ public:
     std::string name;
     int p1, p2, p3, p4, variableIndex;
     Lepton::CompiledExpression forceExpression;
-    mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
-    mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
-    mutable RealOpenMM delta3[ReferenceForce::LastDeltaRIndex];
+    int delta1, delta2, delta3;
     mutable RealOpenMM cross1[3];
     mutable RealOpenMM cross2[3];
-    DihedralTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression, CompiledExpressionSet& set) :
-            name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), p4(atoms[3]), forceExpression(forceExpression) {
-        variableIndex = set.getVariableIndex(name);
-    }
+    DihedralTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression, ThreadData& data);
+};
+
+class CpuCustomManyParticleForce::ThreadData {
+public:
+    CompiledExpressionSet expressionSet;
+    Lepton::CompiledExpression energyExpression;
+    std::vector<std::vector<int> > particleParamIndices;
+    std::vector<std::pair<int, int> > deltaPairs;
+    std::vector<ParticleTermInfo> particleTerms;
+    std::vector<DistanceTermInfo> distanceTerms;
+    std::vector<AngleTermInfo> angleTerms;
+    std::vector<DihedralTermInfo> dihedralTerms;
+    double energy;
+    ThreadData(const CustomManyParticleForce& force, Lepton::ParsedExpression& energyExpr,
+            std::map<std::string, std::vector<int> >& distances, std::map<std::string, std::vector<int> >& angles, std::map<std::string, std::vector<int> >& dihedrals);
+    /**
+     * Request a pair of particles whose distance or displacement vector is needed in the computation.
+     */
+    void requestDeltaPair(int p1, int p2, int& pairIndex, float& pairSign, bool allowReversed);
 };
 
 } // namespace OpenMM

platforms/cpu/src/CpuKernels.cpp

@@ -870,8 +870,6 @@ void CpuCalcCustomManyParticleForceKernel::initialize(const System& system, cons
 
 double CpuCalcCustomManyParticleForceKernel::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]);
@@ -882,7 +880,8 @@ double CpuCalcCustomManyParticleForceKernel::execute(ContextImpl& context, bool
             throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
         ixn->setPeriodic(box);
     }
-    ixn->calculateIxn(data.posq, posData, particleParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
+    double energy = 0;
+    ixn->calculateIxn(data.posq, posData, particleParamArray, globalParameters, data.threadForce, includeForces, includeEnergy, energy);
     return energy;
 }
 

platforms/reference/src/SimTKReference/ReferenceCustomManyParticleIxn.cpp

@@ -192,7 +192,7 @@ void ReferenceCustomManyParticleIxn::calculateOneIxn(const vector<int>& particle
 
     for (int i = 0; i < (int) particleTerms.size(); i++) {
         const ParticleTermInfo& term = particleTerms[i];
-        variables[term.name] = atomCoordinates[term.atom][term.component];
+        variables[term.name] = atomCoordinates[permutedParticles[term.atom]][term.component];
     }
     for (int i = 0; i < (int) distanceTerms.size(); i++) {
         const DistanceTermInfo& term = distanceTerms[i];

platforms/reference/tests/TestReferenceCustomManyParticleForce.cpp

@@ -39,6 +39,7 @@
 #include "openmm/internal/AssertionUtilities.h"
 #include "openmm/Context.h"
 #include "ReferencePlatform.h"
+#include "openmm/CustomCompoundBondForce.h"
 #include "openmm/CustomManyParticleForce.h"
 #include "openmm/System.h"
 #include "openmm/TabulatedFunction.h"
@@ -217,6 +218,72 @@ void testExclusions() {
     validateAxilrodTeller(force, positions, expectedSets, 2.0);
 }
 
+void testAllTerms() {
+    int numParticles = 4;
+    ReferencePlatform 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.
     
@@ -403,6 +470,7 @@ int main() {
         testCutoff();
         testPeriodic();
         testExclusions();
+        testAllTerms();
         testParameters();
         testTabulatedFunctions();
         testTypeFilters();