Simple fix to combining rules allowing AMOEBA VdW parameters to be zero. Zero...

Simple fix to combining rules allowing AMOEBA VdW parameters to be zero.  Zero parameter excludes all interactions involving that particle, in agreement with TINKER.

plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp

@@ -179,7 +179,11 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
             } else {
                 double iSigma2 = iSigma*iSigma;
                 double jSigma2 = jSigma*jSigma;
-                sigma = 2.0f * (iSigma2 * iSigma + jSigma2 * jSigma) / (iSigma2 + jSigma2);
+		if ((iSigma2 + jSigma2) != 0.0) {
+		  sigma = 2.0f * (iSigma2 * iSigma + jSigma2 * jSigma) / (iSigma2 + jSigma2);
+		} else {
+		  sigma = 0.0;
+		}
             }
             // ARITHMETIC = 1
             // GEOMETRIC  = 2
@@ -190,11 +194,19 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
             } else if (epsilonCombiningRule == "GEOMETRIC") {
                 epsilon = std::sqrt(iEpsilon * jEpsilon);
             } else if (epsilonCombiningRule == "HARMONIC") {
-                epsilon = 2.0f * (iEpsilon * jEpsilon) / (iEpsilon + jEpsilon);
+	      if ((iEpsilon + jEpsilon) != 0.0) {
+		  epsilon = 2.0f * (iEpsilon * jEpsilon) / (iEpsilon + jEpsilon);
+		} else {
+		epsilon = 0.0;
+	      }
             } else {
-                double epsilonS = std::sqrt(iEpsilon) + std::sqrt(jEpsilon);
-                epsilon = 4.0f * (iEpsilon * jEpsilon) / (epsilonS * epsilonS);
-            }
+	      double epsilonS = std::sqrt(iEpsilon) + std::sqrt(jEpsilon);
+	      if (epsilonS != 0.0) {
+		epsilon = 4.0f * (iEpsilon * jEpsilon) / (epsilonS * epsilonS);
+	      } else {
+		epsilon = 0.0;
+	      }
+	    }
             int count = class1->second * class2->second;
             // Below is an exact copy of stuff from the previous block.
             double rv = sigma;

plugins/amoeba/platforms/cuda/src/kernels/amoebaVdwForce2.cu

@@ -12,17 +12,28 @@
 #else
     real sigma1_2 = sigmaEpsilon1.x*sigmaEpsilon1.x;
     real sigma2_2 = sigmaEpsilon2.x*sigmaEpsilon2.x;
-    real sigma = 2*(sigmaEpsilon1.x*sigma1_2 + sigmaEpsilon2.x*sigma2_2)/(sigma1_2+sigma2_2);
+    real sigmasum = sigma1_2+sigma2_2;
+    real sigma = 0.0f;
+    if (sigmasum != 0.0f) {
+        sigma = 2*(sigmaEpsilon1.x*sigma1_2 + sigmaEpsilon2.x*sigma2_2)/(sigma1_2+sigma2_2);
+    }
 #endif
 #if EPSILON_COMBINING_RULE == 1
     real epsilon = 0.5f*(sigmaEpsilon1.y + sigmaEpsilon2.y);
 #elif EPSILON_COMBINING_RULE == 2
     real epsilon = SQRT(sigmaEpsilon1.y*sigmaEpsilon2.y);
 #elif EPSILON_COMBINING_RULE == 3
-    real epsilon = 2*(sigmaEpsilon1.y*sigmaEpsilon2.y)/(sigmaEpsilon1.y+sigmaEpsilon2.y);
+    real epssum = sigmaEpsilon1.y+sigmaEpsilon2.y;
+    real epsilon = 0.0f;
+    if (epssum != 0.0f) {
+        real epsilon = 2*(sigmaEpsilon1.y*sigmaEpsilon2.y)/(sigmaEpsilon1.y+sigmaEpsilon2.y);
+    }
 #else
     real epsilon_s = SQRT(sigmaEpsilon1.y) + SQRT(sigmaEpsilon2.y);
-    real epsilon = 4*sigmaEpsilon1.y*sigmaEpsilon2.y/(epsilon_s*epsilon_s);
+    real epsilon = 0.0f;
+    if (epsilon_s != 0.0f) {
+        epsilon = 4*sigmaEpsilon1.y*sigmaEpsilon2.y/(epsilon_s*epsilon_s);
+    }
 #endif
     real r6 = r2*r2*r2;
     real r7 = r6*r;