radialdistribution.cpp

Go to the documentation of this file.

/***************************************************************************************************
*****           Copyright (C) 2005  John Schneiderman <JohnMS@member.fsf.org>                  *****
*****                                                                                          *****
*****           This program is free software; you can redistribute it and/or modify           *****
*****           it under the terms of the GNU General Public License as published by           *****
*****           the Free Software Foundation; either version 2 of the License, or              *****
*****           (at your option) any later version.                                            *****
*****                                                                                          *****
*****           This program is distributed in the hope that it will be useful,                *****
*****           but WITHOUT ANY WARRANTY; without even the implied warranty of                 *****
*****           MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                  *****
*****           GNU General Public License for more details.                                   *****
*****                                                                                          *****
*****           You should have received a copy of the GNU General Public License              *****
*****           along with this program; if not, write to the Free Software                    *****
*****           Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA     *****
***************************************************************************************************/
#include "radialdistribution.h"
#include "cluster.h"
#include "configurationdatabase.h"

#include <cmath>

#include <fstream>
using std::ofstream;

#include <iomanip>
using std::ios;
using std::scientific;
using std::showpos;

// class constructor
RadialDistributionFunction::RadialDistributionFunction(const ConfigurationDatabase &conf):BaseHistogram(conf.doubleValue("RadialDistributionFunction", "min"), conf.doubleValue("RadialDistributionFunction", "max"), conf.intValue("RadialDistributionFunction", "bins"))
{
    m_pullingInterval=conf.intValue("RadialDistributionFunction", "pull");
    m_pi=conf.doubleValue("RadialDistributionFunction", "PI");
}

void RadialDistributionFunction::acquire(const Cluster& cluster)
{
    Coordinate ri,rij;
    double del = (m_max - m_min)/m_bins, rijsq;
    ++m_numberHistogramsTaken;

    int n = int(cluster.size());

    //  cout << "entering acquiregr\n";
    //  cout << "acquiring gr, ngr = " << ngr << "\n";

    /* Outer loop begins here */

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

        ri = cluster[i].properties.position;

        /* Inner loop begins here */

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

            rij = ri - cluster[j].properties.position;
            rijsq = rij.x*rij.x + rij.y*rij.y + rij.z*rij.z;

            int nbin=( (int)( (sqrt(rijsq)-m_min)/del ) );
            if(nbin < m_bins && nbin >= 0)
                m_histogramBin.at(nbin)++;
        }

        /* End of inner loop */

    }

    /* End of outer loop */

    //  cout << "leaving acquiregr\n";
}

void RadialDistributionFunction::write(int numAtoms, string base_name)
{
    //  cout << "\n opening output file:\n";
    ofstream out((base_name+".rdf").c_str());
    out.precision(Cluster::PRECISION);
    out.setf(ios::scientific);
    out.setf(ios::showpos);

    double aconst=(4.*m_pi/3.);
    double del = (m_max - m_min)/m_bins;

    //  cout << "delr = " << delr << "\n";

    double rlower,rupper,ravge,rnideal,gr;
    for(int i = 0; i < m_bins; ++i)
    {
        rlower=i*del + m_min;
        rupper=rlower + del;
        ravge=(rupper + rlower)/2.;
        rnideal=aconst*(pow(rupper,3) - pow(rlower,3));
        gr=m_histogramBin.at(i)/( 0.5*numAtoms*rnideal*1.0*m_numberHistogramsTaken );
        out << ravge << " " << gr << "\n";
    }

    out.close();

    //  cout << "\n leaving writegr:\n";
}

---