Added a "replace" flag to NonbondedForce::addException()

openmmapi/include/openmm/NonbondedForce.h

@@ -36,6 +36,7 @@
 #include "Vec3.h"
 #include <map>
 #include <set>
+#include <utility>
 #include <vector>
 #include "internal/windowsExport.h"
 
@@ -215,9 +216,11 @@ public:
      * @param chargeProd the scaled product of the atomic charges (i.e. the strength of the Coulomb interaction), measured in units of the proton charge squared
      * @param sigma      the sigma parameter of the Lennard-Jones potential (corresponding to the van der Waals radius of the particle), measured in nm
      * @param epsilon    the epsilon parameter of the Lennard-Jones potential (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
+     * @param replace    determines the behavior if there is already an exception for the same two particles.  If true, the existing one is replaced.  If false,
+     *                   an exception is thrown.
      * @return the index of the exception that was added
      */
-    int addException(int particle1, int particle2, double chargeProd, double sigma, double epsilon);
+    int addException(int particle1, int particle2, double chargeProd, double sigma, double epsilon, bool replace = false);
     /**
      * Get the force field parameters for an interaction that should be calculated differently from others.
      * 
@@ -275,6 +278,7 @@ private:
 #endif
     std::vector<ParticleInfo> particles;
     std::vector<ExceptionInfo> exceptions;
+    std::map<std::pair<int, int>, int> exceptionMap;
 #if defined(_MSC_VER)
 #pragma warning(pop)
 #endif

openmmapi/src/NonbondedForce.cpp

@@ -33,12 +33,18 @@
 #include "openmm/OpenMMException.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/internal/NonbondedForceImpl.h"
+#include "openmm/NonBondedForce.h"
 #include <cmath>
+#include <map>
+#include <sstream>
 #include <utility>
 
 using namespace OpenMM;
+using std::map;
 using std::pair;
 using std::set;
+using std::string;
+using std::stringstream;
 using std::vector;
 
 NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), ewaldErrorTol(1e-4) {
@@ -107,9 +113,30 @@ void NonbondedForce::setParticleParameters(int index, double charge, double sigm
     particles[index].epsilon = epsilon;
 }
 
-int NonbondedForce::addException(int particle1, int particle2, double chargeProd, double sigma, double epsilon) {
-    exceptions.push_back(ExceptionInfo(particle1, particle2, chargeProd, sigma, epsilon));
-    return exceptions.size()-1;
+int NonbondedForce::addException(int particle1, int particle2, double chargeProd, double sigma, double epsilon, bool replace) {
+    map<pair<int, int>, int>::iterator iter = exceptionMap.find(pair<int, int>(particle1, particle2));
+    int newIndex;
+    if (iter == exceptionMap.end())
+        iter = exceptionMap.find(pair<int, int>(particle2, particle1));
+    if (iter != exceptionMap.end()) {
+        if (!replace) {
+            stringstream msg;
+            msg << "NonbondedForce: There is already an exception for particles ";
+            msg << particle1;
+            msg << " and ";
+            msg << particle2;
+            throw OpenMMException(msg.str());
+        }
+        exceptions[iter->second] = ExceptionInfo(particle1, particle2, chargeProd, sigma, epsilon);
+        newIndex = iter->second;
+        exceptionMap.erase(iter->first);
+    }
+    else {
+        exceptions.push_back(ExceptionInfo(particle1, particle2, chargeProd, sigma, epsilon));
+        newIndex = exceptions.size()-1;
+    }
+    exceptionMap[pair<int, int>(particle1, particle2)] = newIndex;
+    return newIndex;
 }
 void NonbondedForce::getExceptionParameters(int index, int& particle1, int& particle2, double& chargeProd, double& sigma, double& epsilon) const {
     particle1 = exceptions[index].particle1;

tests/TestFindExclusions.cpp

@@ -29,15 +29,6 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-/**
- * This tests the createExceptionsFromBonds() method of NonbondedForce, which identifies pairs of atoms
- * whose nonbonded atoms are either excluded or decreased.  The test system is a chain with branches:
- * 
- * 1  3  5  7  9  11  13  15  17  19
- * |  |  |  |  |  |   |   |   |   |
- * 0--2--4--6--8--10--12--14--16--18
- */
-
 #include "AssertionUtilities.h"
 #include "openmm/NonbondedForce.h"
 #include <iostream>
@@ -61,66 +52,124 @@ void addAtomsToExclusions(int atom1, int atom2, vector<set<int> >& exclusions, i
     }
 }
 
