Use ASIN instead of ACOS where ACOS may be unstable

This copies the fix for the AmoebaTorsionTorsionForce over to the other two
forces that calculate angles via the inverse cosine -- specifically the in-plane
angle force and the stretch-bend force.  The TorsionTorsion fix did not
completely solve the occasional simulation blowup, although it did make it a
*lot* less frequent (it went from happening within 5 minutes to happening ~once
per day).

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

@@ -43,11 +43,19 @@ real dot = xap*xcp + yap*ycp + zap*zcp;
 real product = SQRT(rap2*rcp2);
 real cosine = (product > 0 ? (dot/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
+    real 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);
+real3 force4 = make_real3(-dedxid, -dedyid, -dedzid);
\ No newline at end of file

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

@@ -21,14 +21,22 @@ real dot = xab*xcb + yab*ycb + zab*zcb;
 real cosine = rab*rcb > 0 ? (dot / (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) {
+    real3 cross_prod = cross(make_real3(xab, yab, zab), make_real3(xcb, ycb, zcb));
+    angle = ASIN(SQRT(dot(cross_prod, cross_prod)/(rap2*rcp2)))*RAD_TO_DEG;
+    if (cosine < 0.0f)
+        angle = 180-angle;
+}
+else
+    real 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;