CustomNonbondedForce supports named parameters, global parameters, more flexible combining rules

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -89,7 +89,13 @@ public:
      * Get the number of per-particle parameters that the interaction depends on.
      */
     int getNumParameters() const {
-        return combiningRules.size();
+        return parameters.size();
+    }
+    /**
+     * Get the number of global parameters that the interaction depends on.
+     */
+    int getNumGlobalParameters() const {
+        return globalParameters.size();
     }
     /**
      * Get the algebraic expression that gives the interaction energy between two particles
@@ -144,10 +150,25 @@ public:
     /**
      * Add a new per-particle parmeter that the interaction may depend on.
      *
+     * @param name             the name of the parameter
      * @param combiningRule    an algebraic expression giving the combining rule for this parameter
      * @return the index of the parameter that was added
      */
-    int addParameter(const std::string& combiningRule);
+    int addParameter(const std::string& name, const std::string& combiningRule);
+    /**
+     * Get the name of a per-particle parameter.
+     *
+     * @param index     the index of the parameter for which to get the name
+     * @return the parameter name
+     */
+    const std::string& getParameterName(int index) const;
+    /**
+     * Set the name of a per-particle parameter.
+     *
+     * @param index          the index of the parameter for which to set the name
+     * @param name           the name of the parameter
+     */
+    void setParameterName(int index, const std::string& name);
     /**
      * Get the combining rule for a per-particle parameter.
      *
@@ -162,6 +183,27 @@ public:
      * @param combiningRule  an algebraic expression giving the combining rule for the parameter
      */
     void setParameterCombiningRule(int index, const std::string& combiningRule);
