Fixed bug where angles at 180 degrees would sometimes produce nans

platforms/cuda/src/kernels/kCalculateCustomAngleForces.cu

@@ -112,7 +112,7 @@ void kCalculateCustomAngleForces_kernel()
         float r21 = DOT3(v0, v0);
         float r23 = DOT3(v1, v1);
         float dot = DOT3(v0, v1);
-        float cosine = dot/sqrt(r21*r23);
+        float cosine = max(-1.0f, min(1.0f, dot/sqrt(r21*r23)));
         VARIABLE(0) = acos(cosine);
         VARIABLE(1) = params.x;
         VARIABLE(2) = params.y;

platforms/cuda/src/kernels/kCalculateLocalForces.cu

@@ -198,7 +198,7 @@ void kCalculateLocalForces_kernel()
             float r21               = DOT3(A->v0, A->v0); // dx1 * dx1 + dy1 * dy1 + dz1 * dz1;
             float r23               = DOT3(A->v1, A->v1); // dx2 * dx2 + dy2 * dy2 + dz2 * dz2;
             float dot               = DOT3(A->v0, A->v1); // dx1 * dx2 + dy1 * dy2 + dz1 * dz2;
-            float cosine            = dot / sqrt(r21 * r23);
+            float cosine            = max(-1.0f, min(1.0f, dot / sqrt(r21 * r23)));
 
             float angle_energy;
 /* E */     GETENERGYGIVENANGLECOSINE(cosine, bond_angle, angle_energy);

platforms/opencl/src/kernels/customAngleForce.cl

@@ -22,7 +22,7 @@ __kernel void computeCustomAngleForces(int numAtoms, int numAngles, __global flo
         float r21 = v0.x*v0.x + v0.y*v0.y + v0.z*v0.z;
         float r23 = v1.x*v1.x + v1.y*v1.y + v1.z*v1.z;
         float dot = v0.x*v1.x + v0.y*v1.y + v0.z*v1.z;
-        float cosine = dot/sqrt(r21*r23);
+        float cosine = clamp(dot/sqrt(r21*r23), -1.0f, 1.0f);
         float theta = acos(cosine);
         COMPUTE_FORCE
         float4 c21 = cross(v0, cp)*(dEdAngle/(r21*rp));

platforms/opencl/src/kernels/harmonicAngleForce.cl

@@ -24,7 +24,7 @@ __kernel void calcHarmonicAngleForce(int numAtoms, int numAngles, __global float
         float r21 = v0.x*v0.x + v0.y*v0.y + v0.z*v0.z;
         float r23 = v1.x*v1.x + v1.y*v1.y + v1.z*v1.z;
         float dot = v0.x*v1.x + v0.y*v1.y + v0.z*v1.z;
-        float cosine = dot*RSQRT(r21*r23);
+        float cosine = clamp(dot*RSQRT(r21*r23), -1.0f, 1.0f);
         float deltaIdeal = acos(cosine)-angleParams.x;
         energy += 0.5f*angleParams.y*deltaIdeal*deltaIdeal;
         float dEdR = angleParams.y*deltaIdeal;