Created API and reference implementation of restricting CustomNonbondedForce...

Created API and reference implementation of restricting CustomNonbondedForce to selected interaction groups

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -80,6 +80,25 @@ namespace OpenMM {
  * if the labels 1 and 2 are reversed.  In contrast, if it depended on the difference sigma1-sigma2, the results would
  * be undefined, because reversing the labels 1 and 2 would change the energy.
  * 
+ * CustomNonbondedForce can operate in two modes.  By default, it computes the interaction of every particle in the System
+ * with every other particle.  Alternatively, you can restrict it to only a subset of particle pairs.  To do this, specify
+ * one or more "interaction groups".  An interaction group consists of two sets of particles that should interact with
+ * each other.  Every particle in the first set interacts with every particle in the second set.  For example, you might use
+ * this feature to compute a solute-solvent interaction energy, while omitting all interactions between two solute atoms
+ * or two solvent atoms.
+ * 
+ * To create an interaction group, call addInteractionGroup().  You may add as many interaction groups as you want.
+ * Be aware of the following:
+ * 
+ * <ul>
+ * <li>Exclusions are still taken into account, so the interactions between excluded pairs are omitted.</li>
+ * <li>Likewise, a particle will never interact with itself, even if it appears in both sets of an interaction group.</li>
+ * <li>If a particle pair appears in two different interaction groups, its interaction will be computed twice.  This is
+ * sometimes useful, but be aware of it so you do not accidentally create unwanted duplicate interactions.</li>
+ * <li>If you do not add any interaction groups to a CustomNonbondedForce, it operates in the default mode where every
+ * particle interacts with every other particle.</li>
+ * </ul>
+ * 
  * When using a cutoff, by default the interaction is sharply truncated at the cutoff distance.
  * Optionally you can instead use a switching function to make the interaction smoothly go to zero over a finite
  * distance range.  To enable this, call setUseSwitchingFunction().  You must also call setSwitchingDistance()
@@ -167,6 +186,12 @@ public:
     int getNumFunctions() const {
         return functions.size();
     }
+    /**
+     * Get the number of interaction groups that have been defined.
+     */
+    int getNumInteractionGroups() const {
+        return interactionGroups.size();
+    }
     /**
      * Get the algebraic expression that gives the interaction energy between two particles
      */
@@ -363,10 +388,26 @@ public:
      * @param max            the value of the independent variable corresponding to the last element of values
      */
     void setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max);
+    /**
+     * Add an interaction group.  An interaction will be computed between every particle in set1 and every particle in set2.
+     * 
+     * @param set1    the first set of particles forming the interaction group
+     * @param set2    the second set of particles forming the interaction group
+     * @return the index of the interaction group that was added
+     */
+    int addInteractionGroup(const std::set<int>& set1, const std::set<int>& set2);
+    /**
+     * Get the parameters for an interaction group.
+     * 
+     * @param index   the index of the interaction group for which to get parameters
+     * @param set1    the first set of particles forming the interaction group
+     * @param set2    the second set of particles forming the interaction group
+     */
+    void getInteractionGroupParameters(int index, std::set<int>& set1, std::set<int>& set2) const;
     /**
      * Update the per-particle parameters in a Context to match those stored in this Force object.  This method provides
      * an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
-     * Simply call setParticleParameters() to modify this object's parameters, then call updateParametersInState()
+     * Simply call setParticleParameters() to modify this object's parameters, then call updateParametersInContext()
      * to copy them over to the Context.
      * 
      * This method has several limitations.  The only information it updates is the values of per-particle parameters.
@@ -383,6 +424,7 @@ private:
     class GlobalParameterInfo;
     class ExclusionInfo;
     class FunctionInfo;
+    class InteractionGroupInfo;
     NonbondedMethod nonbondedMethod;
     double cutoffDistance, switchingDistance;
     bool useSwitchingFunction, useLongRangeCorrection;
@@ -392,6 +434,7 @@ private:
     std::vector<ParticleInfo> particles;
     std::vector<ExclusionInfo> exclusions;
     std::vector<FunctionInfo> functions;
+    std::vector<InteractionGroupInfo> interactionGroups;
 };
 
 /**
@@ -465,6 +508,20 @@ public:
     }
 };
 
+/**
+ * This is an internal class used to record information about an interaction group.
+ * @private
+ */
+class CustomNonbondedForce::InteractionGroupInfo {
+public:
+    std::set<int> set1, set2;
+    InteractionGroupInfo() {
+    }
+    InteractionGroupInfo(const std::set<int>& set1, const std::set<int>& set2) :
+        set1(set1), set2(set2) {
+    }
+};
+
 } // namespace OpenMM
 
 #endif /*OPENMM_CUSTOMNONBONDEDFORCE_H_*/

openmmapi/src/CustomNonbondedForce.cpp

