Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Radial distribution function - weird result

  1. Jun 10, 2010 #1
    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;
        }
     
    Last edited by a moderator: May 4, 2017
  2. jcsd
  3. Jun 14, 2010 #2
    Re: Radial distribution function -- weird result

    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.
     
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook