mdljsimulation.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 "mdljsimulation.h"
#include "configurationdatabase.h"

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

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

MDLJ_Simulation::MDLJ_Simulation(const ConfigurationDatabase &conf):MD_Simulation(conf)
{
    m_reducedPotentialEnergy = 0.0, m_reducedKineticEnergy = 0.0, m_reducedEnergyTotal = 0.0;
}

double MDLJ_Simulation::reducedPotentialEnergy() const
{
    return m_reducedPotentialEnergy;
}

double MDLJ_Simulation::reducedKineticEnergy() const
{
    return m_reducedKineticEnergy;
}

double MDLJ_Simulation::reducedTotalEnergy() const
{
    return m_reducedEnergyTotal;
}

void MDLJ_Simulation::run(Cluster &cluster, string currentRun, string prevRun)
{
    ofstream runFile;//This is the run file output.
    ofstream mscFile;//The name of the current miscellaneous file.
    string runFilename = currentRun + ".run";//This is the name of the run file.
    string mscFilename = currentRun + ".mis";//This is the miscellaneous file.

    m_stepStart = findLastStep(prevRun);//This is the starting time.
    //Write initial informaion to the miscellaneous file.
    mscFile.open(mscFilename.c_str());
    mscFile.precision(Cluster::PRECISION);
    mscFile.setf(ios::scientific);
    mscFile.setf(ios::showpos);
    mscFile << "    Number of steps: " << m_runSteps << endl;
    mscFile << "Change in time step: " << m_changeInStep << endl;
    mscFile << "     Print Interval: " << m_printInterval << endl;
    mscFile << "Number of particles: " << cluster.size() << endl;

    //Calculate Cluster information
    cluster.calculateCentreMass();
    cluster.calculatePotentialEnergy();
    cluster.calculateKineticEnergy();
    m_reducedEnergyTotal = (cluster.kineticEnergyTotal() + cluster.potentialEnergyTotal())/cluster.size();
    m_reducedPotentialEnergy = cluster.potentialEnergyTotal()/cluster.size();
    m_reducedKineticEnergy = cluster.kineticEnergyTotal()/cluster.size();
    m_temperature = m_reducedKineticEnergy * (2.0 * cluster.size())/(3.0 * cluster.size() - 6.0);

    //Write Cluster information to the miscellaneous file.
    mscFile << "Energies (over epsilon, per particle): \n";
    mscFile << "     Initial reduced potential energy: " << m_reducedPotentialEnergy << endl;
    mscFile << "       Initial reduced kinetic energy: " << m_reducedKineticEnergy << endl;
    mscFile << "         Initial reduced total energy: " << m_reducedEnergyTotal << endl;
    mscFile << "                  Initial temperature: " << m_temperature << endl;

    //Start run
    runFile.open(runFilename.c_str());
    runFile.precision(Cluster::PRECISION);
    runFile.setf(ios::scientific);
    runFile.setf(ios::showpos);
    m_timeStep = m_stepStart;
    calculateAccelerations(cluster);
    cout << "Time Step" << " " << setw(Cluster::WIDTH) << "Potential" << setw(Cluster::WIDTH) << "Kinetic" << setw(Cluster::WIDTH) << "Energy" << setw(Cluster::WIDTH) << "Temp" << endl;
    for (int intervalStep=0; intervalStep < m_runSteps; intervalStep++)
    {
        m_timeStep = intervalStep * m_changeInStep + m_stepStart;
        if ((intervalStep%m_printInterval)==0)
        {
            runFile << m_timeStep << " " << setw(Cluster::WIDTH) << m_reducedPotentialEnergy << setw(Cluster::WIDTH)<< m_reducedKineticEnergy <<  setw(Cluster::WIDTH) << m_reducedEnergyTotal << setw(Cluster::WIDTH)  << m_temperature << endl;
            cout<< m_timeStep << " " << setw(Cluster::WIDTH) << m_reducedPotentialEnergy << setw(Cluster::WIDTH)<< m_reducedKineticEnergy <<  setw(Cluster::WIDTH) << m_reducedEnergyTotal << setw(Cluster::WIDTH)  << m_temperature << endl;
        }
        calculatePositions(cluster);
        calculateAccelerations(cluster);
        calculateVelocities(cluster);
        m_reducedEnergyTotal = (cluster.kineticEnergyTotal() + cluster.potentialEnergyTotal())/cluster.size();
        m_reducedPotentialEnergy = cluster.potentialEnergyTotal()/cluster.size();
        m_reducedKineticEnergy = cluster.kineticEnergyTotal()/cluster.size();
        m_temperature = m_reducedKineticEnergy * (2.0 * cluster.size())/(3.0 * cluster.size() - 6.0);
        m_accumPotentialTotal += m_reducedPotentialEnergy;
        m_accumKineticTotal += m_reducedKineticEnergy;
        m_accumEnergyTotal += m_reducedEnergyTotal;
        m_accumPotentialSQ += (m_reducedPotentialEnergy*m_reducedPotentialEnergy);
        m_accumKineticSQ += (m_reducedKineticEnergy*m_reducedKineticEnergy);
        m_accumEnergySQ += (m_reducedEnergyTotal*m_reducedEnergyTotal);
        if ((intervalStep % m_radialDensityHistogram.pullingInterval()) == 0)
            m_radialDensityHistogram.acquire(cluster);
        if ((intervalStep % m_radialDistributionFunction.pullingInterval()) == 0)
            m_radialDistributionFunction.acquire(cluster);
        if ((intervalStep % m_potentialEnergyHistogram.pullingInterval()) == 0)
        {
            cluster.calculatePotentialEnergy();
            m_potentialEnergyHistogram.acquire(cluster);
        }
        if ((intervalStep % m_kineticEnergyHistogram.pullingInterval()) == 0)
        {
            cluster.calculateKineticEnergy();
            m_kineticEnergyHistogram.acquire(cluster);
        }
    }
    m_timeStep = m_runSteps * m_changeInStep + m_stepStart;
    runFile << m_timeStep << " " << setw(Cluster::WIDTH) << m_reducedPotentialEnergy << setw(Cluster::WIDTH) << m_reducedKineticEnergy <<  setw(Cluster::WIDTH) << m_reducedEnergyTotal << setw(Cluster::WIDTH)  << m_temperature << endl;
    cout << m_timeStep << " " << setw(Cluster::WIDTH) << m_reducedPotentialEnergy << setw(Cluster::WIDTH) << m_reducedKineticEnergy <<  setw(Cluster::WIDTH) << m_reducedEnergyTotal << setw(Cluster::WIDTH)  << m_temperature << endl;
    runFile.close();

    //Write out average energies
    writeAverages(mscFile, cluster);
    mscFile.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);
}