@@ -199,6 +199,17 @@ void CustomNonbondedForce::setFunctionParameters(int index, const std::string& n
     functions[index].max = max;
 }
 
+int CustomNonbondedForce::addInteractionGroup(const std::set<int>& set1, const std::set<int>& set2) {
+    interactionGroups.push_back(InteractionGroupInfo(set1, set2));
+    return interactionGroups.size()-1;
+}
+
+void CustomNonbondedForce::getInteractionGroupParameters(int index, std::set<int>& set1, std::set<int>& set2) const {
+    ASSERT_VALID_INDEX(index, interactionGroups);
+    set1 = interactionGroups[index].set1;
+    set2 = interactionGroups[index].set2;
+}
+
 ForceImpl* CustomNonbondedForce::createImpl() const {
     return new CustomNonbondedForceImpl(*this);
 }

platforms/reference/include/ReferenceCustomNonbondedIxn.h

@@ -29,6 +29,8 @@
 #include "ReferenceNeighborList.h"
 #include "lepton/ExpressionProgram.h"
 #include <map>
+#include <set>
+#include <utility>
 #include <vector>
 
 // ---------------------------------------------------------------------------------------
@@ -47,6 +49,7 @@ class ReferenceCustomNonbondedIxn {
       Lepton::ExpressionProgram forceExpression;
       std::vector<std::string> paramNames;
       std::vector<std::string> particleParamNames;
+      std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups;
 
       /**---------------------------------------------------------------------------------------
 
@@ -97,6 +100,17 @@ class ReferenceCustomNonbondedIxn {
 
       void setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors );
 
+      /**---------------------------------------------------------------------------------------
+
+         Restrict the force to a list of interaction groups.
+
+         @param distance            the cutoff distance
+         @param neighbors           the neighbor list to use
+
+         --------------------------------------------------------------------------------------- */
+
+      void setInteractionGroups(const std::vector<std::pair<std::set<int>, std::set<int> > >& groups);
+
       /**---------------------------------------------------------------------------------------
       
          Set the force to use a switching function.
@@ -126,10 +140,8 @@ class ReferenceCustomNonbondedIxn {
          @param numberOfAtoms    number of atoms
          @param atomCoordinates  atom coordinates
          @param atomParameters   atom parameters (charges, c6, c12, ...)     atomParameters[atomIndex][paramterIndex]
-         @param exclusions       atom exclusion indices                      exclusions[atomIndex][atomToExcludeIndex]
-                                 exclusions[atomIndex][0] = number of exclusions
-                                 exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
-                                 interacting w/ atom atomIndex
+         @param exclusions       atom exclusion indices
+                                 exclusions[atomIndex] contains the list of exclusions for that atom
          @param fixedParameters  non atom parameters (not currently used)
          @param globalParameters the values of global parameters
          @param forces           force array (forces added)
@@ -139,7 +151,7 @@ class ReferenceCustomNonbondedIxn {
          --------------------------------------------------------------------------------------- */
 
       void calculatePairIxn( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
-                            RealOpenMM** atomParameters, int** exclusions,
+                            RealOpenMM** atomParameters, std::vector<std::set<int> >& exclusions,
                             RealOpenMM* fixedParameters, const std::map<std::string, double>& globalParameters,
                             std::vector<OpenMM::RealVec>& forces, RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
 

platforms/reference/include/ReferenceKernels.h

@@ -615,7 +615,6 @@ public:
     void copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force);
 private:
     int numParticles;
-    int **exclusionArray;
     RealOpenMM **particleParamArray;
     RealOpenMM nonbondedCutoff, switchingDistance, periodicBoxSize[3], longRangeCoefficient;
     bool useSwitchingFunction, hasInitializedLongRangeCorrection;
@@ -624,6 +623,7 @@ private:
     std::vector<std::set<int> > exclusions;
     Lepton::ExpressionProgram energyExpression, forceExpression;
     std::vector<std::string> parameterNames, globalParameterNames;
+    std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups;
     NonbondedMethod nonbondedMethod;
     NeighborList* neighborList;
 };

platforms/reference/src/ReferenceKernels.cpp

@@ -1002,7 +1002,6 @@ public:
 
 ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {
     disposeRealArray(particleParamArray, numParticles);
-    disposeIntArray(exclusionArray, numParticles);
     if (neighborList != NULL)
         delete neighborList;
     if (forceCopy != NULL)
@@ -1032,14 +1031,6 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
         for (int j = 0; j < numParameters; j++)
             particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
     }
