Hmm that’s odd, I should only need to #include math.h right?

Maybe I did something wrong then, I will check again.

My code is pretty large out of where it could be going wrong. But So far I have tracked the error down to one function.

```
__device__ double getInteraction(int isoN, int run, double xtemp, double ytemp, double ztemp, double thxt, double thyt, double thzt,
double Isomer1_arr[100][17], double Isomer2_arr[100][17], double SP_arr[100][17], double Shift_arr[100][17]){
//rotation variables
double cx, cy, cz; //cos xyz
double sx, sy, sz; //sin xyz
double xt, yt, zt; //temp/prime
double R11, R12, R13, R21, R22, R23, R31, R32, R33; //temp values
double thx, thy, thz; //shift angle values
double xshift, yshift, zshift; //x,y,z shift values
int iMax, jMax;
double x, y, z; //holds values from isomer 1 or 2
double temp_arr[100][17]; //temp array (replaces sheet 5) 100 rows of 17 columns
//energy variables
double xs, ys, zs; //coordinates of SP, reuse x,y,z from before for coord of isomers
double r, r2; //Inter atome (mol-SP) distance and square of same
double Eq, Eqt; //electrostatic energy
double Enb, Enbt; //non-bonded energy
double A, B, rStar, epsilon; //non-bonded parameters
double rStarS, epsilonS; //non-bonded parameters
int line; //non-bonded parameters
double Et; //total energy
double Ett;
double q, qs; //charges for isomer and SP
//////////////////////ROTATION/////////////////////////
//check if first run or not, if first use shift_arr[][] else use new xyz, txyz from input
//these first values should NOT be set, they should be obtained somehow **fix**
if(run == 1){
thx = -66.25976621 * pi / 180.000;
thy = 120.893644 * pi / 180.000;
thz = 261.891781 * pi / 180.000;
xshift = 0.18333688;
yshift = 9.258853604;
zshift = -0.041139504;
}
else{
thx = thxt * pi / 180.000;
thy = thyt * pi / 180.000;
thz = thzt * pi / 180.000;
xshift = xtemp;
yshift = ytemp;
zshift = ztemp;
}
xt = 0;
yt = 0;
zt = 0;
//do trig
cx = cos(thx);
sx = sin(thx);
cy = cos(thy);
sy = sin(thy);
cz = cos(thz);
sz = sin(thz);
//get matrix elements
R11 = cy * cz;
R12 = cx * sz + sx * sy * cz;
R13 = sx * sz - cx * sy * cz;
R21 = -cy * sz;
R22 = cx * cz - sx * sy * sz;
R23 = sx * cz + cx * sy * sz;
R31 = sy;
R32 = -sx * cy;
R33 = cx * cy;
//get number of atoms in Isomer1 (iMax) and stationary phase (jMax)
if(isoN == 1)
iMax = Isomer1_arr[0][2]; //num rows [0,2] Isomer1_arr ----- iMax = Isomer1_arr_1D[0*100+2];
else
iMax = Isomer2_arr[0][2]; //num rows [0,2] Isomer2_arr ----- iMax = Isomer2_arr_1D[0*100+2];
jMax = 17; //num cols
//copy isomer1 or 2 to temp_arr (sheet 5)
for(int i = 0; i < iMax; i++){ //Y
for(int j = 0; j < jMax; j++){ //X
if(isoN == 1)
temp_arr[i][j] = Isomer1_arr[i][j]; //[i,j]
else
temp_arr[i][j] = Isomer2_arr[i][j]; //[i,j]
}
}
//make transformation
for(int i = 2; i < iMax + 2; i++){
if(isoN == 1){
x = Isomer1_arr[i][2]; //[i,2]
y = Isomer1_arr[i][3]; //[i,3]
z = Isomer1_arr[i][4]; //[i,4]
}
else{
x = Isomer2_arr[i][2]; //[i,2]
y = Isomer2_arr[i][3]; //[i,3]
z = Isomer2_arr[i][4]; //[i,4]
}
xt = x * R11 + y * R12 + z * R13;
yt = x * R21 + y * R22 + z * R23;
zt = x * R31 + y * R32 + z * R33;
temp_arr[i][2] = xt + xshift; //[i,2]
temp_arr[i][3] = yt + yshift; //[i,3]
temp_arr[i][4] = zt + zshift; //[i,4]
}
////////////////////END ROTATION --> ENERGY///////////////////////
//reset x,y,z values to 0.0 for reuse
x = 0.0;
y = 0.0;
z = 0.0;
//get number of atoms in Isomer 1 and SP
iMax = Isomer1_arr[0][2]; //[0,2]
jMax = SP_arr[0][2]; //[0,2]
//initialize energies
Eqt = 0.0; //sum of the electrostatic interactions to date
Enbt = 0.0; //sum of the non-bonded (6-12) interactions to date
Ett = 0.0; //sum of the total energies to date
line = 2;
for(int i = 2; i < iMax + 2; i++){
x = temp_arr[i][2]; //x for isomer //[i + 2, 2]
y = temp_arr[i][3]; //y for isomer //[i + 2, 3]
z = temp_arr[i][4]; //z for isomer //[i + 2, 4]
rStar = temp_arr[i][5]; //rStar for isomer //[i + 2, 5]
epsilon = temp_arr[i][6]; //epsilon for isomer//[i + 2, 6]
q = temp_arr[i][7]; //charge for isomer //[i + 2, 7]
for(int j = 2; j < jMax + 2; j++){
xs = SP_arr[j][2]; //x for SP //[j + 2, 2]
ys = SP_arr[j][3]; //y for SP //[j + 2, 3]
zs = SP_arr[j][4]; //z for SP //[j + 2, 4]
rStarS = SP_arr[j][5]; //rStar for SP //[j + 2, 5]
epsilonS = SP_arr[j][6]; //epsilon for SP //[j + 2, 6]
qs = SP_arr[j][7]; //charge for SP //[j + 2, 7]
//calculate distance from atoms on SP
r2 = ((x - xs) * (x - xs)) + ((y - ys) * (y - ys)) + ((z - zs) * (z - zs));
r = sqrt(r2); //distance between atome j on stationary phase and atom i on isomer
//printf(i);
//printf(j);
//printf(r);
if(r2 == 0.0){
Eq = 1000000.0;
}
else{
Eq = q * qs / r * 0.529177249 * 627.509; //A to Bohr gives hartrees which then needs conversion to kcal/mol
}
Eqt = Eqt + Eq;
//printf(Eq);
//calculate non bonded (6-12) potentioals and ouput
if(r == 0.0){
Enb = 1000000.0;
}
else{
A = pow((rStar + rStarS), 12) * sqrt(epsilon * epsilonS);
B = 2 * pow((rStar + rStarS), 6) * sqrt(epsilon * epsilonS);
Enb = (A / pow(r, 12)) - (B / pow(r, 6));
}
Enbt = Enbt + Enb;
//printf(Enb);
//calculate total interaction and output
Et = Enb + Eq;
Ett = Et + Ett;
//printf(Et);
line = line + 1;
}
}
//printout total energy associated with interaction of isomer with stationary phase
//printf("total");
//printf(Eqt);
//printf(Enbt);
//printf(Ett); //this is the total interaction energy
return Ett;
}
```

The arrays read in data in .csv files to be used. I have checked all the input and it is correct. But somewhere in here there is an error causing a NaN in the end.