Fixed an edge case bug for Gyrcuk Radii. If the radius (ri) of the atom being...

Fixed an edge case bug for Gyrcuk Radii. If the radius (ri) of the atom being descreened is larger than sum of the separation distance (r) and scaled descreening radius (sk), then there is no descreening. This could happen for carbon atom (ri ~= 2) engulfing its bonded hydrogens (sk < 1).

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;