-    exclusionArray = new int*[numParticles];
-    for (int i = 0; i < numParticles; ++i) {
-        exclusionArray[i] = new int[exclusions[i].size()+1];
-        exclusionArray[i][0] = exclusions[i].size();
-        int index = 0;
-        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
-            exclusionArray[i][++index] = *iter;
-    }
     nonbondedMethod = CalcCustomNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
     nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
     if (nonbondedMethod == NoCutoff) {
@@ -1090,6 +1081,14 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
         longRangeCoefficient = 0.0;
         hasInitializedLongRangeCorrection = true;
     }
+    
+    // Record the interaction groups.
+    
+    for (int i = 0; i < force.getNumInteractionGroups(); i++) {
+        set<int> set1, set2;
+        force.getInteractionGroupParameters(i, set1, set2);
+        interactionGroups.push_back(make_pair(set1, set2));
+    }
 }
 
 double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -1109,6 +1108,8 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
             throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
         ixn.setPeriodic(box);
     }
+    if (interactionGroups.size() > 0)
+        ixn.setInteractionGroups(interactionGroups);
     bool globalParamsChanged = false;
     for (int i = 0; i < (int) globalParameterNames.size(); i++) {
         double value = context.getParameter(globalParameterNames[i]);
@@ -1118,7 +1119,7 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
     }
     if (useSwitchingFunction)
         ixn.setUseSwitchingFunction(switchingDistance);
-    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParamValues, forceData, 0, includeEnergy ? &energy : NULL);
+    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusions, 0, globalParamValues, forceData, 0, includeEnergy ? &energy : NULL);
     
     // Add in the long range correction.
     

platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp

@@ -32,8 +32,10 @@
 #include "ReferenceCustomNonbondedIxn.h"
 
 using std::map;
+using std::pair;
 using std::string;
 using std::stringstream;
+using std::set;
 using std::vector;
 using OpenMM::RealVec;
 
@@ -94,6 +96,19 @@ ReferenceCustomNonbondedIxn::~ReferenceCustomNonbondedIxn( ){
     neighborList = &neighbors;
   }
 
+/**---------------------------------------------------------------------------------------
+
+   Restrict the force to a list of interaction groups.
+
+   @param distance            the cutoff distance
+   @param neighbors           the neighbor list to use
+
+   --------------------------------------------------------------------------------------- */
+
+void ReferenceCustomNonbondedIxn::setInteractionGroups(const vector<pair<set<int>, set<int> > >& groups) {
+    interactionGroups = groups;
+}
+
 /**---------------------------------------------------------------------------------------
 
    Set the force to use a switching function.
@@ -138,10 +153,8 @@ void ReferenceCustomNonbondedIxn::setUseSwitchingFunction( RealOpenMM distance )
    @param numberOfAtoms    number of atoms
    @param atomCoordinates  atom coordinates
    @param atomParameters   atom parameters                             atomParameters[atomIndex][paramterIndex]
-   @param exclusions       atom exclusion indices                      exclusions[atomIndex][atomToExcludeIndex]
-                           exclusions[atomIndex][0] = number of exclusions
-                           exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
-                           interacting w/ atom atomIndex
+   @param exclusions       atom exclusion indices
+                           exclusions[atomIndex] contains the list of exclusions for that atom
    @param fixedParameters  non atom parameters (not currently used)
    @param globalParameters the values of global parameters
    @param forces           force array (forces added)
@@ -151,53 +164,57 @@ void ReferenceCustomNonbondedIxn::setUseSwitchingFunction( RealOpenMM distance )
    --------------------------------------------------------------------------------------- */
 
 void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates,
-                                             RealOpenMM** atomParameters, int** exclusions,
+                                             RealOpenMM** atomParameters, vector<set<int> >& exclusions,
                                              RealOpenMM* fixedParameters, const map<string, double>& globalParameters, vector<RealVec>& forces,
                                              RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
 
