Simplified reference platform code for processing nonbonded exclusions

platforms/reference/include/ReferenceCustomHbondIxn.h

@@ -29,6 +29,7 @@
 #include "lepton/ExpressionProgram.h"
 #include "lepton/ParsedExpression.h"
 #include <map>
+#include <set>
 #include <vector>
 
 // ---------------------------------------------------------------------------------------
@@ -143,10 +144,8 @@ class ReferenceCustomHbondIxn : public ReferenceBondIxn {
          @param atomCoordinates    atom coordinates
          @param donorParameters    donor parameters values       donorParameters[donorIndex][parameterIndex]
          @param acceptorParameters acceptor parameters values    acceptorParameters[acceptorIndex][parameterIndex]
-         @param exclusions         exclusion indices             exclusions[donorIndex][acceptorToExcludeIndex]
-                                   exclusions[donorIndex][0] = number of exclusions
-                                   exclusions[donorIndex][no.-1] = indices of acceptors to excluded from
-                                   interacting w/ donor donorIndex
+         @param exclusions         exclusion indices
+                                   exclusions[donorIndex] contains the list of excluded acceptors for that donor
          @param globalParameters   the values of global parameters
          @param forces             force array (forces added)
          @param totalEnergy        total energy
@@ -154,7 +153,7 @@ class ReferenceCustomHbondIxn : public ReferenceBondIxn {
          --------------------------------------------------------------------------------------- */
 
       void calculatePairIxn(std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** donorParameters, RealOpenMM** acceptorParameters,
-                            int** exclusions, const std::map<std::string, double>& globalParameters,
+                            std::vector<std::set<int> >& exclusions, const std::map<std::string, double>& globalParameters,
                             std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy) const;
 
 // ---------------------------------------------------------------------------------------

platforms/reference/include/ReferenceKernels.h

@@ -572,7 +572,7 @@ public:
     void copyParametersToContext(ContextImpl& context, const NonbondedForce& force);
 private:
     int numParticles, num14;
-    int **exclusionArray, **bonded14IndexArray;
+    int **bonded14IndexArray;
     RealOpenMM **particleParamArray, **bonded14ParamArray;
     RealOpenMM nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, dispersionCoefficient;
     int kmax[3], gridSize[3];
@@ -819,7 +819,6 @@ public:
 private:
     int numDonors, numAcceptors, numParticles;
     bool isPeriodic;
-    int **exclusionArray;
     RealOpenMM **donorParamArray, **acceptorParamArray;
     RealOpenMM nonbondedCutoff;
     ReferenceCustomHbondIxn* ixn;

platforms/reference/include/ReferenceLJCoulombIxn.h

@@ -156,10 +156,8 @@ class ReferenceLJCoulombIxn {
          @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 forces           force array (forces added)
          @param energyByAtom     atom energy
@@ -170,7 +168,7 @@ class ReferenceLJCoulombIxn {
          --------------------------------------------------------------------------------------- */
           
       void calculatePairIxn(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
-                            RealOpenMM** atomParameters, int** exclusions,
+                            RealOpenMM** atomParameters, std::vector<std::set<int> >& exclusions,
                             RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
                             RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, bool includeDirect, bool includeReciprocal) const;
 
@@ -182,10 +180,8 @@ private:
          @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 forces           force array (forces added)
          @param energyByAtom     atom energy
@@ -196,7 +192,7 @@ private:
          --------------------------------------------------------------------------------------- */
           
       void calculateEwaldIxn(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
-                            RealOpenMM** atomParameters, int** exclusions,
+                            RealOpenMM** atomParameters, std::vector<std::set<int> >& exclusions,
                             RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
                             RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, bool includeDirect, bool includeReciprocal) const;
 };

platforms/reference/src/ReferenceKernels.cpp

@@ -805,7 +805,6 @@ void ReferenceCalcCustomTorsionForceKernel::copyParametersToContext(ContextImpl&
 
 ReferenceCalcNonbondedForceKernel::~ReferenceCalcNonbondedForceKernel() {
     disposeRealArray(particleParamArray, numParticles);
-    disposeIntArray(exclusionArray, numParticles);
     disposeIntArray(bonded14IndexArray, num14);
     disposeRealArray(bonded14ParamArray, num14);
     if (neighborList != NULL)
@@ -843,14 +842,6 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
         particleParamArray[i][2] = static_cast<RealOpenMM>(charge);
     }
     this->exclusions = exclusions;
-    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;
-    }
     for (int i = 0; i < num14; ++i) {
         int particle1, particle2;
         double charge, radius, depth;
@@ -914,7 +905,7 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
         clj.setUsePME(ewaldAlpha, gridSize);
     if (useSwitchingFunction)
         clj.setUseSwitchingFunction(switchingDistance);
-    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, includeEnergy ? &energy : NULL, includeDirect, includeReciprocal);
+    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusions, 0, forceData, 0, includeEnergy ? &energy : NULL, includeDirect, includeReciprocal);
     if (includeDirect) {
         ReferenceBondForce refBondForce;
         ReferenceLJCoulomb14 nonbonded14;
@@ -1500,7 +1491,6 @@ void ReferenceCalcCustomExternalForceKernel::copyParametersToContext(ContextImpl
 ReferenceCalcCustomHbondForceKernel::~ReferenceCalcCustomHbondForceKernel() {
     disposeRealArray(donorParamArray, numDonors);
     disposeRealArray(acceptorParamArray, numAcceptors);
-    disposeIntArray(exclusionArray, numDonors);
     if (ixn != NULL)
         delete ixn;
 }
@@ -1547,14 +1537,6 @@ void ReferenceCalcCustomHbondForceKernel::initialize(const System& system, const
         for (int j = 0; j < numAcceptorParameters; j++)
             acceptorParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
     }
-    exclusionArray = new int*[numDonors];
-    for (int i = 0; i < numDonors; ++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 nonbondedMethod = CalcCustomHbondForceKernel::NonbondedMethod(force.getNonbondedMethod());
     nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
 
@@ -1603,7 +1585,7 @@ double ReferenceCalcCustomHbondForceKernel::execute(ContextImpl& context, bool i
     map<string, double> globalParameters;
     for (int i = 0; i < (int) globalParameterNames.size(); i++)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
-    ixn->calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusionArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
+    ixn->calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusions, globalParameters, forceData, includeEnergy ? &energy : NULL);
     return energy;
 }
 

platforms/reference/src/SimTKReference/ReferenceCustomHbondIxn.cpp

@@ -34,6 +34,7 @@
 
 using std::map;
 using std::pair;
+using std::set;
 using std::string;
 using std::stringstream;
 using std::vector;
@@ -109,10 +110,8 @@ void ReferenceCustomHbondIxn::setPeriodic(RealVec& boxSize) {
    @param atomCoordinates    atom coordinates
    @param donorParameters    donor parameters values       donorParameters[donorIndex][parameterIndex]
    @param acceptorParameters acceptor parameters values    acceptorParameters[acceptorIndex][parameterIndex]
-   @param exclusions         exclusion indices             exclusions[donorIndex][acceptorToExcludeIndex]
-                             exclusions[donorIndex][0] = number of exclusions
-                             exclusions[donorIndex][no.-1] = indices of acceptors to excluded from
-                             interacting w/ donor donorIndex
+   @param exclusions         exclusion indices
+                             exclusions[donorIndex] contains the list of excluded acceptors for that donor
    @param globalParameters   the values of global parameters
    @param forces             force array (forces added)
    @param totalEnergy        total energy
@@ -120,7 +119,7 @@ void ReferenceCustomHbondIxn::setPeriodic(RealVec& boxSize) {
    --------------------------------------------------------------------------------------- */
 
 void ReferenceCustomHbondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates, RealOpenMM** donorParameters, RealOpenMM** acceptorParameters,
-                                             int** exclusions, const map<string, double>& globalParameters, vector<RealVec>& forces,
+                                             vector<set<int> >& exclusions, const map<string, double>& globalParameters, vector<RealVec>& forces,
                                              RealOpenMM* totalEnergy) const {
 
    map<string, double> variables = globalParameters;
@@ -129,18 +128,8 @@ void ReferenceCustomHbondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates,
 
    int numDonors = donorAtoms.size();
    int numAcceptors = acceptorAtoms.size();
-   int* exclusionIndices = new int[numAcceptors];
-   for( int ii = 0; ii < numAcceptors; ii++ ){
-      exclusionIndices[ii] = -1;
-   }
 
    for( int donor = 0; donor < numDonors; donor++ ){
-
-      // set exclusions
-
-      for (int j = 1; j <= exclusions[donor][0]; j++)
-         exclusionIndices[exclusions[donor][j]] = donor;
-
       // Initialize per-donor parameters.
 
       for (int j = 0; j < (int) donorParamNames.size(); j++)
@@ -149,16 +138,13 @@ void ReferenceCustomHbondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates,
       // loop over atom pairs
 
       for( int acceptor = 0; acceptor < numAcceptors; acceptor++ ){
-
-         if( exclusionIndices[acceptor] != donor ){
+         if (exclusions[donor].find(acceptor) == exclusions[donor].end()) {
              for (int j = 0; j < (int) acceptorParamNames.size(); j++)
                  variables[acceptorParamNames[j]] = acceptorParameters[acceptor][j];
              calculateOneIxn(donor, acceptor, atomCoordinates, variables, forces, totalEnergy);
          }
       }
    }
-
-   delete[] exclusionIndices;
 }
 
   /**---------------------------------------------------------------------------------------

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

@@ -37,6 +37,7 @@
 // make sure that erf() and erfc() are defined.
 #include "openmm/internal/MSVC_erfc.h"
 
+using std::set;
 using std::vector;
 using OpenMM::RealVec;
 
@@ -169,10 +170,8 @@ void ReferenceLJCoulombIxn::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 forces           force array (forces added)
    @param energyByAtom     atom energy
@@ -183,7 +182,7 @@ void ReferenceLJCoulombIxn::setUseSwitchingFunction( RealOpenMM distance ) {
    --------------------------------------------------------------------------------------- */
 
 void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>& atomCoordinates,
-                                             RealOpenMM** atomParameters, int** exclusions,
+                                             RealOpenMM** atomParameters, vector<set<int> >& exclusions,
                                              RealOpenMM* fixedParameters, vector<RealVec>& forces,
                                              RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, bool includeDirect, bool includeReciprocal) const {
     typedef std::complex<RealOpenMM> d_complex;
@@ -423,10 +422,10 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
 
     RealOpenMM totalExclusionEnergy = 0.0f;
     for (int i = 0; i < numberOfAtoms; i++)
-        for (int j = 1; j <= exclusions[i][0]; j++)
-            if (exclusions[i][j] > i) {
+        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter) {
+            if (*iter > i) {
                int ii = i;
-               int jj = exclusions[i][j];
+               int jj = *iter;
 
                RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
                ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR[0] );
@@ -454,6 +453,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
                    energyByAtom[jj] -= realSpaceEwaldEnergy;
                }
             }
+        }
 
     if( totalEnergy )
         *totalEnergy -= totalExclusionEnergy;
@@ -467,10 +467,8 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
    @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 forces           force array (forces added)
    @param energyByAtom     atom energy
@@ -481,7 +479,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
    --------------------------------------------------------------------------------------- */
 
 void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<RealVec>& atomCoordinates,
-                                             RealOpenMM** atomParameters, int** exclusions,
+                                             RealOpenMM** atomParameters, vector<set<int> >& exclusions,
                                              RealOpenMM* fixedParameters, vector<RealVec>& forces,
                                              RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, bool includeDirect, bool includeReciprocal) const {
 
@@ -499,32 +497,13 @@ void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<RealVec>&
        }
    }
    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 ){
-                 calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
-             }
-          }
+          for( int jj = ii+1; jj < numberOfAtoms; jj++ )
+              if (exclusions[jj].find(ii) == exclusions[jj].end())
+                  calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
        }
-
-       delete[] exclusionIndices;
    }
 }