mcljsimulation.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 "mcljsimulation.h"
#include "atom.h"
#include "cluster.h"
#include "configurationdatabase.h"
#include "coordinate.h"

#include <fstream>
using std::ifstream;
using std::ofstream;
using std::cout;

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

#include <cmath>

MCLJ_Simulation::MCLJ_Simulation(const ConfigurationDatabase &conf):MC_Simulation(conf)
{
    time_t seconds;

    m_potentialEnergy = 0.0;
    time(&seconds);
    srand((unsigned int) seconds);
}

void MCLJ_Simulation::run(Cluster& cluster, string currentRun, string previousRun)
{
    int numTried = 0, numAccepted = 0;
    double potEngDiff, potentialPart, eDV, ranc, acceptRatio;
    Coordinate distanceToMove, testPosition;
    ofstream runFile, mscFile;
    string currentRunFile = currentRun + ".run", mscFilename = currentRun + ".mis";

    // Open files.
    runFile.open(currentRunFile.c_str());
    runFile.precision(Cluster::PRECISION);
    runFile.setf(ios::scientific);
    runFile.setf(ios::showpos);
    mscFile.open(mscFilename.c_str());
    mscFile.precision(Cluster::PRECISION);
    mscFile.setf(ios::scientific);
    mscFile.setf(ios::showpos);

    // Pre-run operations and manipulations follow.
    m_stepStart = findLastStep(previousRun);
    cluster.calculateCentreMass();
    cluster.calculatePotentialEnergy();
    m_potentialEnergy = cluster.potentialEnergyTotal();

    /* The big outer loop is the "step" loop.  For each step
      we make n attempted moves (n = # atoms).  Hence the double loop.
      We choose the particles sequentially, going through each
      one in turn and trying to move it.
    */
    mscFile << "Number of steps = " << m_runSteps << " at T* = " << m_temperature << endl;
    mscFile << "Current run name = " << currentRun << endl;
    mscFile << "Previous run name = " << previousRun << endl;
    mscFile <<"For the input Cluster-energy (over epsilon, per particle): " << m_potentialEnergy/cluster.size() << endl;
    cout << "Step" << setw(Cluster::WIDTH) << "Potential Energy" << setw(Cluster::WIDTH) << "Accept Ratio" << endl;
    for(int step = 0; step < m_runSteps; ++step)
    {
        //timeStep = step * changeInTime + m_stepStart;
        //Try to move an atom.
        for(int i = 0; i < cluster.size(); ++i)
        {
            numTried++;
            distanceToMove=generateDistanceToMove();
            testPosition = distanceToMove + cluster[i].properties.position;
            potEngDiff=potentialEnergyDifference(cluster, cluster[i], testPosition);
            if (isAcceptableMove(potEngDiff))
            {
                cluster[i].properties.position = testPosition;
                ++numAccepted;
                m_potentialEnergy+=potEngDiff;
                potentialPart = m_potentialEnergy/cluster.size();
            }
            else
            {
                eDV = exp(-potEngDiff/m_temperature);
                ranc = double(rand())/RAND_MAX;
                if (eDV > ranc)
                {
                    cluster[i].properties.position = testPosition;
                    ++numAccepted;
                    m_potentialEnergy+=potEngDiff;
                    potentialPart = m_potentialEnergy/cluster.size();
                }
            }
        }
        potentialPart = m_potentialEnergy/cluster.size();
        cluster.calculatePotentialEnergy();
        acceptRatio = double(numAccepted)/double(numTried);
        m_accumPotentialTotal += potentialPart;
        m_accumPotentialSQ += (potentialPart * potentialPart);
        if ((step % m_radialDensityHistogram.pullingInterval()) == 0)
            m_radialDensityHistogram.acquire(cluster);
        if ((step % m_radialDistributionFunction.pullingInterval()) == 0)
            m_radialDistributionFunction.acquire(cluster);
        if ((step % m_potentialEnergyHistogram.pullingInterval()) == 0)
        {
            cluster.calculatePotentialEnergy();
            m_potentialEnergyHistogram.acquire(cluster);
        }
        if( (step % m_printInterval) == 0)
        {
            runFile << step << setw(Cluster::WIDTH) << potentialPart << setw(Cluster::WIDTH) << acceptRatio << endl;
            cout << step << setw(Cluster::WIDTH) << potentialPart << setw(Cluster::WIDTH) << acceptRatio << endl;
        }
    }
    runFile << m_runSteps << setw(Cluster::WIDTH) << potentialPart << setw(Cluster::WIDTH) << acceptRatio << endl;
    cout << m_runSteps << setw(Cluster::WIDTH) << potentialPart << setw(Cluster::WIDTH) << acceptRatio << endl;
    writeAverages(mscFile, cluster);
    mscFile.close();
    runFile.close();

    //Write out histograms
    m_radialDensityHistogram.write(cluster.size(), currentRun);
    m_radialDistributionFunction.write(cluster.size(), currentRun);
    m_kineticEnergyHistogram.write(cluster.size(), currentRun);
    m_potentialEnergyHistogram.write(cluster.size(), currentRun);
}

double MCLJ_Simulation::findLastStep(string previousRun) const
{
    double timeTemp=0.0, temp=0.0, temp0=0.0;//Temperary values of time.
    ifstream inputFile;//The input run file.
    string filename = previousRun + ".run";

    inputFile.open(filename.c_str());
    while(!inputFile.fail())
        inputFile >> timeTemp >> temp >> temp0;
    inputFile.close();
    return timeTemp;
}

Coordinate MCLJ_Simulation::generateDistanceToMove() const
{
    Coordinate distance;

    distance.x = m_changeInStep * ((double(rand()) / RAND_MAX) - 0.5);
    distance.y = m_changeInStep * ((double(rand()) / RAND_MAX) - 0.5);
    distance.z = m_changeInStep * ((double(rand()) / RAND_MAX) - 0.5);
    return distance;
}

bool MCLJ_Simulation::isAcceptableMove(double potentialEnergyDifference) const
{
    if( potentialEnergyDifference <= 0.0)
        return true;
    else
        return false;
}

double MCLJ_Simulation::potentialEnergyDifference(const Cluster &cluster, const Atom &examinedAtom, const Coordinate &testPosition) const
{
    double potEngOld=0.0, potEngNew=0.0, rsqold, sr2, sr6, sr12, rsqnew;
    Coordinate rijOld, rijNew;

    for(int i=0; i < cluster.size(); ++i)
    {
        if(cluster[i] != examinedAtom)
        {
            //Calculate the potential energy of the Cluster in the old position.
            rijOld = cluster[i].properties.position - examinedAtom.properties.position;
            rsqold = (rijOld.x * rijOld.x) + (rijOld.y * rijOld.y) + (rijOld.z * rijOld.z);
            sr2 = 1.0 / rsqold;
            sr6 = sr2 * sr2 * sr2;
            sr12 = sr6 * sr6;
            potEngOld += (4.0) * (sr12 - sr6);

            //Calculate the potential energy of the Cluster in the new position.
            rijNew = cluster[i].properties.position - testPosition;
            rsqnew = (rijNew.x * rijNew.x) + (rijNew.y * rijNew.y) + (rijNew.z * rijNew.z);
            sr2 = 1.0 / rsqnew;
            sr6 = sr2 * sr2 * sr2;
            sr12= sr6 * sr6;
            potEngNew += (4.0)*(sr12 - sr6);
        }
    }
    return potEngNew - potEngOld;
}

---