Radial Distribution Function g(r)
see code at
"gr.c", after complied by "icc -o gr gr.c", uses input .xyz configuration file "configuration.xyz", to calculate g(r) and output results in "gr.dat". Periodic boundary condition is used here for the simulation box. The program asks for on screen the following input parameters:
number of configurations
number of equilibration
particle number N
simulation box size L
number of bins of the histogram
L-J particles
gas and liquid at T = 2.0, rho = 0.1 and 0.8
perfect and equilibrium bcc lattice of 1024 sites at T = 1.0, rho = 1.2
peak location (in unit of first peak) = 1: 2/sqrt(3) : 2sqrt(2)/sqrt(3) : 2 = 1 : 1.15 : 1.63 : 2
perfect and equilibrium fcc lattice of 1372 sites at T = 1.0, rho = 1.2
peak location (in unit of first peak) = 1: sqrt(2) : sqrt(3) : 2 = 1 : 1.41 : 1.73 : 2
Multicomponent version, e.g. AA,BB,CC,AB,BC,CA pairs
***********************
void RadialDis(void)
{
int i,j;
double rij,rij2;
for(i=0;i<NumberOfParticles-1;i++)
for(j=i+1;j<NumberOfParticles;j++)
{
rij2 = Distance(i,j);
rij = sqrt(rij2);
if(rij < L/2.)
{
g[(int)(rij/dradial)] += 2.;
if(identity[i] == identity[j])
{
if(identity[i] == 1) // AA
gAA[(int)(rij/dradial)] += 2.;
else if(identity[i] == 2) //BB
gBB[(int)(rij/dradial)] += 2.;
else // CC
gCC[(int)(rij/dradial)] += 2.;
}// endif same particle
else //AB BC CA
{
if(identity[i]+identity[j] == 3) // AB (1,2)
gAB[(int)(rij/dradial)] += 2.;
else if(identity[i]+identity[j] == 5) // BC (2,3)
gBC[(int)(rij/dradial)] += 2.;
else // CA (3,1)
gCA[(int)(rij/dradial)] += 2.;
}
}//endif r<L/2
}
return;
}
***********************
Normalization (e.g. gCC(r) and gAB(r))
*************************
sprintf(filename,"grCC_%d.dat",JobIndex);
fp=fopen(filename,"w");
for(i=0;i<drBins;i++)
{
gCC[i] /= Count;
gCC[i] /= 4./3.*M_PI*(CUBIC(i+1)-CUBIC(i))*CUBIC(dradial)*rhoC; // rhoC = NC/V
fprintf(fp,"r = %lf\tgCC(r) = %lf\n",(i+0.5)*dradial,gCC[i]/NC);
}
fclose(fp);
sprintf(filename,"grAB_%d.dat",JobIndex);
fp=fopen(filename,"w");
for(i=0;i<drBins;i++)
{
gAB[i] /= Count;
gAB[i] /= 4./3.*M_PI*(CUBIC(i+1)-CUBIC(i))*CUBIC(dradial); // can also use 4*pi*r^2*dr to estimate shell volume
fprintf(fp,"r = %lf\tgAB(r) = %lf\n",(i+0.5)*dradial,gAB[i]*V/(NA*NB)/2.);
}
fclose(fp);
*************************