Radial distribution function - weird result

[jhsu]

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 anyone 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
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
(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!

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

 

[jhsu]

Sorry I figured it out where the bug was. My NINT function was not working properly for negative numbers, making the distribution dies off 1. Sorry for your time.