Bug fixes to interaction groups

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -404,6 +404,14 @@ public:
      * @param set2    the second set of particles forming the interaction group
      */
     void getInteractionGroupParameters(int index, std::set<int>& set1, std::set<int>& set2) const;
+    /**
+     * Set the parameters for an interaction group.
+     * 
+     * @param index   the index of the interaction group for which to set parameters
+     * @param set1    the first set of particles forming the interaction group
+     * @param set2    the second set of particles forming the interaction group
+     */
+    void setInteractionGroupParameters(int index, const std::set<int>& set1, const std::set<int>& set2);
     /**
      * 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.

openmmapi/src/CustomNonbondedForce.cpp

@@ -210,6 +210,12 @@ void CustomNonbondedForce::getInteractionGroupParameters(int index, std::set<int
     set2 = interactionGroups[index].set2;
 }
 
+void CustomNonbondedForce::setInteractionGroupParameters(int index, const std::set<int>& set1, const std::set<int>& set2) {
+    ASSERT_VALID_INDEX(index, interactionGroups);
+    interactionGroups[index].set1 = set1;
+    interactionGroups[index].set2 = set2;
+}
+
 ForceImpl* CustomNonbondedForce::createImpl() const {
     return new CustomNonbondedForceImpl(*this);
 }

platforms/cuda/src/CudaKernels.cpp

@@ -2059,6 +2059,7 @@ void CudaCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNonbo
     
     vector<vector<int> > atomLists;
     vector<pair<int, int> > tiles;
+    map<pair<int, int>, int> duplicateInteractions;
     for (int group = 0; group < force.getNumInteractionGroups(); group++) {
         // Get the list of atoms in this group and sort them.
         
@@ -2100,6 +2101,23 @@ void CudaCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNonbo
                 atoms.push_back(atoms2[j]);
             atomLists.push_back(atoms);
         }
+        
+        // If this group contains duplicate interactions, record that we need to skip them once.
+        
+        for (int i = 0; i < (int) atoms1.size(); i++) {
+            int a1 = atoms1[i];
+            if (set2.find(a1) == set2.end())
+                continue;
+            for (int j = 0; j < (int) atoms2.size() && atoms2[j] < a1; j++) {
+                int a2 = atoms2[j];
+                if (set1.find(a2) != set1.end()) {
+                    pair<int, int> key = make_pair(a2, a1);
+                    if (duplicateInteractions.find(key) == duplicateInteractions.end())
+                        duplicateInteractions[key] = 0;
+                    duplicateInteractions[key]++;
+                }
+            }
+        }
     }
     
     // Build a lookup table for quickly identifying excluded interactions.
@@ -2108,7 +2126,7 @@ void CudaCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNonbo
     for (int i = 0; i < force.getNumExclusions(); i++) {
         int p1, p2;
         force.getExclusionParticles(i, p1, p2);
-        exclusions.insert(make_pair(p1, p2));
+        exclusions.insert(make_pair(min(p1, p2), max(p1, p2)));
     }
     
     // Build the exclusion flags for each tile.  While we're at it, filter out tiles
@@ -2126,13 +2144,23 @@ void CudaCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNonbo
         }
         vector<int>& atoms1 = atomLists[tiles[tile].first];
         vector<int>& atoms2 = atomLists[tiles[tile].second];
-        vector<int> flags(atoms1.size(), (1<<atoms2.size())-1);
+        vector<int> flags(atoms1.size(), (1L<<atoms2.size())-1);
         int numExcluded = 0;
         for (int i = 0; i < (int) atoms1.size(); i++)
             for (int j = 0; j < (int) atoms2.size(); j++) {
                 int a1 = atoms1[i];
                 int a2 = atoms2[j];
-                if (a1 == a2 || exclusions.find(make_pair(a1, a2)) != exclusions.end() || exclusions.find(make_pair(a2, a1)) != exclusions.end()) {
+                bool isExcluded = false;
+                pair<int, int> key = make_pair(min(a1, a2), max(a1, a2));
+                if (a1 == a2 || exclusions.find(key) != exclusions.end())
+                    isExcluded = true; // This is an excluded interaction.
+                else if (duplicateInteractions.find(key) != duplicateInteractions.end() && duplicateInteractions[key] > 0) {
+                    // Both atoms are in both sets, so skip duplicate interactions.
+                    
+                    isExcluded = true;
+                    duplicateInteractions[key]--;
+                }
+                if (isExcluded) {
                     flags[i] &= -1-(1<<j);
                     numExcluded++;
                 }
@@ -2140,8 +2168,7 @@ void CudaCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNonbo
         if (numExcluded == atoms1.size()*atoms2.size())
             continue; // All interactions are excluded.
         tileOrder.push_back(make_pair((int) -atoms2.size(), tile));
-        if (numExcluded > 0)
-            exclusionFlags[tile] = flags;
+        exclusionFlags[tile] = flags;
     }
     sort(tileOrder.begin(), tileOrder.end());
     

platforms/cuda/src/kernels/customNonbondedGroups.cu

@@ -80,10 +80,10 @@ extern "C" __global__ void computeInteractionGroups(
             atomicAdd(&forceBuffers[atom1], static_cast<unsigned long long>((long long) (force.x*0x100000000)));
             atomicAdd(&forceBuffers[atom1+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (force.y*0x100000000)));
             atomicAdd(&forceBuffers[atom1+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (force.z*0x100000000)));
-            atomicAdd(&forceBuffers[atom2], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fx*0x100000000)));
-            atomicAdd(&forceBuffers[atom2+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fy*0x100000000)));
-            atomicAdd(&forceBuffers[atom2+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fz*0x100000000)));
         }
+        atomicAdd(&forceBuffers[atom2], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fx*0x100000000)));
+        atomicAdd(&forceBuffers[atom2+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fy*0x100000000)));
+        atomicAdd(&forceBuffers[atom2+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fz*0x100000000)));
     }
     energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += energy;
 }
\ No newline at end of file

platforms/cuda/tests/TestCudaCustomNonbondedForce.cpp

@@ -582,6 +582,7 @@ void testLargeInteractionGroup() {
     // Create a large system.
     
     System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     for (int i = 0; i < numParticles; i++)
         system.addParticle(1.0);
     CustomNonbondedForce* nonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935456*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
@@ -639,6 +640,18 @@ void testLargeInteractionGroup() {
     // The force on that one particle should be the same.
     
     ASSERT_EQUAL_VEC(state1.getForces()[151], state2.getForces()[151], 1e-4);
+    
+    // Modify the interaction group so it includes all interactions.  This should now reproduce the original forces
+    // on all atoms.
+
+    for (int i = 0; i < numParticles; i++)
+        set1.insert(i);
+    nonbonded->setInteractionGroupParameters(0, set1, set2);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state3 = context.getState(State::Forces);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state3.getForces()[i], 1e-4);
 }
 
 int main(int argc, char* argv[]) {

platforms/opencl/src/OpenCLKernels.cpp

@@ -2070,6 +2070,7 @@ void OpenCLCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNon
     
     vector<vector<int> > atomLists;
     vector<pair<int, int> > tiles;
+    map<pair<int, int>, int> duplicateInteractions;
     for (int group = 0; group < force.getNumInteractionGroups(); group++) {
         // Get the list of atoms in this group and sort them.
         
@@ -2111,6 +2112,23 @@ void OpenCLCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNon
                 atoms.push_back(atoms2[j]);
             atomLists.push_back(atoms);
         }
+        
+        // If this group contains duplicate interactions, record that we need to skip them once.
+        
+        for (int i = 0; i < (int) atoms1.size(); i++) {
+            int a1 = atoms1[i];
+            if (set2.find(a1) == set2.end())
+                continue;
+            for (int j = 0; j < (int) atoms2.size() && atoms2[j] < a1; j++) {
+                int a2 = atoms2[j];
+                if (set1.find(a2) != set1.end()) {
+                    pair<int, int> key = make_pair(a2, a1);
+                    if (duplicateInteractions.find(key) == duplicateInteractions.end())
+                        duplicateInteractions[key] = 0;
+                    duplicateInteractions[key]++;
+                }
+            }
+        }
     }
     
     // Build a lookup table for quickly identifying excluded interactions.
@@ -2119,7 +2137,7 @@ void OpenCLCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNon
     for (int i = 0; i < force.getNumExclusions(); i++) {
         int p1, p2;
         force.getExclusionParticles(i, p1, p2);
-        exclusions.insert(make_pair(p1, p2));
+        exclusions.insert(make_pair(min(p1, p2), max(p1, p2)));
     }
     
     // Build the exclusion flags for each tile.  While we're at it, filter out tiles
@@ -2137,13 +2155,23 @@ void OpenCLCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNon
         }
         vector<int>& atoms1 = atomLists[tiles[tile].first];
         vector<int>& atoms2 = atomLists[tiles[tile].second];
-        vector<int> flags(atoms1.size(), (1<<atoms2.size())-1);
+        vector<int> flags(atoms1.size(), (1L<<atoms2.size())-1);
         int numExcluded = 0;
         for (int i = 0; i < (int) atoms1.size(); i++)
             for (int j = 0; j < (int) atoms2.size(); j++) {
                 int a1 = atoms1[i];
                 int a2 = atoms2[j];
-                if (a1 == a2 || exclusions.find(make_pair(a1, a2)) != exclusions.end() || exclusions.find(make_pair(a2, a1)) != exclusions.end()) {
+                bool isExcluded = false;
+                pair<int, int> key = make_pair(min(a1, a2), max(a1, a2));
+                if (a1 == a2 || exclusions.find(key) != exclusions.end())
+                    isExcluded = true; // This is an excluded interaction.
+                else if (duplicateInteractions.find(key) != duplicateInteractions.end() && duplicateInteractions[key] > 0) {
+                    // Both atoms are in both sets, so skip duplicate interactions.
+                    
+                    isExcluded = true;
+                    duplicateInteractions[key]--;
+                }
+                if (isExcluded) {
                     flags[i] &= -1-(1<<j);
                     numExcluded++;
                 }
@@ -2151,8 +2179,7 @@ void OpenCLCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNon
         if (numExcluded == atoms1.size()*atoms2.size())
             continue; // All interactions are excluded.
         tileOrder.push_back(make_pair((int) -atoms2.size(), tile));
-        if (numExcluded > 0)
-            exclusionFlags[tile] = flags;
+        exclusionFlags[tile] = flags;
     }
     sort(tileOrder.begin(), tileOrder.end());
     

platforms/opencl/src/kernels/customNonbondedGroups.cl

@@ -78,32 +78,34 @@ __kernel void computeInteractionGroups(
         int tj = tgx;
         SYNC_WARPS;
         for (int j = rangeStart; j < rangeEnd; j++) {
-            bool isExcluded = (((exclusions>>tj)&1) == 0);
-            int localIndex = tbx+tj;
-            posq2 = (real4) (localData[localIndex].x, localData[localIndex].y, localData[localIndex].z, localData[localIndex].q);
-            real4 delta = (real4) (posq2.xyz - posq1.xyz, 0);
+            if (tj < rangeEnd) {
+                bool isExcluded = (((exclusions>>tj)&1) == 0);
+                int localIndex = tbx+tj;
+                posq2 = (real4) (localData[localIndex].x, localData[localIndex].y, localData[localIndex].z, localData[localIndex].q);
+                real4 delta = (real4) (posq2.xyz - posq1.xyz, 0);
 #ifdef USE_PERIODIC
-            delta.xyz -= floor(delta.xyz*invPeriodicBoxSize.xyz+0.5f)*periodicBoxSize.xyz;
+                delta.xyz -= floor(delta.xyz*invPeriodicBoxSize.xyz+0.5f)*periodicBoxSize.xyz;
 #endif
-            real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
+                real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 #ifdef USE_CUTOFF
-            if (!isExcluded && r2 < CUTOFF_SQUARED) {
+                if (!isExcluded && r2 < CUTOFF_SQUARED) {
 #endif
-                real invR = RSQRT(r2);
-                real r = RECIP(invR);
-                LOAD_ATOM2_PARAMETERS
-                real dEdR = 0.0f;
-                real tempEnergy = 0.0f;
-                COMPUTE_INTERACTION
-                energy += tempEnergy;
-                delta *= dEdR;
-                force.xyz -= delta.xyz;
-                localData[localIndex].fx += delta.x;
-                localData[localIndex].fy += delta.y;
-                localData[localIndex].fz += delta.z;
+                    real invR = RSQRT(r2);
+                    real r = RECIP(invR);
+                    LOAD_ATOM2_PARAMETERS
+                    real dEdR = 0.0f;
+                    real tempEnergy = 0.0f;
+                    COMPUTE_INTERACTION
+                    energy += tempEnergy;
+                    delta *= dEdR;
+                    force.xyz -= delta.xyz;
+                    localData[localIndex].fx += delta.x;
+                    localData[localIndex].fy += delta.y;
+                    localData[localIndex].fz += delta.z;
 #ifdef USE_CUTOFF
-            }
+                }
 #endif
+            }
             tj = (tj == rangeEnd-1 ? rangeStart : tj+1);
             SYNC_WARPS;
         }
@@ -112,10 +114,10 @@ __kernel void computeInteractionGroups(
             atom_add(&forceBuffers[atom1], (long) (force.x*0x100000000));
             atom_add(&forceBuffers[atom1+PADDED_NUM_ATOMS], (long) (force.y*0x100000000));
             atom_add(&forceBuffers[atom1+2*PADDED_NUM_ATOMS], (long) (force.z*0x100000000));
-            atom_add(&forceBuffers[atom2], (long) (localData[get_local_id(0)].fx*0x100000000));
-            atom_add(&forceBuffers[atom2+PADDED_NUM_ATOMS], (long) (localData[get_local_id(0)].fy*0x100000000));
-            atom_add(&forceBuffers[atom2+2*PADDED_NUM_ATOMS], (long) (localData[get_local_id(0)].fz*0x100000000));
         }
+        atom_add(&forceBuffers[atom2], (long) (localData[get_local_id(0)].fx*0x100000000));
+        atom_add(&forceBuffers[atom2+PADDED_NUM_ATOMS], (long) (localData[get_local_id(0)].fy*0x100000000));
+        atom_add(&forceBuffers[atom2+2*PADDED_NUM_ATOMS], (long) (localData[get_local_id(0)].fz*0x100000000));
 #else
         writeForces(forceBuffers, localData, atom2);
         localData[get_local_id(0)].fx = force.x;

platforms/opencl/tests/TestOpenCLCustomNonbondedForce.cpp

@@ -582,6 +582,7 @@ void testLargeInteractionGroup() {
     // Create a large system.
     
     System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     for (int i = 0; i < numParticles; i++)
         system.addParticle(1.0);
     CustomNonbondedForce* nonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935456*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
@@ -639,6 +640,18 @@ void testLargeInteractionGroup() {
     // The force on that one particle should be the same.
     
     ASSERT_EQUAL_VEC(state1.getForces()[151], state2.getForces()[151], 1e-4);
+    
+    // Modify the interaction group so it includes all interactions.  This should now reproduce the original forces
+    // on all atoms.
+
+    for (int i = 0; i < numParticles; i++)
+        set1.insert(i);
+    nonbonded->setInteractionGroupParameters(0, set1, set2);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state3 = context.getState(State::Forces);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state3.getForces()[i], 1e-4);
 }
 
 int main(int argc, char* argv[]) {

platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp

@@ -177,13 +177,15 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re
             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);
+                    if (*atom1 == *atom2 || exclusions[*atom1].find(*atom2) != exclusions[*atom1].end())
+                        continue; // This is an excluded interaction.
+                    if (*atom1 > *atom2 && set1.find(*atom2) != set1.end() && set2.find(*atom1) != set2.end())
+                        continue; // Both atoms are in both sets, so skip duplicate interactions.
+                    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);
                 }
             }
         }

platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp

@@ -508,6 +508,87 @@ void testInteractionGroups() {
     ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), TOL);
 }
 
+void testLargeInteractionGroup() {
+    const int numMolecules = 300;
+    const int numParticles = numMolecules*2;
+    const double boxSize = 20.0;
+    
+    // Create a large system.
+    
+    ReferencePlatform platform;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    CustomNonbondedForce* nonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935456*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
+    nonbonded->addPerParticleParameter("q");
+    nonbonded->addPerParticleParameter("sigma");
+    nonbonded->addPerParticleParameter("eps");
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<double> params(3);
+    for (int i = 0; i < numMolecules; i++) {
+        if (i < numMolecules/2) {
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.1;
+            nonbonded->addParticle(params);
+            params[0] = -1.0;
+            params[1] = 0.1;
+            nonbonded->addParticle(params);
+        }
+        else {
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.2;
+            nonbonded->addParticle(params);
+            params[0] = -1.0;
+            params[1] = 0.1;
+            nonbonded->addParticle(params);
+        }
+        positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
+        nonbonded->addExclusion(2*i, 2*i+1);
+    }
+    nonbonded->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
+    system.addForce(nonbonded);
+    
+    // Compute the forces.
+    
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    State state1 = context.getState(State::Forces);
+    
+    // Modify the force so only one particle interacts with everything else.
+    
+    set<int> set1, set2;
+    set1.insert(151);
+    for (int i = 0; i < numParticles; i++)
+        set2.insert(i);
+    nonbonded->addInteractionGroup(set1, set2);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state2 = context.getState(State::Forces);
+    
+    // The force on that one particle should be the same.
+    
+    ASSERT_EQUAL_VEC(state1.getForces()[151], state2.getForces()[151], 1e-4);
+    
+    // Modify the interaction group so it includes all interactions.  This should now reproduce the original forces
+    // on all atoms.
+
+    for (int i = 0; i < numParticles; i++)
+        set1.insert(i);
+    nonbonded->setInteractionGroupParameters(0, set1, set2);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state3 = context.getState(State::Forces);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state3.getForces()[i], 1e-4);
+}
+
 int main() {
     try {
         testSimpleExpression();
@@ -520,6 +601,7 @@ int main() {
         testSwitchingFunction();
         testLongRangeCorrection();
         testInteractionGroups();
+        testLargeInteractionGroup();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;