void MDLJ_Simulation::clear()
{
    MD_Simulation::clear();
    m_reducedPotentialEnergy = 0.0, m_reducedKineticEnergy = 0.0, m_reducedEnergyTotal = 0.0;
}

void MDLJ_Simulation::calculatePositions(Cluster &cluster) const
{
    double changeInTimeSquared = 0.0;//The change in the time squared.

    changeInTimeSquared = m_changeInStep*m_changeInStep;
    for(int i = 0; i < cluster.size(); i++)
    {
        cluster[i].properties.position.x += (m_changeInStep * cluster[i].properties.velocity.x + 0.5 * changeInTimeSquared * cluster[i].properties.acceleration.x);
        cluster[i].properties.position.y += (m_changeInStep * cluster[i].properties.velocity.y + 0.5 * changeInTimeSquared * cluster[i].properties.acceleration.y);
        cluster[i].properties.position.z += (m_changeInStep * cluster[i].properties.velocity.z + 0.5 * changeInTimeSquared * cluster[i].properties.acceleration.z);
        cluster[i].properties.velocity.x += (0.5 * m_changeInStep * cluster[i].properties.acceleration.x);
        cluster[i].properties.velocity.y += (0.5 * m_changeInStep * cluster[i].properties.acceleration.y);
        cluster[i].properties.velocity.z += (0.5 * m_changeInStep * cluster[i].properties.acceleration.z);
    }
}

