Merge pull request #110 from peastman/master

Fixed numeric problems using PME with Drude particles

platforms/cuda/src/CudaContext.cpp

@@ -224,6 +224,8 @@ CudaContext::CudaContext(const System& system, int deviceIndex, bool useBlocking
     compilationDefines["ACOS"] = useDoublePrecision ? "acos" : "acosf";
     compilationDefines["ASIN"] = useDoublePrecision ? "asin" : "asinf";
     compilationDefines["ATAN"] = useDoublePrecision ? "atan" : "atanf";
+    compilationDefines["ERF"] = useDoublePrecision ? "erf" : "erff";
+    compilationDefines["ERFC"] = useDoublePrecision ? "erfc" : "erfcf";
     
     // Create the work thread used for parallelization when running on multiple devices.
     

platforms/cuda/src/kernels/coulombLennardJones.cu

@@ -23,8 +23,9 @@ if ((!isExcluded && r2 < CUTOFF_SQUARED) || needCorrection) {
         // Subtract off the part of this interaction that was included in the reciprocal space contribution.
 
         if (1-erfcAlphaR > 1e-6) {
-            tempForce = -prefactor*((1.0f-erfcAlphaR)-alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
-            tempEnergy += -prefactor*(1.0f-erfcAlphaR);
+            real erfAlphaR = ERF(alphaR); // Our erfc approximation is not accurate enough when r is very small, which happens with Drude particles.
+            tempForce = -prefactor*(erfAlphaR-alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
+            tempEnergy += -prefactor*erfAlphaR;
         }
     }
     else {

platforms/opencl/src/kernels/coulombLennardJones.cl

@@ -23,8 +23,9 @@ if ((!isExcluded && r2 < CUTOFF_SQUARED) || needCorrection) {
         // Subtract off the part of this interaction that was included in the reciprocal space contribution.
 
         if (1-erfcAlphaR > 1e-6) {
-            tempForce = -prefactor*((1.0f-erfcAlphaR)-alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
-            tempEnergy += -prefactor*(1.0f-erfcAlphaR);
+            real erfAlphaR = erf(alphaR); // Our erfc approximation is not accurate enough when r is very small, which happens with Drude particles.
+            tempForce = -prefactor*(erfAlphaR-alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
+            tempEnergy += -prefactor*erfAlphaR;
         }
     }
     else {