Continuing to create reference implementation of CustomManyParticleForce

openmmapi/src/CustomManyParticleForce.cpp

@@ -44,7 +44,7 @@ using namespace OpenMM;
 using namespace std;
 
 CustomManyParticleForce::CustomManyParticleForce(int particlesPerSet, const string& energy) :
-        particlesPerSet(particlesPerSet), energyExpression(energy), typeFilters(particlesPerSet) {
+        particlesPerSet(particlesPerSet), energyExpression(energy), nonbondedMethod(NoCutoff), typeFilters(particlesPerSet) {
 }
 
 CustomManyParticleForce::~CustomManyParticleForce() {

platforms/reference/include/ReferenceCustomManyParticleIxn.h

@@ -26,9 +26,11 @@
 #define __ReferenceCustomManyParticleIxn_H__
 
 #include "ReferenceBondIxn.h"
+#include "openmm/CustomManyParticleForce.h"
 #include "lepton/ExpressionProgram.h"
 #include "lepton/ParsedExpression.h"
 #include <map>
+#include <set>
 #include <vector>
 
 // ---------------------------------------------------------------------------------------
@@ -47,6 +49,7 @@ class ReferenceCustomManyParticleIxn {
       RealOpenMM periodicBoxSize[3];
       Lepton::ExpressionProgram energyExpression;
       std::vector<std::vector<std::string> > particleParamNames;
+      std::vector<std::set<int> > exclusions;
       std::vector<ParticleTermInfo> particleTerms;
       std::vector<DistanceTermInfo> distanceTerms;
       std::vector<AngleTermInfo> angleTerms;
@@ -85,9 +88,7 @@ class ReferenceCustomManyParticleIxn {
 
          --------------------------------------------------------------------------------------- */
 
-       ReferenceCustomManyParticleIxn(int numParticlesPerSet, const Lepton::ParsedExpression& energyExpression,
-                               const std::vector<std::string>& particleParameterNames, const std::map<std::string, std::vector<int> >& distances,
-                               const std::map<std::string, std::vector<int> >& angles, const std::map<std::string, std::vector<int> >& dihedrals);
+       ReferenceCustomManyParticleIxn(const OpenMM::CustomManyParticleForce& force);
 
       /**---------------------------------------------------------------------------------------
 

platforms/reference/src/ReferenceKernels.cpp

@@ -68,7 +68,6 @@
 #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/CMAPTorsionForceImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
@@ -1637,34 +1636,11 @@ void ReferenceCalcCustomManyParticleForceKernel::initialize(const System& system
         for (int j = 0; j < numParticleParameters; j++)
             particleParamArray[i][j] = parameters[j];
     }
-
-    // Create custom functions for the tabulated functions.
-
-    map<string, Lepton::CustomFunction*> functions;
-    for (int i = 0; i < force.getNumTabulatedFunctions(); i++)
-        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
-
-    // Parse the expression and create the object used to calculate the interaction.
-
-    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);
-    vector<string> particleParameterNames;
-    for (int i = 0; i < numParticleParameters; i++)
-        particleParameterNames.push_back(force.getPerParticleParameterName(i));
     for (int i = 0; i < force.getNumGlobalParameters(); i++)
         globalParameterNames.push_back(force.getGlobalParameterName(i));
-    ixn = new ReferenceCustomManyParticleIxn(force.getNumParticlesPerSet(), energyExpression, particleParameterNames, distances, angles, dihedrals);
+    ixn = new ReferenceCustomManyParticleIxn(force);
     nonbondedMethod = CalcCustomManyParticleForceKernel::NonbondedMethod(force.getNonbondedMethod());
     cutoffDistance = force.getCutoffDistance();
-    if (nonbondedMethod != NoCutoff)
-        ixn->setUseCutoff(cutoffDistance);
-
-    // Delete the custom functions.
-
-    for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
-        delete iter->second;
 }
 
 double ReferenceCalcCustomManyParticleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {

platforms/reference/src/SimTKReference/ReferenceCustomManyParticleIxn.cpp

@@ -31,40 +31,72 @@
 #include "SimTKOpenMMUtilities.h"
 #include "ReferenceForce.h"
 #include "ReferenceCustomManyParticleIxn.h"
+#include "ReferenceTabulatedFunction.h"
+#include "openmm/internal/CustomManyParticleForceImpl.h"
+#include "lepton/CustomFunction.h"
+
+using namespace OpenMM;
+using namespace std;
+
+ReferenceCustomManyParticleIxn::ReferenceCustomManyParticleIxn(const CustomManyParticleForce& force) : useCutoff(false), usePeriodic(false) {
+    numParticlesPerSet = force.getNumParticlesPerSet();
+    numPerParticleParameters = force.getNumPerParticleParameters();
+    
+    // Create custom functions for the tabulated functions.
+
+    map<string, Lepton::CustomFunction*> functions;
+    for (int i = 0; i < (int) force.getNumTabulatedFunctions(); i++)
+        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
+
+    // Parse the expression and create the object used to calculate the interaction.
+
+    map<string, vector<int> > distances;
+    map<string, vector<int> > angles;
+    map<string, vector<int> > dihedrals;
+    Lepton::ParsedExpression energyExpr = CustomManyParticleForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
+    energyExpression = energyExpr.createProgram();
+    vector<string> particleParameterNames;
+    if (force.getNonbondedMethod() != CustomManyParticleForce::NoCutoff)
+        setUseCutoff(force.getCutoffDistance());
+
+    // Delete the custom functions.
+
+    for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
+        delete iter->second;
+
+    // Differentiate the energy to get expressions for the force.
 
-using std::map;
-using std::pair;
-using std::string;
-using std::stringstream;
-using std::vector;
-using OpenMM::RealVec;
-
-ReferenceCustomManyParticleIxn::ReferenceCustomManyParticleIxn(int numParticlesPerSet,
-            const Lepton::ParsedExpression& energyExpression, const vector<string>& particleParameterNames,
-            const map<string, vector<int> >& distances, const map<string, vector<int> >& angles, const map<string, vector<int> >& dihedrals) :
-            numParticlesPerSet(numParticlesPerSet), energyExpression(energyExpression.createProgram()), useCutoff(false), usePeriodic(false) {
     particleParamNames.resize(numParticlesPerSet);
-    numPerParticleParameters = particleParameterNames.size();
     for (int i = 0; i < numParticlesPerSet; i++) {
         stringstream xname, yname, zname;
         xname << 'x' << (i+1);
         yname << 'y' << (i+1);
         zname << 'z' << (i+1);
-        particleTerms.push_back(ReferenceCustomManyParticleIxn::ParticleTermInfo(xname.str(), i, 0, energyExpression.differentiate(xname.str()).optimize().createProgram()));
-        particleTerms.push_back(ReferenceCustomManyParticleIxn::ParticleTermInfo(yname.str(), i, 1, energyExpression.differentiate(yname.str()).optimize().createProgram()));
-        particleTerms.push_back(ReferenceCustomManyParticleIxn::ParticleTermInfo(zname.str(), i, 2, energyExpression.differentiate(zname.str()).optimize().createProgram()));
+        particleTerms.push_back(ReferenceCustomManyParticleIxn::ParticleTermInfo(xname.str(), i, 0, energyExpr.differentiate(xname.str()).optimize().createProgram()));
+        particleTerms.push_back(ReferenceCustomManyParticleIxn::ParticleTermInfo(yname.str(), i, 1, energyExpr.differentiate(yname.str()).optimize().createProgram()));
+        particleTerms.push_back(ReferenceCustomManyParticleIxn::ParticleTermInfo(zname.str(), i, 2, energyExpr.differentiate(zname.str()).optimize().createProgram()));
         for (int j = 0; j < numPerParticleParameters; j++) {
             stringstream paramname;
-            paramname << particleParameterNames[j] << (i+1);
+            paramname << force.getPerParticleParameterName(j) << (i+1);
             particleParamNames[i].push_back(paramname.str());
         }
     }
     for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter)
-        distanceTerms.push_back(ReferenceCustomManyParticleIxn::DistanceTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
+        distanceTerms.push_back(ReferenceCustomManyParticleIxn::DistanceTermInfo(iter->first, iter->second, energyExpr.differentiate(iter->first).optimize().createProgram()));
     for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter)
-        angleTerms.push_back(ReferenceCustomManyParticleIxn::AngleTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
+        angleTerms.push_back(ReferenceCustomManyParticleIxn::AngleTermInfo(iter->first, iter->second, energyExpr.differentiate(iter->first).optimize().createProgram()));
     for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter)
-        dihedralTerms.push_back(ReferenceCustomManyParticleIxn::DihedralTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
+        dihedralTerms.push_back(ReferenceCustomManyParticleIxn::DihedralTermInfo(iter->first, iter->second, energyExpr.differentiate(iter->first).optimize().createProgram()));
+    
+    // Record exclusions.
+    
+    exclusions.resize(force.getNumParticles());
+    for (int i = 0; i < (int) force.getNumExclusions(); i++) {
+        int p1, p2;
+        force.getExclusionParticles(i, p1, p2);
+        exclusions[p1].insert(p2);
+        exclusions[p2].insert(p1);
+    }
 }
 
 ReferenceCustomManyParticleIxn::~ReferenceCustomManyParticleIxn( ){
@@ -84,7 +116,7 @@ void ReferenceCustomManyParticleIxn::setUseCutoff(RealOpenMM distance) {
 }
 
 void ReferenceCustomManyParticleIxn::setPeriodic(RealVec& boxSize) {
-    assert(cutoff);
+    assert(useCutoff);
     assert(boxSize[0] >= 2.0*cutoffDistance);
     assert(boxSize[1] >= 2.0*cutoffDistance);
     assert(boxSize[2] >= 2.0*cutoffDistance);
@@ -112,6 +144,23 @@ void ReferenceCustomManyParticleIxn::loopOverInteractions(vector<int>& particles
 
 void ReferenceCustomManyParticleIxn::calculateOneIxn(const vector<int>& particles, vector<RealVec>& atomCoordinates,
                         map<string, double>& variables, vector<RealVec>& forces, RealOpenMM* totalEnergy) const {
+    // Decide whether to include this interaction.
+    
+    for (int i = 0; i < (int) particles.size(); i++) {
+        int p1 = particles[i];
+        for (int j = i+1; j < (int) particles.size(); j++) {
+            int p2 = particles[j];
+            if (exclusions[p1].find(p2) != exclusions[p1].end())
+                return;
+            if (useCutoff) {
+                RealOpenMM delta[ReferenceForce::LastDeltaRIndex];
+                computeDelta(p1, p2, delta, atomCoordinates);
+                if (delta[ReferenceForce::RIndex] >= cutoffDistance)
+                    return;
+            }
+        }
+    }
+    
     // Compute all of the variables the energy can depend on.
 
     for (int i = 0; i < (int) particleTerms.size(); i++) {
@@ -122,8 +171,6 @@ void ReferenceCustomManyParticleIxn::calculateOneIxn(const vector<int>& particle
         const DistanceTermInfo& term = distanceTerms[i];
         computeDelta(particles[term.p1], particles[term.p2], term.delta, atomCoordinates);
         variables[term.name] = term.delta[ReferenceForce::RIndex];
-        if (useCutoff && term.delta[ReferenceForce::RIndex] > cutoffDistance)
-            return;
     }
     for (int i = 0; i < (int) angleTerms.size(); i++) {
         const AngleTermInfo& term = angleTerms[i];

platforms/reference/tests/TestReferenceCustomManyParticleForce.cpp

@@ -41,6 +41,7 @@
 #include "ReferencePlatform.h"
 #include "openmm/CustomManyParticleForce.h"
 #include "openmm/System.h"
+#include "openmm/TabulatedFunction.h"
 #include "openmm/VerletIntegrator.h"
 #include "sfmt/SFMT.h"
 #include <iostream>
@@ -191,6 +192,31 @@ void testPeriodic() {
     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 testParameters() {
     // Create a system.
     
@@ -217,7 +243,7 @@ void testParameters() {
     
     // See if the energy is correct.
 
-    State state1 = context.getState(State::Energy);
+    State state = context.getState(State::Energy);
     double expectedEnergy = 0;
     for (int i = 0; i < numParticles; i++)
         for (int j = i+1; j < numParticles; j++)
@@ -230,7 +256,93 @@ void testParameters() {
                 double r23 = sqrt(d23.dot(d23));
                 expectedEnergy += 2.0*(i+1)*(j+1)*(k+1)*(r12+r13+r23);
             }
-    ASSERT_EQUAL_TOL(expectedEnergy, state1.getPotentialEnergy(), 1e-5);
+    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);
+    ReferencePlatform 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);
 }
 
 int main() {
@@ -238,7 +350,9 @@ int main() {
         testNoCutoff();
         testCutoff();
         testPeriodic();
+        testExclusions();
         testParameters();
+        testTabulatedFunctions();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;