-1. #IND00 appearing without division?

Question

While I know that -1.#IND00 means I am trying to divide something by zero, my program gave out -1.#IND00 without me doing any sort of division.

I am trying to make a code which records the distance of particles in a lattice, and I told the program to ignore the cases where it is trying to subtract the position of the same particle by ignoring that iteration and move on to the next. However, with the formula:

comparative_distance[par_num][current_par]=sqrt((position[current_par][0]-position[par_num][0])*(position[current_par][0]-position[par_num][0])+(position[current_par][0]-position[par_num][0])*(position[current_par][1]-position[par_num][1])+(position[current_par][2]-position[par_num][2])*(position[current_par][2]-position[par_num][2]));, 

which is simply sqrt(x^2+y^2+z^2), it gave me -1.#IND00 for several results, and I think we cannot ignore these cases since this information will be carried on to calculate inter-particle potentials and forces (when adding the potentials up, I think the computer will not know what to do when it needs to add an integer with a -1.#IND00 result).

Maybe I am misinformed in what -1.#IND00 means, but what does the output want to say about -1.#IND00?

Here's the whole code if you need to take a look (it is still not finished):

#include <stdlib.h>
#include <math.h>
#define Kb 8.6173324e-5 //Boltzmann constant (in eV/K)
#define PI 3.14159265
#define C 299792.458  // in angstrom/ps
#define temp 298 // in K
#define mass 3.72112e10 //in eV/c^2
#define upsilon 1.03e-2 //in eV
#define rho 3.4//in Angstrom

//Initialization
double Box_Muller_Transformationx(double box1) //this is the algorithm for Box-Muller transformation
{
    double vx = (sqrt(-2*log(box1))*cos(2*PI*box1))*sqrt(Kb*temp/mass);
    return vx;
}
double Box_Muller_Transformationy(double box2) //this is the algorithm for Box-Muller transformation
{
    double vy = (sqrt(-2*log(box2))*cos(2*PI*box2))*sqrt(Kb*temp/mass);
    return vy;
}
double Box_Muller_Transformationz(double box3) //this is the algorithm for Box-Muller transformation
{
    double vz = (sqrt(-2*log(box3))*cos(2*PI*box3))*sqrt(Kb*temp/mass);
    return vz;
}

int par_num = 125, b = 3, a=0, coordinates=3, current_par=124; //this is the lattice size, boundary is set as 10 in each direction
double t=0, box1, box2, box3;
int main()
{
    double i=0, j=0, k=0, vx=0, vy=0, vz=0; //generating the 125 particles as a 3D lattice with random velocities
    double velocity[125][3];
    double position[125][3];
    double potential[125][125];
    double comparative_distance[125][125];
    double v_new[125][3];
    double force[125][3];
    double rho_term[125][125];
    double each_total_V[par_num][1];
    double vx_sum;
    double vy_sum;
    double vz_sum;
    double vx_final;
    double vy_final;
    double vz_final;
    double r;
    //finding the altered velocities after initializing velocities
    for (par_num=0;par_num<125;par_num++)
    {
       vx_sum+=velocity[par_num][0];
       vy_sum+=velocity[par_num][1];
       vz_sum+=velocity[par_num][2];
    }
    vx_final=vx_sum/125;
    vy_final=vy_sum/125;
    vz_final=vz_sum/125;

    FILE *fp;
    fp = fopen ("lab2testing.txt", "w+");
    fprintf(fp, "%d \n%s \n", 125, "particle");

    for (i=0, par_num=0; par_num<125; i++, par_num++)
    {
            double box1 = (rand()/ (double)RAND_MAX);
            double box2 = (rand()/ (double)RAND_MAX);
            double box3 = (rand()/ (double)RAND_MAX);
            position[par_num][0]=i; //position components, i=-1 because if we want the particle to be placed at 0 again, but the big for-loop will generate i=1 unless we tune it down.
            position[par_num][1]=j;
            position[par_num][2]=k;
            double vx = Box_Muller_Transformationx(box1); //velocity components
            double vy = Box_Muller_Transformationy(box2);
            double vz = Box_Muller_Transformationz(box3);
            velocity[par_num][0]=vx;
            velocity[par_num][1]=vy;
            velocity[par_num][2]=vz;
            v_new[par_num][0]=velocity[par_num][0]-vx_final;
            v_new[par_num][1]=velocity[par_num][1]-vy_final;
            v_new[par_num][2]=velocity[par_num][2]-vz_final;
            if(par_num%5==4)
            {
                i=-1;
                j++;
            }

            if(par_num%25==24)
            {
                i=-1;
                j=0;
                k++;
            }

            if(par_num==125)
            {
                position[par_num][2]=5;
            }
            fprintf(fp, "%s %lf %lf %lf\n ","Ar", position[par_num][0], position[par_num][1], position[par_num][2]);
            printf("%s %lf %lf %lf\n ","Ar", position[par_num][0], position[par_num][1], position[par_num][2]);
    }

//Iteration......
double p,q;
        fprintf(fp, "%d \n%s \n", 125, "particle");
        for (coordinates=0; coordinates<3; coordinates++) //coordinates loop
        {
            for (par_num=0; par_num<125; par_num++)  //calculating for each particle
            {
                for (current_par=0; current_par<125; current_par++) //counting from 1st until 125th particle
                {
                    comparative_distance[par_num][current_par]=sqrt((position[current_par][0]-position[par_num][0])*(position[current_par][0]-position[par_num][0])+(position[current_par][0]-position[par_num][0])*(position[current_par][1]-position[par_num][1])+(position[current_par][2]-position[par_num][2])*(position[current_par][2]-position[par_num][2]));
                    if(comparative_distance[par_num][current_par]==0 || par_num==current_par)
                    {
                    continue;
                    }
                    rho_term[par_num][current_par] = rho/comparative_distance[par_num][current_par];
                    p = rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par];
                    q = rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par]*rho_term[par_num][current_par];
                    potential[par_num][current_par]=4*upsilon*(q-p);
                    printf("%lf\n", comparative_distance[par_num][current_par]);
                }
            }
        }


    return(0);
}

Answer 1

The printed value -1.#IND00 typically indicates the result of a floating point operation is either infinity or not-a-number. In this case it happens to be the latter. Where is comes from is the call to sqrt. With some additional formatting the problem becomes more obvious:

comparative_distance[par_num][current_par]=sqrt(
  (position[current_par][0]-position[par_num][0])*(position[current_par][0]-position[par_num][0])+
  (position[current_par][0]-position[par_num][0])*(position[current_par][1]-position[par_num][1])+
  (position[current_par][2]-position[par_num][2])*(position[current_par][2]-position[par_num][2]));

The second term is incorrect. Instead of y^2 you have x*y. Change it to:

comparative_distance[par_num][current_par]=sqrt(
  (position[current_par][0]-position[par_num][0])*(position[current_par][0]-position[par_num][0])+
  (position[current_par][1]-position[par_num][1])*(position[current_par][1]-position[par_num][1])+
  (position[current_par][2]-position[par_num][2])*(position[current_par][2]-position[par_num][2]));