RDF

//
// rdf.C
//
// Calculate radial distribution functions (RDF).
//
// Written by: Yanting Wang             December 12, 2009
//

#include <iostream>
#include <string>
using namespace std;

#include <stdlib.h>
#include <libgen.h>
#include <string.h>
#include <math.h>

#include "vector.h"

const double PI = 2.0 * acos(0.0);

//
// Periodic boundary condition (PBC)
//

Vector pbc( Vector r, double L )
{
//
// PBC correction
//

double hL = L / 2.0;

if( r.x > hL )  r.x -= L;
if( r.x < -hL ) r.x += L;

if( r.y > hL )  r.y -= L;
if( r.y < -hL ) r.y += L;

if( r.z > hL )  r.z -= L;
if( r.z < -hL ) r.z += L;

//
// Check if the vector is now inside the box
//

if( r.x > hL || r.x < -hL )
{
cerr << "r.x = " << r.x << " is out of simulation box." << endl;

exit( -1 );
}

if( r.y > hL || r.y < -hL )
{
cerr << "r.y = " << r.y << " is out of simulation box." << endl;

exit( -1 );
}

if( r.z > hL || r.z < -hL )
{
cerr << "r.z = " << r.z << " is out of simulation box." << endl;

exit( -1 );
}

return r;
}

int main( int argc, char *argv[] )
{
//
// Parse input arguments
//

if( argc != 3 )
{
cerr << "Usage: " << basename( argv[0] ) << " MC/MD ng" << endl;
return -1;
}

bool md = false;
if( !strcmp( argv[1], "MD" ) ) md = true;

const int ng = atoi( argv[2] );         // number of bins for RDF

int g[ng];    // array of RDF

int natom;       // number of atoms

double L;        // simulation box size
double width;    // bin size of RDF

int nconf = 0;   // number of configurations.

for( int i=0; i< ng; ++i ) g[i] = 0;

while( cin >> natom )
{
++ nconf;

//
// Read in the configuration
//

cin >> L;

width = L / 2.0 / ng;

Vector r[natom], v, a, pa;

string name;

for( int i=0; i<natom; ++i )
{
cin >> name >> r[i];

if( md ) cin >> v >> a >> pa;
}

//
// Count pairs
//

for( int i=0; i<natom-1; ++i )
{
for( int j=i+1; j<natom; ++j )
{
double d = pbc( r[i]-r[j], L ).r();

int index = d / width;

if( index < ng ) g[ index ] += 2;      // one pair contributes twice
}
}
}

//
// Normalize
//

const double rho = natom / L / L / L;       // density

const double c = nconf * natom * rho * 4.0 / 3.0 * PI * pow( width, 3.0 );

for( int i=0; i<ng; ++i )
{
cout << width * (i+0.5) << " "
<< g[i] / ( c * ( pow( i+1, 3.0) - pow( i, 3.0 ) ) )
<< endl;
}

return 0;
}