+    /**
+     * Add a new global parmeter that the interaction may depend on.
+     *
+     * @param name             the name of the parameter
+     * @return the index of the parameter that was added
+     */
+    int addGlobalParameter(const std::string& name);
+    /**
+     * Get the name of a global parameter.
+     *
+     * @param index     the index of the parameter for which to get the name
+     * @return the parameter name
+     */
+    const std::string& getGlobalParameterName(int index) const;
+    /**
+     * Set the name of a global parameter.
+     *
+     * @param index          the index of the parameter for which to set the name
+     * @param name           the name of the parameter
+     */
+    void setGlobalParameterName(int index, const std::string& name);
     /**
      * Add the nonbonded force parameters for a particle.  This should be called once for each particle
      * in the System.  When it is called for the i'th time, it specifies the parameters for the i'th particle.
@@ -220,12 +262,14 @@ protected:
     ForceImpl* createImpl();
 private:
     class ParticleInfo;
+    class ParameterInfo;
     class ExceptionInfo;
     NonbondedMethod nonbondedMethod;
     double cutoffDistance;
     Vec3 periodicBoxVectors[3];
     std::string energyExpression;
-    std::vector<std::string> combiningRules;
+    std::vector<ParameterInfo> parameters;
+    std::vector<std::string> globalParameters;
     std::vector<ParticleInfo> particles;
     std::vector<ExceptionInfo> exceptions;
     std::map<std::pair<int, int>, int> exceptionMap;
@@ -240,6 +284,15 @@ public:
     }
 };
 
+class CustomNonbondedForce::ParameterInfo {
+public:
+    std::string name, combiningRule;
+    ParameterInfo() {
+    }
+    ParameterInfo(const std::string& name, const std::string& combiningRule) : name(name), combiningRule(combiningRule) {
+    }
+};
+
 class CustomNonbondedForce::ExceptionInfo {
 public:
     int particle1, particle2;

openmmapi/include/openmm/internal/CustomNonbondedForceImpl.h

@@ -36,6 +36,7 @@
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/Kernel.h"
 #include <utility>
+#include <map>
 #include <string>
 
 namespace OpenMM {
@@ -57,9 +58,7 @@ public:
     }
     void calcForces(ContextImpl& context, Stream& forces);
     double calcEnergy(ContextImpl& context);
-    std::map<std::string, double> getDefaultParameters() {
-        return std::map<std::string, double>(); // This force field doesn't define any parameters.
-    }
+    std::map<std::string, double> getDefaultParameters();
     std::vector<std::string> getKernelNames();
 private:
     CustomNonbondedForce& owner;

openmmapi/src/CustomNonbondedForce.cpp

@@ -94,17 +94,38 @@ void CustomNonbondedForce::setPeriodicBoxVectors(Vec3 a, Vec3 b, Vec3 c) {
     periodicBoxVectors[2] = c;
 }
 
-int CustomNonbondedForce::addParameter(const string& combiningRule) {
-    combiningRules.push_back(combiningRule);
-    return combiningRules.size()-1;
+int CustomNonbondedForce::addParameter(const string& name, const string& combiningRule) {
+    parameters.push_back(ParameterInfo(name, combiningRule));
+    return parameters.size()-1;
+}
+
+const string& CustomNonbondedForce::getParameterName(int index) const {
+    return parameters[index].name;
+}
+
+void CustomNonbondedForce::setParameterName(int index, const string& name) {
+    parameters[index].name = name;
 }
 
 const string& CustomNonbondedForce::getParameterCombiningRule(int index) const {
-    return combiningRules[index];
+    return parameters[index].combiningRule;
 }
 
 void CustomNonbondedForce::setParameterCombiningRule(int index, const string& combiningRule) {
-    combiningRules[index] = combiningRule;
+    parameters[index].combiningRule = combiningRule;
+}
+
+int CustomNonbondedForce::addGlobalParameter(const string& name) {
+    globalParameters.push_back(name);
+    return globalParameters.size()-1;
+}
+
+const string& CustomNonbondedForce::getGlobalParameterName(int index) const {
+    return globalParameters[index];
+}
+
+void CustomNonbondedForce::setGlobalParameterName(int index, const string& name) {
+    globalParameters[index] = name;
 }
 
 int CustomNonbondedForce::addParticle(const vector<double>& parameters) {

openmmapi/src/CustomNonbondedForceImpl.cpp

@@ -36,9 +36,11 @@
 #include <sstream>
 
 using namespace OpenMM;
+using std::map;
 using std::pair;
 using std::vector;
 using std::set;
+using std::string;
 using std::stringstream;
 
 CustomNonbondedForceImpl::CustomNonbondedForceImpl(CustomNonbondedForce& owner) : owner(owner) {
@@ -110,9 +112,15 @@ double CustomNonbondedForceImpl::calcEnergy(ContextImpl& context) {
     return dynamic_cast<CalcCustomNonbondedForceKernel&>(kernel.getImpl()).executeEnergy(context);
 }
 
-std::vector<std::string> CustomNonbondedForceImpl::getKernelNames() {
-    std::vector<std::string> names;
+vector<string> CustomNonbondedForceImpl::getKernelNames() {
+    vector<string> names;
     names.push_back(CalcCustomNonbondedForceKernel::Name());
     return names;
 }
 
+map<string, double> CustomNonbondedForceImpl::getDefaultParameters() {
+    map<string, double> parameters;
+    for (int i = 0; i < owner.getNumGlobalParameters(); i++)
+        parameters[owner.getGlobalParameterName(i)] = 0.0;
+    return parameters;
+}

platforms/reference/src/ReferenceKernels.cpp

@@ -545,14 +545,18 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
     Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
     energyExpression = expression.createProgram();
     forceExpression = expression.differentiate("r").optimize().createProgram();
-    for (int i = 0; i < numParameters; i++)
+    for (int i = 0; i < numParameters; i++) {
+        parameterNames.push_back(force.getParameterName(i));
         combiningRules.push_back(Lepton::Parser::parse(force.getParameterCombiningRule(i)).optimize().createProgram());
+    }
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        globalParameterNames.push_back(force.getGlobalParameterName(i));
 }
 
 void ReferenceCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
     RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
     RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
-    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, combiningRules);
+    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames, combiningRules);
     bool periodic = (nonbondedMethod == CutoffPeriodic);
     if (nonbondedMethod != NoCutoff) {
         computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
@@ -560,15 +564,18 @@ void ReferenceCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context
     }
     if (periodic)
         ixn.setPeriodic(periodicBoxSize);
-    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
-    ixn.calculateExceptionIxn(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0);
+    map<string, double> globalParameters;
+    for (int i = 0; i < globalParameterNames.size(); i++)
+        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, 0);
+    ixn.calculateExceptionIxn(num14, bonded14IndexArray, posData, bonded14ParamArray, globalParameters, forceData, 0, 0);
 }
 
 double ReferenceCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
     RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
     RealOpenMM** forceData = allocateRealArray(numParticles, 3);
     RealOpenMM energy = 0;
-    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, combiningRules);
+    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames, combiningRules);
     bool periodic = (nonbondedMethod == CutoffPeriodic);
     if (nonbondedMethod != NoCutoff) {
         computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
@@ -576,8 +583,11 @@ double ReferenceCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& conte
     }
     if (periodic)
         ixn.setPeriodic(periodicBoxSize);
-    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
-    ixn.calculateExceptionIxn(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, &energy);
+    map<string, double> globalParameters;
+    for (int i = 0; i < globalParameterNames.size(); i++)
+        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, &energy);
+    ixn.calculateExceptionIxn(num14, bonded14IndexArray, posData, bonded14ParamArray, globalParameters, forceData, 0, &energy);
     disposeRealArray(forceData, numParticles);
     return energy;
 }

platforms/reference/src/ReferenceKernels.h

@@ -310,6 +310,7 @@ private:
     RealOpenMM nonbondedCutoff, periodicBoxSize[3];
     std::vector<std::set<int> > exclusions;
     Lepton::ExpressionProgram energyExpression, forceExpression;
+    std::vector<std::string> parameterNames, globalParameterNames;
     std::vector<Lepton::ExpressionProgram> combiningRules;
     NonbondedMethod nonbondedMethod;
     NeighborList* neighborList;

platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp

@@ -43,9 +43,9 @@ using std::vector;
    --------------------------------------------------------------------------------------- */
 
 ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::ExpressionProgram& energyExpression,
