Created CustomNonbondedForce::createExclusionsFromBonds()

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -343,6 +343,8 @@ public:
     void setParticleParameters(int index, const std::vector<double>& parameters);
     /**
      * Add a particle pair to the list of interactions that should be excluded.
+     * 
+     * In many cases, you can use createExclusionsFromBonds() rather than adding each exclusion explicitly.
      *
      * @param particle1  the index of the first particle in the pair
      * @param particle2  the index of the second particle in the pair
@@ -365,6 +367,15 @@ public:
      * @param particle2  the index of the second particle in the pair
      */
     void setExclusionParticles(int index, int particle1, int particle2);
+    /**
+     * Identify exclusions based on the molecular topology.  Particles which are separated by up to a specified number of
+     * bonds are added as exclusions.
+     *
+     * @param bonds       the set of bonds based on which to construct exclusions.  Each element specifies the indices of
+     *                    two particles that are bonded to each other.
+     * @param bondCutoff  pairs of particles that are separated by this many bonds or fewer are added to the list of exclusions
+     */
+    void createExclusionsFromBonds(const std::vector<std::pair<int, int> >& bonds, int bondCutoff);
     /**
      * Add a tabulated function that may appear in the energy expression.
      *

openmmapi/src/CustomNonbondedForce.cpp

@@ -174,6 +174,31 @@ void CustomNonbondedForce::setExclusionParticles(int index, int particle1, int p
     exclusions[index].particle1 = particle1;
     exclusions[index].particle2 = particle2;
 }
+
+void CustomNonbondedForce::createExclusionsFromBonds(const vector<pair<int, int> >& bonds, int bondCutoff) {
+    vector<set<int> > exclusions(particles.size());
+    vector<set<int> > bonded12(exclusions.size());
+    for (int i = 0; i < (int) bonds.size(); ++i) {
+        int p1 = bonds[i].first;
+        int p2 = bonds[i].second;
+        exclusions[p1].insert(p2);
+        exclusions[p2].insert(p1);
+        bonded12[p1].insert(p2);
+        bonded12[p2].insert(p1);
+    }
+    for (int level = 0; level < bondCutoff-1; level++) {
+        vector<set<int> > currentExclusions = exclusions;
+        for (int i = 0; i < (int) particles.size(); i++) {
+            for (set<int>::const_iterator iter = currentExclusions[i].begin(); iter != currentExclusions[i].end(); ++iter)
+                exclusions[*iter].insert(bonded12[i].begin(), bonded12[i].end());
+        }
+    }
+    for (int i = 0; i < (int) exclusions.size(); ++i)
+        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
+            if (*iter < i)
+                addExclusion(*iter, i);
+}
+
 int CustomNonbondedForce::addTabulatedFunction(const std::string& name, TabulatedFunction* function) {
     functions.push_back(FunctionInfo(name, function));
     return functions.size()-1;

tests/TestFindExclusions.cpp

@@ -30,6 +30,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "openmm/internal/AssertionUtilities.h"
+#include "openmm/CustomNonbondedForce.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/OpenMMException.h"
 #include <iostream>
@@ -166,10 +167,62 @@ void testReplaceExceptions() {
     ASSERT(charge == 5.0);
 }
 
+/**
+ * This is the same as testFindExceptions(), except it tests adding exclusions to a CustomNonbondedForce.
+ */
+
+void testFindCustomExclusions() {
+    CustomNonbondedForce nonbonded("r");
+    vector<pair<int, int> > bonds;
+    vector<double> params;
+    for (int i = 0; i < NUM_ATOMS; i++)
+        nonbonded.addParticle(params);
+    // loop over all main-chain atoms (even numbered atoms)
+    for (int i = 0; i < NUM_ATOMS-1; i += 2)
+    {
+        // side-chain bonds
+        bonds.push_back(pair<int, int>(i, i+1));
+        // main-chain bonds
+        if (i < NUM_ATOMS-2) // penultimate atom (NUM_ATOMS-2) has no subsequent main-chain atom
+            bonds.push_back(pair<int, int>(i, i+2));
+    }
+    nonbonded.createExclusionsFromBonds(bonds, 3);
+
+    // Build lists of the expected exclusions.
+
+    vector<set<int> > expectedExclusions(NUM_ATOMS);
+    int totalExclusions = 0;
+    for (int i = 0; i < NUM_ATOMS; i += 2) {
+        addAtomsToExclusions(i, i+1, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i, i+2, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i, i+3, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i, i+4, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i+1, i+2, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i, i+5, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i, i+6, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i+1, i+3, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i+1, i+4, expectedExclusions, totalExclusions);
+    }
+    for (int i = 0; i < nonbonded.getNumExclusions(); i++) {
+        int particle1, particle2;
+        nonbonded.getExclusionParticles(i, particle1, particle2);
+    }
+
+    // Compare them to the exceptions that were generated.
+
+    ASSERT_EQUAL(totalExclusions, nonbonded.getNumExclusions());
+    for (int i = 0; i < nonbonded.getNumExclusions(); i++) {
+        int particle1, particle2;
+        nonbonded.getExclusionParticles(i, particle1, particle2);
+        ASSERT(expectedExclusions[particle1].find(particle2) != expectedExclusions[particle1].end());
+    }
+}
+
 int main() {
     try {
         testFindExceptions();
         testReplaceExceptions();
+        testFindCustomExclusions();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;