Merge pull request #1336 from swails/acos

More acos protections in Amoeba

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

@@ -39,15 +39,23 @@ real ym = zcp*xap - xcp*zap;
 real zm = xcp*yap - ycp*xap;
 
 real rm = max(SQRT(xm*xm + ym*ym + zm*zm), (real) 1e-6f);
-real dot = xap*xcp + yap*ycp + zap*zcp;
+real dotp = xap*xcp + yap*ycp + zap*zcp;
 real product = SQRT(rap2*rcp2);
-real cosine = (product > 0 ? (dot/product) : 0);
+real cosine = (product > 0 ? (dotp/product) : 0);
 cosine = max(min(cosine, (real) 1), (real) -1);
-real angle = ACOS(cosine);
+real angle;
+if (cosine > 0.99f || cosine < -0.99f) {
+    real3 cross_prod = cross(make_real3(xap, yap, zap), make_real3(xcp, ycp, zcp));
+    angle = ASIN(SQRT(dot(cross_prod, cross_prod)/(rap2*rcp2)))*RAD_TO_DEG;
+    if (cosine < 0.0f)
+        angle = 180-angle;
+}
+else
+    angle = ACOS(cosine)*RAD_TO_DEG;
 
 // if product == 0, set force/energy to 0
 
-real deltaIdeal = (product > 0 ? (angle*RAD_TO_DEG - angleParams.x) : 0);
+real deltaIdeal = (product > 0 ? (angle - angleParams.x) : 0);
 real deltaIdeal2 = deltaIdeal*deltaIdeal;
 real deltaIdeal3 = deltaIdeal*deltaIdeal2;
 real deltaIdeal4 = deltaIdeal2*deltaIdeal2;
@@ -111,4 +119,4 @@ real dedzid = -dedzia - dedzib - dedzic;
 real3 force1 = make_real3(-dedxia, -dedyia, -dedzia);
 real3 force2 = make_real3(-dedxib, -dedyib, -dedzib);
 real3 force3 = make_real3(-dedxic, -dedyic, -dedzic);
-real3 force4 = make_real3(-dedxid, -dedyid, -dedzid);
\ No newline at end of file
+real3 force4 = make_real3(-dedxid, -dedyid, -dedzid);

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

@@ -17,18 +17,27 @@ real zp = xcb*yab - ycb*xab;
 
 real rp = SQRT(xp*xp + yp*yp + zp*zp);
 
-real dot = xab*xcb + yab*ycb + zab*zcb;
-real cosine = rab*rcb > 0 ? (dot / (rab*rcb)) : (real) 1;
+real dotp = xab*xcb + yab*ycb + zab*zcb;
+real cosine = rab*rcb > 0 ? (dotp / (rab*rcb)) : (real) 1;
 cosine = (cosine > 1 ? (real) 1 : cosine);
 cosine = (cosine < -1 ? -(real) 1 : cosine);
-real angle = ACOS(cosine);
+real angle;
+if (cosine > 0.99f || cosine < -0.99f) {
+    // Highly unlikely a stretch-bend angle will be near 0 or 180, but just in case...
+    real3 cross_prod = cross(make_real3(xab, yab, zab), make_real3(xcb, ycb, zcb));
+    angle = ASIN(SQRT(dot(cross_prod, cross_prod))/(rab*rcb))*RAD_TO_DEG;
+    if (cosine < 0.0f)
+        angle = 180-angle;
+}
+else
+    angle = ACOS(cosine)*RAD_TO_DEG;
 
 // find chain rule terms for the bond angle deviation
 
 float3 parameters = PARAMS[index];
 float2 force_constants = FORCE_CONSTANTS[index];
 
-real dt = RAD_TO_DEG*(angle - parameters.z);
+real dt = angle - RAD_TO_DEG*parameters.z;
 real terma = rab*rp != 0 ? (-RAD_TO_DEG/(rab*rab*rp)) : (real) 0;
 real termc = rcb*rp != 0 ? (RAD_TO_DEG/(rcb*rcb*rp)) : (real) 0;