-   map<string, double> variables = globalParameters;
-   if (cutoff) {
-       for (int i = 0; i < (int) neighborList->size(); i++) {
-           OpenMM::AtomPair pair = (*neighborList)[i];
-           for (int j = 0; j < (int) paramNames.size(); j++) {
-               variables[particleParamNames[j*2]] = atomParameters[pair.first][j];
-               variables[particleParamNames[j*2+1]] = atomParameters[pair.second][j];
-           }
-           calculateOneIxn(pair.first, pair.second, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
-       }
-   }
-   else {
-       // allocate and initialize exclusion array
-
-       int* exclusionIndices = new int[numberOfAtoms];
-       for( int ii = 0; ii < numberOfAtoms; ii++ ){
-          exclusionIndices[ii] = -1;
-       }
-
-       for( int ii = 0; ii < numberOfAtoms; ii++ ){
-
-          // set exclusions
-
-          for( int jj = 1; jj <= exclusions[ii][0]; jj++ ){
-             exclusionIndices[exclusions[ii][jj]] = ii;
-          }
-
-          // loop over atom pairs
-
-          for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
-
-             if( exclusionIndices[jj] != ii ){
-                 for (int j = 0; j < (int) paramNames.size(); j++) {
-                     variables[particleParamNames[j*2]] = atomParameters[ii][j];
-                     variables[particleParamNames[j*2+1]] = atomParameters[jj][j];
-                 }
-                 calculateOneIxn(ii, jj, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
-             }
-          }
-       }
-
-       delete[] exclusionIndices;
-   }
+    map<string, double> variables = globalParameters;
+    if (interactionGroups.size() > 0) {
+        // The user has specified interaction groups, so compute only the requested interactions.
+        
+        for (int group = 0; group < (int) interactionGroups.size(); group++) {
+            const set<int>& set1 = interactionGroups[group].first;
+            const set<int>& set2 = interactionGroups[group].second;
+            for (set<int>::const_iterator atom1 = set1.begin(); atom1 != set1.end(); ++atom1) {
+                for (set<int>::const_iterator atom2 = set2.begin(); atom2 != set2.end(); ++atom2) {
+                    if (*atom1 != *atom2 && exclusions[*atom1].find(*atom2) == exclusions[*atom1].end()) {
+                        for (int j = 0; j < (int) paramNames.size(); j++) {
+                            variables[particleParamNames[j*2]] = atomParameters[*atom1][j];
+                            variables[particleParamNames[j*2+1]] = atomParameters[*atom2][j];
+                        }
+                        calculateOneIxn(*atom1, *atom2, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
+                    }
+                }
+            }
+        }
+    }
+    else if (cutoff) {
+        // We are using a cutoff, so get the interactions from the neighbor list.
+        
+        for (int i = 0; i < (int) neighborList->size(); i++) {
+            OpenMM::AtomPair pair = (*neighborList)[i];
+            for (int j = 0; j < (int) paramNames.size(); j++) {
+                variables[particleParamNames[j*2]] = atomParameters[pair.first][j];
+                variables[particleParamNames[j*2+1]] = atomParameters[pair.second][j];
+            }
+            calculateOneIxn(pair.first, pair.second, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
+        }
+    }
+    else {
+        // Every particle interacts with every other one.
+        
+        for (int ii = 0; ii < numberOfAtoms; ii++) {
+            for (int jj = ii+1; jj < numberOfAtoms; jj++) {
+                if (exclusions[jj].find(ii) == exclusions[jj].end()) {
+                    for (int j = 0; j < (int) paramNames.size(); j++) {
+                        variables[particleParamNames[j*2]] = atomParameters[ii][j];
+                        variables[particleParamNames[j*2+1]] = atomParameters[jj][j];
+                    }
+                    calculateOneIxn(ii, jj, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
+                }
+            }
+        }
+    }
 }
 
   /**---------------------------------------------------------------------------------------

platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp

@@ -43,6 +43,7 @@
 #include "openmm/System.h"
 #include "openmm/VerletIntegrator.h"
 #include <iostream>
+#include <set>
 #include <vector>
 
 using namespace OpenMM;
@@ -470,6 +471,43 @@ void testLongRangeCorrection() {
     ASSERT_EQUAL_TOL(standardEnergy1-standardEnergy2, customEnergy1-customEnergy2, 1e-4);
 }
 
+void testInteractionGroups() {
+    const int numParticles = 6;
+    ReferencePlatform platform;
+    System system;
+    VerletIntegrator integrator(0.01);
+    CustomNonbondedForce* nonbonded = new CustomNonbondedForce("v1+v2");
+    nonbonded->addPerParticleParameter("v");
+    vector<double> params(1, 0.001);
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        nonbonded->addParticle(params);
+        params[0] *= 10;
+    }
+    set<int> set1, set2, set3, set4;
+    set1.insert(2);
+    set2.insert(0);
+    set2.insert(1);
+    set2.insert(2);
+    set2.insert(3);
+    set2.insert(4);
+    set2.insert(5);
+    nonbonded->addInteractionGroup(set1, set2); // Particle 2 interacts with every other particle.
+    set3.insert(0);
+    set3.insert(1);
+    set4.insert(4);
+    set4.insert(5);
+    nonbonded->addInteractionGroup(set3, set4); // Particles 0 and 1 interact with 4 and 5.
+    nonbonded->addExclusion(1, 2); // Add an exclusion to make sure it gets skipped.
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(numParticles);
+    context.setPositions(positions);
+    State state = context.getState(State::Energy);
+    double expectedEnergy = 331.423; // Each digit is the number of interactions a particle particle is involved in.
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), TOL);
+}
+
 int main() {
     try {
         testSimpleExpression();
@@ -481,6 +519,7 @@ int main() {
         testCoulombLennardJones();
         testSwitchingFunction();
         testLongRangeCorrection();
+        testInteractionGroups();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;