Fixed wrapping of dihedral potentials and changed how singularities are handled

src/TabulatedAngle.cu

@@ -59,7 +59,7 @@ TabulatedAnglePotential::TabulatedAnglePotential(String fileName) : fileName(fil
 		pot[size++] = atof(tokenList[1].val());
 	}
 	// units "1/deg" "1/radian"  *57.29578
-	angle_step_inv = 57.29578f / (angle[1]-angle[0]); 
+	angle_step_inv = 57.29578f * (size-1) / (angle[size-1]-angle[0]); 
 	delete[] angle;
 	fclose(inp);
 }

src/TabulatedDihedral.cu

@@ -2,6 +2,8 @@
 // Authors: Justin Dufresne and Terrance Howard, 2013
 
 #include "TabulatedDihedral.h"
+#include <cassert>
+#define BD_PI 3.1415927f
 
 TabulatedDihedralPotential::TabulatedDihedralPotential() :
 		pot(NULL), size(0), fileName("") {}
@@ -53,7 +55,24 @@ TabulatedDihedralPotential::TabulatedDihedralPotential(String fileName) : fileNa
 		pot[size++] = atof(tokenList[1].val());
 	}
 	// units "1/deg" "1/radian"  *57.29578
-	angle_step_inv = 57.29578f / (angle[1]-angle[0]); 
+	float deltaAngle = (angle[size-1]-angle[0])/(size-1); 
+	assert( deltaAngle > 0 );
+	assert( size*deltaAngle >= 360 );
+
+	float tmp[size];
+	for (int i = 0; i < size; ++i) {
+	    // j=0 corresponsds to angle[i] in [-Pi,-Pi+delta)
+	    float a = (angle[i] + 180.0f);
+	    while (a < 0) a += 360.0f;
+	    while (a >= 360) a -= 360.0f;
+	    int j = floor( a / deltaAngle );
+	    if (j >= size) continue;
+	    tmp[j] = pot[i];
+	}
+	for (int i = 0; i < size; ++i) pot[i] = tmp[i];
+
+	angle_step_inv = 57.29578f / deltaAngle;
+		 
 	delete[] angle;
 	fclose(inp);
 }

src/TabulatedMethods.cuh

@@ -126,8 +126,16 @@ __device__ inline void computeDihedral(const TabulatedDihedralPotential* __restr
 	force = -dU * d->angle_step_inv;
 
 	// avoid singularity when one angle is straight 
-	force = (crossABC.rLength() > 1.0f || crossBCD.rLength() > 1.0f) ? 0.0f : force;
+	force = (ab.length2()*bc.length2()*crossABC.rLength2() > 100.0f || (ab.length2()*bc.length2()*crossABC.rLength2())*crossBCD.rLength2() > 100.0f) ? 0.0f : force;
+
+	// if (force >= 10000.0f)
+	//     printf("pot[%d] = %f; pot[%d] = %f\n", home,U0, home1, U0+dU);
+	if ( force > 1000.0f ) 
+	    force = 1000.0f;
+	if ( force < -1000.0f ) 
+	    force = -1000.0f;
 	assert( force < 10000.0f );
+
 	f1 *= force;
 	f2 *= force;
 	f3 *= force;