Merge pull request #1990 from mjschnie/master

Fixed an edge case bug for Gyrcuk Radii. 

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

@@ -56,6 +56,10 @@ __device__ void computeBornSumOneInteraction(AtomData1& atom1, AtomData1& atom2)
     real r2 = dot(delta, delta);
     real r = SQRT(r2);
     float sk = atom2.scaledRadius;
+
+    if (atom1.radius > r + sk)
+        return; // No descreening due to atom1 engulfing atom2.
+
     real sk2 = sk*sk;
     if (atom1.radius+r < sk) {
         real lik = atom1.radius;
@@ -423,6 +427,9 @@ __device__ void computeBornChainRuleInteraction(AtomData3& atom1, AtomData3& ato
     real r = SQRT(r2);
     real de = 0;
 
+    if (atom1.radius > r + sk)
+        return; // No descreening due to atom1 engulfing atom2.
+
     if (atom1.radius+r < sk) {
         real uik = sk-r;
         real uik4 = uik*uik;

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceGeneralizedKirkwoodForce.cpp

@@ -157,6 +157,10 @@ void AmoebaReferenceGeneralizedKirkwoodForce::calculateGrycukBornRadii(const vec
             double r        = sqrt(r2);
 
             double sk       = _atomicRadii[jj]*_scaleFactors[jj];
+
+            // If atom ii engulfs the descreening atom, then continue.
+            if (_atomicRadii[ii] > r + sk) continue;
+
             double sk2      = sk*sk;
 
             if ((_atomicRadii[ii] + r) < sk) {

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp

@@ -4151,6 +4151,9 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateGrycukChainRuleP
     double r                        = sqrt(r2);
     double de                       = 0.0;
 
+    // If atom index engulfs the descreening atom, then there is no descreening. 
+    if (_atomicRadii[iIndex] > r + sk) return;
+
     if ((_atomicRadii[iIndex] + r) < sk) {
         double uik4;
         uik                = sk - r;