Minor optimization to validating exclusions (#4453)

* Minor optimization to validating exclusions

* Optimizations to findMoleculeGroups()

openmmapi/src/CustomGBForceImpl.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -29,20 +29,19 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
+#ifdef WIN32
+  #define _USE_MATH_DEFINES // Needed to get M_PI
+#endif
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/CustomGBForceImpl.h"
 #include "openmm/internal/Messages.h"
 #include "openmm/kernels.h"
+#include <cmath>
 #include <sstream>
 
 using namespace OpenMM;
-using std::map;
-using std::pair;
-using std::vector;
-using std::set;
-using std::string;
-using std::stringstream;
+using namespace std;
 
 CustomGBForceImpl::CustomGBForceImpl(const CustomGBForce& owner) : owner(owner) {
     forceGroup = owner.getForceGroup();
@@ -74,6 +73,8 @@ void CustomGBForceImpl::initialize(ContextImpl& context) {
     for (int i = 0; i < owner.getNumExclusions(); i++) {
         int particle1, particle2;
         owner.getExclusionParticles(i, particle1, particle2);
+        int minp = min(particle1, particle2);
+        int maxp = max(particle1, particle2);
         if (particle1 < 0 || particle1 >= owner.getNumParticles()) {
             stringstream msg;
             msg << "CustomGBForce: Illegal particle index for an exclusion: ";
@@ -86,7 +87,7 @@ void CustomGBForceImpl::initialize(ContextImpl& context) {
             msg << particle2;
             throw OpenMMException(msg.str());
         }
-        if (exclusions[particle1].count(particle2) > 0 || exclusions[particle2].count(particle1) > 0) {
+        if (exclusions[minp].count(maxp) > 0) {
             stringstream msg;
             msg << "CustomGBForce: Multiple exclusions are specified for particles ";
             msg << particle1;
@@ -94,8 +95,7 @@ void CustomGBForceImpl::initialize(ContextImpl& context) {
             msg << particle2;
             throw OpenMMException(msg.str());
         }
-        exclusions[particle1].insert(particle2);
-        exclusions[particle2].insert(particle1);
+        exclusions[minp].insert(maxp);
     }
     if (owner.getNonbondedMethod() == CustomGBForce::CutoffPeriodic) {
         Vec3 boxVectors[3];

openmmapi/src/CustomNonbondedForceImpl.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2022 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -84,6 +84,8 @@ void CustomNonbondedForceImpl::initialize(ContextImpl& context) {
     for (int i = 0; i < owner.getNumExclusions(); i++) {
         int particle1, particle2;
         owner.getExclusionParticles(i, particle1, particle2);
+        int minp = min(particle1, particle2);
+        int maxp = max(particle1, particle2);
         if (particle1 < 0 || particle1 >= owner.getNumParticles()) {
             stringstream msg;
             msg << "CustomNonbondedForce: Illegal particle index for an exclusion: ";
@@ -96,7 +98,7 @@ void CustomNonbondedForceImpl::initialize(ContextImpl& context) {
             msg << particle2;
             throw OpenMMException(msg.str());
         }
-        if (exclusions[particle1].count(particle2) > 0 || exclusions[particle2].count(particle1) > 0) {
+        if (exclusions[minp].count(maxp) > 0) {
             stringstream msg;
             msg << "CustomNonbondedForce: Multiple exclusions are specified for particles ";
             msg << particle1;
@@ -104,8 +106,7 @@ void CustomNonbondedForceImpl::initialize(ContextImpl& context) {
             msg << particle2;
             throw OpenMMException(msg.str());
         }
-        exclusions[particle1].insert(particle2);
-        exclusions[particle2].insert(particle1);
+        exclusions[minp].insert(maxp);
     }
     if (owner.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic) {
         Vec3 boxVectors[3];

openmmapi/src/NonbondedForceImpl.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -81,6 +81,8 @@ void NonbondedForceImpl::initialize(ContextImpl& context) {
         int particle[2];
         double chargeProd, sigma, epsilon;
         owner.getExceptionParameters(i, particle[0], particle[1], chargeProd, sigma, epsilon);
+        int minp = min(particle[0], particle[1]);
+        int maxp = max(particle[0], particle[1]);
         for (int j = 0; j < 2; j++) {
             if (particle[j] < 0 || particle[j] >= owner.getNumParticles()) {
                 stringstream msg;
@@ -89,7 +91,7 @@ void NonbondedForceImpl::initialize(ContextImpl& context) {
                 throw OpenMMException(msg.str());
             }
         }
-        if (exceptions[particle[0]].count(particle[1]) > 0 || exceptions[particle[1]].count(particle[0]) > 0) {
+        if (exceptions[minp].count(maxp) > 0) {
             stringstream msg;
             msg << "NonbondedForce: Multiple exceptions are specified for particles ";
             msg << particle[0];
@@ -97,8 +99,7 @@ void NonbondedForceImpl::initialize(ContextImpl& context) {
             msg << particle[1];
             throw OpenMMException(msg.str());
         }
-        exceptions[particle[0]].insert(particle[1]);
-        exceptions[particle[1]].insert(particle[0]);
+        exceptions[minp].insert(maxp);
         if (sigma < 0)
             throw OpenMMException("NonbondedForce: sigma for an exception cannot be negative");
         if (epsilon < 0)

platforms/common/src/ComputeContext.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2019-2022 Stanford University and the Authors.      *
+ * Portions copyright (c) 2019-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -210,13 +210,14 @@ void ComputeContext::findMoleculeGroups() {
             atomBonds[particle2].push_back(particle1);
         }
         for (auto force : forces) {
+            vector<int> particles;
             for (int j = 0; j < force->getNumParticleGroups(); j++) {
-                vector<int> particles;
                 force->getParticlesInGroup(j, particles);
-                for (int k = 0; k < (int) particles.size(); k++)
-                    for (int m = 0; m < (int) particles.size(); m++)
-                        if (k != m)
-                            atomBonds[particles[k]].push_back(particles[m]);
+                for (int k = 1; k < (int) particles.size(); k++)
+                    for (int m = 0; m < k; m++) {
+                        atomBonds[particles[k]].push_back(particles[m]);
+                        atomBonds[particles[m]].push_back(particles[k]);
+                    }
             }
         }
 
@@ -242,13 +243,14 @@ void ComputeContext::findMoleculeGroups() {
             system.getConstraintParameters(i, particle1, particle2, distance);
             molecules[atomMolecule[particle1]].constraints.push_back(i);
         }
-        for (int i = 0; i < (int) forces.size(); i++)
+        for (int i = 0; i < (int) forces.size(); i++) {
+            vector<int> particles;
             for (int j = 0; j < forces[i]->getNumParticleGroups(); j++) {
-                vector<int> particles;
                 forces[i]->getParticlesInGroup(j, particles);
                 if (particles.size() > 0)
                     molecules[atomMolecule[particles[0]]].groups[i].push_back(j);
             }
+        }
     }
 
     // Sort them into groups of identical molecules.
@@ -293,10 +295,10 @@ void ComputeContext::findMoleculeGroups() {
             for (int i = 0; i < (int) forces.size() && identical; i++) {
                 if (mol.groups[i].size() != mol2.groups[i].size())
                     identical = false;
+                vector<int> p1, p2;
                 for (int k = 0; k < (int) mol.groups[i].size() && identical; k++) {
                     if (!forces[i]->areGroupsIdentical(mol.groups[i][k], mol2.groups[i][k]))
                         identical = false;
-                    vector<int> p1, p2;
                     forces[i]->getParticlesInGroup(mol.groups[i][k], p1);
                     forces[i]->getParticlesInGroup(mol2.groups[i][k], p2);
                     for (int m = 0; m < p1.size(); m++)