Fixed bug processing exclusions with Ewald (bug 1085)

platforms/cuda/tests/TestCudaEwald.cpp

@@ -53,7 +53,7 @@ using namespace std;
 
 const double TOL = 1e-5;
 
-void testEwaldPME() {
+void testEwaldPME(bool includeExceptions) {
 
 //      Use amorphous NaCl system for the tests
 
@@ -79,6 +79,12 @@ void testEwaldPME() {
         nonbonded->addParticle(1.0, 1.0,0.0);
     for (int i = 0; i < numParticles/2; i++)
         nonbonded->addParticle(-1.0, 1.0,0.0);
+    if (includeExceptions) {
+        // Add some exclusions.
+
+        for (int i = 0; i < numParticles-1; i++)
+          nonbonded->addException(i, i+1, 0.0, 1.0, 0.0);
+    }
     system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     system.addForce(nonbonded);
 
@@ -242,7 +248,8 @@ void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
 
 int main() {
     try {
-     testEwaldPME();
+     testEwaldPME(false);
+     testEwaldPME(true);
 //     testEwald2Ions();
      testErrorTolerance(NonbondedForce::Ewald);
      testErrorTolerance(NonbondedForce::PME);

platforms/opencl/tests/TestOpenCLEwald.cpp

@@ -52,7 +52,7 @@ using namespace std;
 
 const double TOL = 1e-5;
 
-void testEwaldPME() {
+void testEwaldPME(bool includeExceptions) {
 
 //      Use amorphous NaCl system for the tests
 
@@ -78,6 +78,12 @@ void testEwaldPME() {
         nonbonded->addParticle(1.0, 1.0,0.0);
     for (int i = 0; i < numParticles/2; i++)
         nonbonded->addParticle(-1.0, 1.0,0.0);
+    if (includeExceptions) {
+        // Add some exclusions.
+
+        for (int i = 0; i < numParticles-1; i++)
+          nonbonded->addException(i, i+1, 0.0, 1.0, 0.0);
+    }
     system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     system.addForce(nonbonded);
 
@@ -241,7 +247,8 @@ void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
 
 int main() {
     try {
-     testEwaldPME();
+     testEwaldPME(false);
+     testEwaldPME(true);
 //     testEwald2Ions();
      testErrorTolerance(NonbondedForce::Ewald);
      testErrorTolerance(NonbondedForce::PME);

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

@@ -419,6 +419,8 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
                RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
                ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
                RealOpenMM r         = deltaR[0][ReferenceForce::RIndex];
+               if (r >= cutoffDistance)
+                   continue;
                RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
                RealOpenMM alphaR    = alphaEwald * r;
                RealOpenMM dEdR      = (RealOpenMM) (ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR);