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;
}