-        const Lepton::ExpressionProgram& forceExpression, const std::vector<Lepton::ExpressionProgram>& combiningRules) :
+        const Lepton::ExpressionProgram& forceExpression, const vector<string>& parameterNames, const vector<Lepton::ExpressionProgram>& combiningRules) :
             cutoff(false), periodic(false), energyExpression(energyExpression), forceExpression(forceExpression),
-            combiningRules(combiningRules) {
+            paramNames(parameterNames), combiningRules(combiningRules) {
 
    // ---------------------------------------------------------------------------------------
 
@@ -53,10 +53,12 @@ ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::Expressio
 
    // ---------------------------------------------------------------------------------------
 
-    for (int i = 0; i < combiningRules.size(); i++) {
-        stringstream name;
-        name << "param" << i;
-        paramNames.push_back(name.str());
+    for (int i = 0; i < paramNames.size(); i++) {
+        for (int j = 1; j < 3; j++) {
+            stringstream name;
+            name << paramNames[i] << j;
+            particleParamNames.push_back(name.str());
+        }
     }
 }
 
@@ -135,6 +137,7 @@ ReferenceCustomNonbondedIxn::~ReferenceCustomNonbondedIxn( ){
                            exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
                            interacting w/ atom atomIndex
    @param fixedParameters  non atom parameters (not currently used)
+   @param globalParameters the values of global parameters
    @param forces           force array (forces added)
    @param energyByAtom     atom energy
    @param totalEnergy      total energy
@@ -145,11 +148,11 @@ ReferenceCustomNonbondedIxn::~ReferenceCustomNonbondedIxn( ){
 
 int ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
                                              RealOpenMM** atomParameters, int** exclusions,
-                                             RealOpenMM* fixedParameters, RealOpenMM** forces,
+                                             RealOpenMM* fixedParameters, map<string, double> globalParameters, RealOpenMM** forces,
                                              RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
 
-   map<string, double> variablesForParams;
-   map<string, double> variablesForForce;
+   map<string, double> variablesForParams = globalParameters;
+   map<string, double> variablesForForce = globalParameters;
    if (cutoff) {
        for (int i = 0; i < (int) neighborList->size(); i++) {
            OpenMM::AtomPair pair = (*neighborList)[i];
@@ -157,10 +160,11 @@ int ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM
            // Apply the combining rules to compute the force field parameters.
 
            for (int j = 0; j < combiningRules.size(); j++) {
-               variablesForParams["particle1"] = atomParameters[pair.first][j];
-               variablesForParams["particle2"] = atomParameters[pair.second][j];
-               variablesForForce[paramNames[j]] = combiningRules[j].evaluate(variablesForParams);
+               variablesForParams[particleParamNames[j*2]] = atomParameters[pair.first][j];
+               variablesForParams[particleParamNames[j*2+1]] = atomParameters[pair.second][j];
            }
+           for (int j = 0; j < combiningRules.size(); j++)
+               variablesForForce[paramNames[j]] = combiningRules[j].evaluate(variablesForParams);
            calculateOneIxn(pair.first, pair.second, atomCoordinates, variablesForForce, forces, energyByAtom, totalEnergy);
        }
    }
@@ -189,10 +193,11 @@ int ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM
                  // Apply the combining rules to compute the force field parameters.
 
                  for (int j = 0; j < combiningRules.size(); j++) {
-                     variablesForParams["particle1"] = atomParameters[ii][j];
-                     variablesForParams["particle2"] = atomParameters[jj][j];
-                     variablesForForce[paramNames[j]] = combiningRules[j].evaluate(variablesForParams);
+                     variablesForParams[particleParamNames[j*2]] = atomParameters[ii][j];
+                     variablesForParams[particleParamNames[j*2+1]] = atomParameters[jj][j];
                  }
+                 for (int j = 0; j < combiningRules.size(); j++)
+                     variablesForForce[paramNames[j]] = combiningRules[j].evaluate(variablesForParams);
                  calculateOneIxn(ii, jj, atomCoordinates, variablesForForce, forces, energyByAtom, totalEnergy);
              }
           }
@@ -212,16 +217,17 @@ int ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM
    @param atomCoordinates      atom coordinates: atomCoordinates[atomIndex][3]
    @param parameters           parameters: parameters[exceptionIndex][*]; contents of array
                                depend on ixn
+   @param globalParameters     the values of global parameters
    @param forces               force array (forces added to current values): forces[atomIndex][3]
    @param totalEnergy          totalEnergy: sum over { energies[atomIndex] }
 
    --------------------------------------------------------------------------------------- */
 
 void ReferenceCustomNonbondedIxn::calculateExceptionIxn( int numberOfExceptions, int** atomIndices, RealOpenMM** atomCoordinates,
-                            RealOpenMM** parameters, RealOpenMM** forces,
+                            RealOpenMM** parameters, map<string, double> globalParameters, RealOpenMM** forces,
                             RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
 
-    map<string, double> variables;
+    map<string, double> variables = globalParameters;
     for (int i = 0; i < numberOfExceptions; i++) {
         for (int j = 0; j < combiningRules.size(); j++)
             variables[paramNames[j]] = parameters[i][j];

platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.h

@@ -33,7 +33,7 @@
 
 // ---------------------------------------------------------------------------------------
 
-class ReferenceCustomNonbondedIxn : public ReferencePairIxn {
+class ReferenceCustomNonbondedIxn {
 
    private:
 
@@ -46,6 +46,7 @@ class ReferenceCustomNonbondedIxn : public ReferencePairIxn {
       Lepton::ExpressionProgram forceExpression;
       std::vector<Lepton::ExpressionProgram> combiningRules;
       std::vector<std::string> paramNames;
+      std::vector<std::string> particleParamNames;
 
       /**---------------------------------------------------------------------------------------
 
@@ -75,7 +76,7 @@ class ReferenceCustomNonbondedIxn : public ReferencePairIxn {
          --------------------------------------------------------------------------------------- */
 
        ReferenceCustomNonbondedIxn(const Lepton::ExpressionProgram& energyExpression, const Lepton::ExpressionProgram& forceExpression,
-                                   const std::vector<Lepton::ExpressionProgram>& combiningRules);
+                                   const std::vector<std::string>& parameterNames, const std::vector<Lepton::ExpressionProgram>& combiningRules);
 
       /**---------------------------------------------------------------------------------------
 
@@ -124,6 +125,7 @@ class ReferenceCustomNonbondedIxn : public ReferencePairIxn {
                                  exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
                                  interacting w/ atom atomIndex
          @param fixedParameters  non atom parameters (not currently used)
+         @param globalParameters the values of global parameters
          @param forces           force array (forces added)
          @param energyByAtom     atom energy
          @param totalEnergy      total energy
@@ -134,8 +136,8 @@ class ReferenceCustomNonbondedIxn : public ReferencePairIxn {
 
       int calculatePairIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
                             RealOpenMM** atomParameters, int** exclusions,
-                            RealOpenMM* fixedParameters, RealOpenMM** forces,
-                            RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
+                            RealOpenMM* fixedParameters, std::map<std::string, double> globalParameters,
+                            RealOpenMM** forces, RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
 
       /**---------------------------------------------------------------------------------------
 
@@ -146,6 +148,7 @@ class ReferenceCustomNonbondedIxn : public ReferencePairIxn {
          @param atomCoordinates      atom coordinates: atomCoordinates[atomIndex][3]
          @param parameters           parameters: parameters[excptionIndex][*]; contents of array
                                      depend on ixn
+         @param globalParameters     the values of global parameters
          @param forces               force array (forces added to current values): forces[atomIndex][3]
          @param energyByAtom         atom energy
          @param totalEnergy          total energy
@@ -153,8 +156,8 @@ class ReferenceCustomNonbondedIxn : public ReferencePairIxn {
          --------------------------------------------------------------------------------------- */
 
       void calculateExceptionIxn( int numberOfExceptions, int** atomIndices, RealOpenMM** atomCoordinates,
-                                    RealOpenMM** parameters, RealOpenMM** forces,
-                                    RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
+                                    RealOpenMM** parameters, std::map<std::string, double> globalParameters,
+                                    RealOpenMM** forces, RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
 
 // ---------------------------------------------------------------------------------------
 

platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp

@@ -77,9 +77,11 @@ void testParameters() {
     system.addParticle(1.0);
     system.addParticle(1.0);
     VerletIntegrator integrator(0.01);
-    CustomNonbondedForce* forceField = new CustomNonbondedForce("param0*(r*param1)^3");
-    forceField->addParameter("particle1*particle2");
-    forceField->addParameter("particle1+particle2");
+    CustomNonbondedForce* forceField = new CustomNonbondedForce("scale*a*(r*b)^3");
+    forceField->addParameter("a", "a1*a2");
+    forceField->addParameter("b", "c+b1+b2");
+    forceField->addGlobalParameter("scale");
+    forceField->addGlobalParameter("c");
     vector<double> params(2);
     params[0] = 1.5;
     params[1] = 2.0;
@@ -93,20 +95,30 @@ void testParameters() {
     positions[0] = Vec3(0, 0, 0);
     positions[1] = Vec3(2, 0, 0);
     context.setPositions(positions);
+    context.setParameter("scale", 1.0);
+    context.setParameter("c", 0.0);
     State state = context.getState(State::Forces | State::Energy);
-    const vector<Vec3>& forces = state.getForces();
+    vector<Vec3> forces = state.getForces();
     double force = -3.0*3*5.0*(10*10);
     ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
     ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
     ASSERT_EQUAL_TOL(3.0*(10*10*10), state.getPotentialEnergy(), TOL);
+    context.setParameter("scale", 1.5);
+    context.setParameter("c", 1.0);
+    state = context.getState(State::Forces | State::Energy);
+    forces = state.getForces();
+    force = -1.5*3.0*3*6.0*(12*12);
+    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_TOL(1.5*3.0*(12*12*12), state.getPotentialEnergy(), TOL);
 }
 
 void testExceptions() {
     ReferencePlatform platform;
     System system;
     VerletIntegrator integrator(0.01);
-    CustomNonbondedForce* nonbonded = new CustomNonbondedForce("param0*r");
-    nonbonded->addParameter("particle1+particle2");
+    CustomNonbondedForce* nonbonded = new CustomNonbondedForce("a*r");
+    nonbonded->addParameter("a", "a1+a2");
     vector<double> params(1);
     vector<Vec3> positions(4);
     for (int i = 0; i < 4; i++) {