Merge pull request #2187 from peastman/offset

Fixed error in parameter offsets for nonbonded exceptions

platforms/cpu/src/CpuKernels.cpp

@@ -505,15 +505,18 @@ void CpuCalcNonbondedForceKernel::initialize(const System& system, const Nonbond
     numParticles = force.getNumParticles();
     exclusions.resize(numParticles);
     vector<int> nb14s;
+    map<int, int> nb14Index;
     for (int i = 0; i < force.getNumExceptions(); i++) {
         int particle1, particle2;
         double chargeProd, sigma, epsilon;
         force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
         exclusions[particle1].insert(particle2);
         exclusions[particle2].insert(particle1);
-        if (chargeProd != 0.0 || epsilon != 0.0 || exceptionsWithOffsets.find(i) != exceptionsWithOffsets.end())
+        if (chargeProd != 0.0 || epsilon != 0.0 || exceptionsWithOffsets.find(i) != exceptionsWithOffsets.end()) {
+            nb14Index[i] = nb14s.size();
             nb14s.push_back(i);
         }
+    }
 
     // Record the particle parameters.
 
@@ -572,7 +575,7 @@ void CpuCalcNonbondedForceKernel::initialize(const System& system, const Nonbond
         }
         else
             paramIndex = paramPos-paramNames.begin();
-        exceptionParamOffsets[exception].push_back(make_tuple(charge, sigma, epsilon, paramIndex));
+        exceptionParamOffsets[nb14Index[exception]].push_back(make_tuple(charge, sigma, epsilon, paramIndex));
     }
     paramValues.resize(paramNames.size(), 0.0);
 

platforms/cuda/src/CudaKernels.cpp

@@ -1584,14 +1584,17 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
     }
     vector<pair<int, int> > exclusions;
     vector<int> exceptions;
+    map<int, int> exceptionIndex;
     for (int i = 0; i < force.getNumExceptions(); i++) {
         int particle1, particle2;
         double chargeProd, sigma, epsilon;
         force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
         exclusions.push_back(pair<int, int>(particle1, particle2));
-        if (chargeProd != 0.0 || epsilon != 0.0 || exceptionsWithOffsets.find(i) != exceptionsWithOffsets.end())
+        if (chargeProd != 0.0 || epsilon != 0.0 || exceptionsWithOffsets.find(i) != exceptionsWithOffsets.end()) {
+            exceptionIndex[i] = exceptions.size();
             exceptions.push_back(i);
         }
+    }
 
     // Initialize nonbonded interactions.
 
@@ -2033,7 +2036,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
         }
         else
             paramIndex = paramPos-paramNames.begin();
-        exceptionOffsetVec[exception].push_back(make_float4(charge, sigma, epsilon, paramIndex));
+        exceptionOffsetVec[exceptionIndex[exception]].push_back(make_float4(charge, sigma, epsilon, paramIndex));
     }
     paramValues.resize(paramNames.size(), 0.0);
     particleParamOffsets.initialize<float4>(cu, max(force.getNumParticleParameterOffsets(), 1), "particleParamOffsets");

platforms/opencl/src/OpenCLKernels.cpp

@@ -1576,14 +1576,17 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
     }
     vector<pair<int, int> > exclusions;
     vector<int> exceptions;
+    map<int, int> exceptionIndex;
     for (int i = 0; i < force.getNumExceptions(); i++) {
         int particle1, particle2;
         double chargeProd, sigma, epsilon;
         force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
         exclusions.push_back(pair<int, int>(particle1, particle2));
-        if (chargeProd != 0.0 || epsilon != 0.0 || exceptionsWithOffsets.find(i) != exceptionsWithOffsets.end())
+        if (chargeProd != 0.0 || epsilon != 0.0 || exceptionsWithOffsets.find(i) != exceptionsWithOffsets.end()) {
+            exceptionIndex[i] = exceptions.size();
             exceptions.push_back(i);
         }
+    }
 
     // Initialize nonbonded interactions.
 
@@ -1963,7 +1966,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
         }
         else
             paramIndex = paramPos-paramNames.begin();
-        exceptionOffsetVec[exception].push_back(mm_float4(charge, sigma, epsilon, paramIndex));
+        exceptionOffsetVec[exceptionIndex[exception]].push_back(mm_float4(charge, sigma, epsilon, paramIndex));
     }
     paramValues.resize(paramNames.size(), 0.0);
     particleParamOffsets.initialize<mm_float4>(cl, max(force.getNumParticleParameterOffsets(), 1), "particleParamOffsets");

platforms/reference/src/ReferenceKernels.cpp

@@ -877,15 +877,18 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
     numParticles = force.getNumParticles();
     exclusions.resize(numParticles);
     vector<int> nb14s;
+    map<int, int> nb14Index;
     for (int i = 0; i < force.getNumExceptions(); i++) {
         int particle1, particle2;
         double chargeProd, sigma, epsilon;
         force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
         exclusions[particle1].insert(particle2);
         exclusions[particle2].insert(particle1);
-        if (chargeProd != 0.0 || epsilon != 0.0 || exceptionsWithOffsets.find(i) != exceptionsWithOffsets.end())
+        if (chargeProd != 0.0 || epsilon != 0.0 || exceptionsWithOffsets.find(i) != exceptionsWithOffsets.end()) {
+            nb14Index[i] = nb14s.size();
             nb14s.push_back(i);
         }
+    }
 
     // Build the arrays.
 
@@ -916,7 +919,7 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
         int exception;
         double charge, sigma, epsilon;
         force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
-        exceptionParamOffsets[make_pair(param, exception)] = {charge, sigma, epsilon};
+        exceptionParamOffsets[make_pair(param, nb14Index[exception])] = {charge, sigma, epsilon};
     }
     nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
     nonbondedCutoff = force.getCutoffDistance();

tests/TestNonbondedForce.h

@@ -743,8 +743,9 @@ void testParameterOffsets() {
     force->addParticle(1.0, 0.5, 0.6);
     force->addParticle(-1.0, 2.0, 0.7);
     force->addParticle(0.5, 2.0, 0.8);
-    force->addException(0, 1, 1.0, 1.5, 1.0);
+    force->addException(0, 3, 0.0, 1.0, 0.0);
     force->addException(2, 3, 0.5, 1.0, 1.5);
+    force->addException(0, 1, 1.0, 1.5, 1.0);
     force->addGlobalParameter("p1", 0.0);
     force->addGlobalParameter("p2", 1.0);
     force->addParticleParameterOffset("p1", 0, 3.0, 0.5, 0.5);
@@ -775,12 +776,15 @@ void testParameterOffsets() {
             pairSigma[i][j] = 0.5*(particleSigma[i]+particleSigma[j]);
             pairEpsilon[i][j] = sqrt(particleEpsilon[i]*particleEpsilon[j]);
         }
-    pairChargeProd[0][1] = 1.0;
-    pairSigma[0][1] = 1.5;
-    pairEpsilon[0][1] = 1.0;
+    pairChargeProd[0][3] = 0.0;
+    pairSigma[0][3] = 1.0;
+    pairEpsilon[0][3] = 0.0;
     pairChargeProd[2][3] = 0.5+0.5*0.5;
     pairSigma[2][3] = 1.0+0.5*0.5;
     pairEpsilon[2][3] = 1.5+1.5*0.5;
+    pairChargeProd[0][1] = 1.0;
+    pairSigma[0][1] = 1.5;
+    pairEpsilon[0][1] = 1.0;
     
     // Compute the expected energy.