-int main() {
-    try {
-        NonbondedForce nonbonded;
-        vector<pair<int, int> > bonds;
-        for (int i = 0; i < NUM_ATOMS; i++)
-            nonbonded.addParticle(1.0, 1.0, 2.0);
-        // 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.createExceptionsFromBonds(bonds, 0.2, 0.4);
-
-        // Build lists of the expected exclusions and 1-4s.
-
-        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);
-        }
-        vector<set<int> > expected14(NUM_ATOMS);
-        int total14 = 0;
-        for (int i = 0; i < NUM_ATOMS; i += 2) {
-            addAtomsToExclusions(i, i+5, expected14, total14);
-            addAtomsToExclusions(i, i+6, expected14, total14);
-            addAtomsToExclusions(i+1, i+3, expected14, total14);
-            addAtomsToExclusions(i+1, i+4, expected14, total14);
+/**
+ * This tests the createExceptionsFromBonds() method of NonbondedForce, which identifies pairs of atoms
+ * whose nonbonded atoms are either excluded or decreased.  The test system is a chain with branches:
+ *
+ * 1  3  5  7  9  11  13  15  17  19
+ * |  |  |  |  |  |   |   |   |   |
+ * 0--2--4--6--8--10--12--14--16--18
+ */
+
+void testFindExceptions() {
+    NonbondedForce nonbonded;
+    vector<pair<int, int> > bonds;
+    for (int i = 0; i < NUM_ATOMS; i++)
+        nonbonded.addParticle(1.0, 1.0, 2.0);
+    // 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.createExceptionsFromBonds(bonds, 0.2, 0.4);
+
+    // Build lists of the expected exclusions and 1-4s.
+
+    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);
+    }
+    vector<set<int> > expected14(NUM_ATOMS);
+    int total14 = 0;
+    for (int i = 0; i < NUM_ATOMS; i += 2) {
+        addAtomsToExclusions(i, i+5, expected14, total14);
+        addAtomsToExclusions(i, i+6, expected14, total14);
+        addAtomsToExclusions(i+1, i+3, expected14, total14);
+        addAtomsToExclusions(i+1, i+4, expected14, total14);
+    }
+
+    // Compare them to the exceptions that were generated.
+
+    ASSERT_EQUAL(totalExclusions+total14, nonbonded.getNumExceptions());
+    for (int i = 0; i < nonbonded.getNumExceptions(); i++) {
+        int particle1, particle2;
+        double chargeProd, sigma, epsilon;
+        nonbonded.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
+        if (chargeProd == 0) {
+            // This is an exclusion.
+
+            ASSERT_EQUAL(0.0, epsilon);
+            ASSERT(expectedExclusions[particle1].find(particle2) != expectedExclusions[particle1].end());
         }
+        else {
+            // This is a 1-4.
 
-        // Compare them to the exceptions that were generated.
-
-        ASSERT_EQUAL(totalExclusions+total14, nonbonded.getNumExceptions());
-        for (int i = 0; i < nonbonded.getNumExceptions(); i++) {
-            int particle1, particle2;
-            double chargeProd, sigma, epsilon;
-            nonbonded.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
-            if (chargeProd == 0) {
-                // This is an exclusion.
-
-                ASSERT_EQUAL(0.0, epsilon);
-                ASSERT(expectedExclusions[particle1].find(particle2) != expectedExclusions[particle1].end());
-            }
-            else {
-                // This is a 1-4.
-
-                ASSERT_EQUAL_TOL(0.2, chargeProd, 1e-10);
-                ASSERT_EQUAL_TOL(1.0, sigma, 1e-10);
-                ASSERT_EQUAL_TOL(0.8, epsilon, 1e-10);
-                ASSERT(expected14[particle1].find(particle2) != expected14[particle1].end());
-            }
+            ASSERT_EQUAL_TOL(0.2, chargeProd, 1e-10);
+            ASSERT_EQUAL_TOL(1.0, sigma, 1e-10);
+            ASSERT_EQUAL_TOL(0.8, epsilon, 1e-10);
+            ASSERT(expected14[particle1].find(particle2) != expected14[particle1].end());
         }
     }
+}
+
+/**
+ * Test replacing existing exclusions.
+ */
+
+void testReplaceExceptions() {
+    NonbondedForce nonbonded;
+    for (int i = 0; i < 10; i++)
+        nonbonded.addParticle(1.0, 1.0, 0.0);
+    nonbonded.addException(0, 1, 1.0, 0.0, 0.0);
+    nonbonded.addException(0, 2, 2.0, 0.0, 0.0);
+    nonbonded.addException(0, 3, 3.0, 0.0, 0.0);
+    try {
+        nonbonded.addException(0, 3, 4.0, 0.0, 0.0);
+        throw std::exception();
+    }
+    catch (const OpenMMException ex) {
+        // This should have thrown an exception.
+    }
+    int p1, p2;
+    double charge, sigma, epsilon;
+    nonbonded.getExceptionParameters(2, p1, p2, charge, sigma, epsilon);
+    ASSERT(p1 == 0);
+    ASSERT(p2 == 3);
+    ASSERT(charge == 3.0);
+    ASSERT(nonbonded.addException(0, 3, 4.0, 0.0, 0.0, true) == 2);
+    nonbonded.getExceptionParameters(2, p1, p2, charge, sigma, epsilon);
+    ASSERT(p1 == 0);
+    ASSERT(p2 == 3);
+    ASSERT(charge == 4.0);
+    ASSERT(nonbonded.addException(1, 3, 4.0, 0.0, 0.0) == 3);
+    try {
+        nonbonded.addException(3, 0, 5.0, 0.0, 0.0);
+        throw std::exception();
+    }
+    catch (const OpenMMException ex) {
+        // This should have thrown an exception.
+    }
+    ASSERT(nonbonded.addException(3, 0, 5.0, 0.0, 0.0, true) == 2);
+    nonbonded.getExceptionParameters(2, p1, p2, charge, sigma, epsilon);
+    ASSERT(p1 == 3);
+    ASSERT(p2 == 0);
+    ASSERT(charge == 5.0);
+}
+
+int main() {
+    try {
+        testFindExceptions();
+        testReplaceExceptions();
+    }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;
         return 1;