(adsbygoogle = window.adsbygoogle || []).push({}); Radial distribution function -- weird result

Dear all,

I think I'm having trouble with my radial distribution function (RDF) calculation and would greatly appreciate it if any one could give me some insight.

I was testing my RDF calculation subroutine to get RDF for certain colloidal crystal configuration (FCC or random distribution). The system I have is a cubic box with length "boxLength_" in the following code. "nParticle_" of colloid are put in the box. Position coordinates of each particle are recorded so I can use "boxLength_", "nParticle_" and positions "x_", "y_" and "z_" as inputs to calculate RDF.

However the first few result I got look like this:

[PLAIN]http://www.pitt.edu/~jeh114/gr_N512_5000steps_cubicStart_sideL800.gif [Broken]

The trouble I felt is that at large r the function doesn't oscillate around 1 as all the literature and textbook show, but some lower value. (Here the curve ends at r = L/2, L the boxLength.)

That could be either my RDF algorithm or my Monte-Carlo algorithm to get the positions. So I tried just randomly spilling all the particles in the box and see what would my RDF be (I expected it to be a constant 1, for all r. Am I right?) and it looks like:

[PLAIN]http://www.pitt.edu/~jeh114/gr_random.gif [Broken]

(again 0 < r < L/2)

And it's around 1 for small r but still decreases for large r. I looked and tested over and over again with my RDF subroutine and couldn't find the bug. I'm attaching it as follows (basically following Frenkel and Smit). Any help would be really appreciated. Thank you very much!

Code (Text):

int RDF::initialize(double length, double diameter, int nPart )

{

boxLength_ = length; cout << "boxLength = " << boxLength_ << endl;

diameter_ = diameter;

delr_ = (boxLength_ / 2.) / 300.;

nParticle_ = nPart;

rho0_ = nParticle_ / pow( boxLength_ ,3);

return 0;

}

int RDF::get_Position( vector<double> x, vector<double> y, vector<double> z )

{

x_ = x; y_ = y; z_ = z;

return 0;

}

int RDF::get_gr(int swicz)

{

switch (swicz){

case 0 : // Initialization

cout << "RDF initialization." << endl;

ngr_ = 0;

nhis_ = NINT ((boxLength_ / 2.) / delr_); cout << "# of bins: " << nhis_ << endl;

for(int i=0; i<=nhis_; i++)

g_.push_back(0.);

break;

case 1 : // Sampling (collecting data through many passing).

ngr_ = ngr_ + 1; cout << "RDF " << ngr_ << "th sampling." << endl;

for( int i = 0; i < nParticle_ - 1; i++){

for( int j = i + 1; j < nParticle_ ; j++){

double xr = x_.at(i) - x_.at(j); double yr = y_.at(i) - y_.at(j); double zr = z_.at(i) - z_.at(j);

xr = xr - boxLength_ * NINT ( xr / boxLength_ ); // PBC

yr = yr - boxLength_ * NINT ( yr / boxLength_ );

zr = zr - boxLength_ * NINT ( zr / boxLength_ );

double r = sqrt( xr*xr + yr*yr + zr*zr );

if( r < (boxLength_ / 2.) ){

int ig = INT (r / delr_);

g_.at(ig) += 2.;

} // if r

//else cout << "Outside!!!!" << endl;

} // for j

} // for i

break;

case 2 : // Calculating g(r)

cout << "RDF calculating." << endl;

for (int i = 0; i < nhis_; i++){

double r = delr_ * (i + 0.5);

double nid = 4. * pi * (r * r) * delr_ * rho0_;

g_.at(i) /= (ngr_ * nParticle_ * nid );

}

break;

default :

cout << "TROUBLE: gr switch not set appropriately!!" << endl;

}

**Physics Forums | Science Articles, Homework Help, Discussion**

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Radial distribution function - weird result

**Physics Forums | Science Articles, Homework Help, Discussion**