void MDLJ_Simulation::calculateVelocities(Cluster &cluster) const
{
    double sumPotentialSQ = 0.0;//The sum of the potential energy squared.
    int size = cluster.size();//The size of the Cluster.

    cluster.setKineticEnergyTotal(0.5 * sumPotentialSQ);
    for (int i = 0; i < size; i++)
    {
        cluster[i].properties.velocity.x += (0.5 * m_changeInStep * cluster[i].properties.acceleration.x);
        cluster[i].properties.velocity.y += (0.5 * m_changeInStep * cluster[i].properties.acceleration.y);
        cluster[i].properties.velocity.z += (0.5 * m_changeInStep * cluster[i].properties.acceleration.z);
        sumPotentialSQ += (cluster[i].properties.velocity.x * cluster[i].properties.velocity.x +  cluster[i].properties.velocity.y * cluster[i].properties.velocity.y + cluster[i].properties.velocity.z * cluster[i].properties.velocity.z);
    }
    cluster.setKineticEnergyTotal(0.5 * sumPotentialSQ);
}

void MDLJ_Simulation::calculateAccelerations(Cluster &cluster) const
{
    Coordinate examAtomPos, examAtomAcc, neighbourAtomAcc, adjustedAccAtom;
    int size = cluster.size();//The size of the Cluster.
    //Lennard-Jones modelling numbers.
    double sr2 = 0.0, sr6 = 0.0, sr12 = 0.0, fij = 0.0;
    double rijsq = 0.0, atomPotentialEnergy = 0.0, wij = 0.0;

    cluster.setPotentialEnergyTotal(0.0);
    for(int i = 0; i < size; i++)
        cluster[i].properties.acceleration=0.0;
    for(int i=0; i < size-1; i++)
    {
        examAtomPos = cluster[i].properties.position;
        examAtomAcc = cluster[i].properties.acceleration;
        for(int j = i+1; j < size; j++)
        {
            neighbourAtomAcc = examAtomPos - cluster[j].properties.position;
            rijsq = (neighbourAtomAcc.x * neighbourAtomAcc.x) + (neighbourAtomAcc.y * neighbourAtomAcc.y) + (neighbourAtomAcc.z * neighbourAtomAcc.z);
            sr2 = 1.0/rijsq;
            sr6 = sr2*sr2*sr2;
            sr12 = sr6*sr6;
            atomPotentialEnergy = sr12-sr6;
            cluster.setPotentialEnergyTotal(cluster.potentialEnergyTotal() + atomPotentialEnergy);
            wij = atomPotentialEnergy + sr12;
            fij = wij * sr2;
            adjustedAccAtom = Coordinate(fij) * neighbourAtomAcc;
            examAtomAcc += adjustedAccAtom;
            cluster[j].properties.acceleration -= adjustedAccAtom;
        }
        cluster[i].properties.acceleration = examAtomAcc;
    }
    cluster.setPotentialEnergyTotal(cluster.potentialEnergyTotal()*4.0);
    for (int i = 0; i < size; i++)
        cluster[i].properties.acceleration *= 24.0;
}

void MDLJ_Simulation::adjustPosCenMass(Cluster& cluster) const
{
    Coordinate centre = 0.0;//The positions of the centre of mass.

    for (int i = 0; i < cluster.size(); i++)
        centre += cluster[i].properties.position;
    centre /= cluster.size();
    for (int i = 0; i < cluster.size(); i++)
        cluster[i].properties.position -= (centre + cluster.positionCentreMass());
}

double MDLJ_Simulation::findLastStep(string prevRun) const
{
    double timeTemp=0.0, temp=0.0, temp0=0.0, temp1=0.0, temp2=0.0;//Temperary values of time.
    ifstream inputFile;//The input run file.
    string filename = "";//The name of the last run file.

    filename = prevRun + ".run";
    inputFile.open(filename.c_str());
    while(!inputFile.fail())
        inputFile >> timeTemp >> temp >> temp0 >> temp1 >> temp2;
    inputFile.close();
    return timeTemp;
}

---