cluster.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 "cluster.h"

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

//Output file includes
using std::ifstream;
using std::ofstream;
using std::endl;

Cluster::Cluster():m_atoms(), m_positionCentreMass(), m_velocityCentreMass()
{
    m_potentialEnergyTotal = 0.0, m_kineticEnergyTotal = 0.0;
}

Cluster::Cluster(const Cluster &cluster)
{
    *this = cluster;
}

double Cluster::potentialEnergyTotal() const
{
    return m_potentialEnergyTotal;
}

double Cluster::kineticEnergyTotal() const
{
    return m_kineticEnergyTotal;
}

Coordinate Cluster::positionCentreMass() const
{
    return m_positionCentreMass;
}

Coordinate Cluster::velocityCentreMass() const
{
    return m_velocityCentreMass;
}

void Cluster::setPotentialEnergyTotal(double potentialEnergyTotal)
{
    m_potentialEnergyTotal = potentialEnergyTotal;
}

void Cluster::setKineticEnergyTotal(double kineticEnergyTotal)
{
    m_kineticEnergyTotal = kineticEnergyTotal;
}

void Cluster::setPositionCentreMass(Coordinate positionCentreMass)
{
    m_positionCentreMass = positionCentreMass;
}

void Cluster::setVelocityCentreMass(Coordinate velocityCentreMass)
{
    m_velocityCentreMass = velocityCentreMass;
}

void Cluster::calculatePotentialEnergy()
{
    Coordinate examAtomPos;
    Coordinate examAtomsDistance;
    //Parts of calculating the Lennard-Jones energies.
    double sqrAtomDistance=0.0, sr2=0.0, sr6=0.0, sr12=0.0;
    //The potential energy between two m_atoms.
    double atomPotentialEnergy=0.0;
    //The partial sums of the potential energy.
    double potentialEnergyPart=0.0;
    //The average of the potential energy.
    double potentialEnergyAverage=0.0;

    for (int i = 0; i < int(m_atoms.size()); i++)
        m_atoms[i].properties.potentialEnergy=0.0;
    for (int i = 0; i < int(m_atoms.size()); i++)
    {
        potentialEnergyPart = 0.0;
        examAtomPos = m_atoms[i].properties.position;
        for (int j = 0; j < int(m_atoms.size()); j++)
            if (j != i)
            {
                examAtomsDistance = examAtomPos-m_atoms[j].properties.position;
                sqrAtomDistance = (examAtomsDistance.x * examAtomsDistance.x) + (examAtomsDistance.y * examAtomsDistance.y) + (examAtomsDistance.z * examAtomsDistance.z);
                sr2 = 1.0/sqrAtomDistance;
                sr6 = sr2*sr2*sr2;
                sr12 = sr6*sr6;
                atomPotentialEnergy = sr12-sr6;
                potentialEnergyPart += (4.0*atomPotentialEnergy);
            }
        m_atoms[i].properties.potentialEnergy = potentialEnergyPart/2;
        potentialEnergyAverage += m_atoms[i].properties.potentialEnergy;
    }
    m_potentialEnergyTotal = potentialEnergyAverage;
}

void Cluster::calculateKineticEnergy()
{
    double kineticEnergyAverage=0.0, kineticEnergyPart=0.0;

    for (int i = 0; i < int(m_atoms.size()); i++)
    {
        kineticEnergyPart = 0.5*((m_atoms[i].properties.velocity.x * m_atoms[i].properties.velocity.x) + (m_atoms[i].properties.velocity.y * m_atoms[i].properties.velocity.y) + (m_atoms[i].properties.velocity.z * m_atoms[i].properties.velocity.z));
        m_atoms[i].properties.kineticEnergy = kineticEnergyPart;
        kineticEnergyAverage += kineticEnergyPart;
    }
    m_kineticEnergyTotal = kineticEnergyAverage;
}

void Cluster::calculateCentreMass()
{
    int size = int(m_atoms.size());//The size of the Cluster.

    m_positionCentreMass=0.0;
    m_velocityCentreMass=0.0;
    for (int i = 0; i < size; i++)
    {
        m_positionCentreMass += m_atoms[i].properties.position;
        m_velocityCentreMass += m_atoms[i].properties.velocity;
    }
    m_positionCentreMass /= size;
    m_velocityCentreMass /= size;
}

void Cluster::read(string baseName)
{
    ifstream positionFile, velocityFile;//The position and velocity files.
    string positionFilename, velocityFilename;//The names of the position and velocity files.
    Atom newAtom;//The values of the atom to placed into the Cluster.

    positionFilename = baseName + ".pos";
    velocityFilename = baseName + ".vel";
    positionFile.open(positionFilename.c_str());
    velocityFile.open(velocityFilename.c_str());
    if (positionFile.fail() || velocityFile.fail())
        throw("position or velocity file not found!");
    m_atoms.clear();
    positionFile >> newAtom.properties.position;
    velocityFile >> newAtom.properties.velocity;
    while((!positionFile.fail()) && (!velocityFile.fail()))
    {
        m_atoms.push_back(newAtom);
        positionFile >> newAtom.properties.position;
        velocityFile >> newAtom.properties.velocity;
    }
    positionFile.close();
    velocityFile.close();
}

void Cluster::write(string baseName)
{
    ofstream outputFile;//This is the file to write to.
    string filename;//This is the name of the position and velocity files.

    filename = baseName + ".pos";
    outputFile.open(filename.c_str());
    outputFile.precision(PRECISION);
    outputFile.setf(ios::scientific);
    outputFile.setf(ios::showpos);
    for (int i = 0; i < int(m_atoms.size()); i++)
        outputFile << m_atoms[i].properties.position << endl;
    if (outputFile.fail())
        throw("Could not write position file.");
    outputFile.close();
    outputFile.clear();
    filename = baseName + ".vel";
    outputFile.open(filename.c_str());
    outputFile.precision(PRECISION);
    outputFile.setf(ios::scientific);
    outputFile.setf(ios::showpos);
    for (int i = 0; i < int(m_atoms.size()); i++)
        outputFile << m_atoms[i].properties.velocity << endl;
    if (outputFile.fail())
        throw("Could not write velocity file.");
    outputFile.close();
}

void Cluster::clear()
{
    m_atoms.clear();
    m_positionCentreMass = 0.0, m_velocityCentreMass= 0.0, m_potentialEnergyTotal = 0.0, m_kineticEnergyTotal = 0.0;
}

int Cluster::size() const
{
    return int(m_atoms.size());
}

Cluster& Cluster::operator=(const Cluster &rhs)
{
    if (this != &rhs)
    {
        this->clear();
        for (int i = 0; i < rhs.size(); i++)
            m_atoms.push_back(rhs[i]);
        m_positionCentreMass = rhs.positionCentreMass();
        m_velocityCentreMass = rhs.velocityCentreMass();
        m_potentialEnergyTotal = rhs.potentialEnergyTotal();
        m_kineticEnergyTotal = rhs.kineticEnergyTotal();
    }
    return *this;
}

Atom& Cluster::operator[](int i)
{
    return m_atoms[i];
}

const Atom& Cluster::operator[](int i) const
{
    return m_atoms[i];
}

void Cluster::addAtom(const Atom &atom)
{
    m_atoms.push_back(atom);
}

ostream& operator<<(ostream &outStream, const Cluster &cluster)
{
    outStream.precision(Cluster::PRECISION);
    outStream.setf(ios::scientific);
    outStream.setf(ios::showpos);
    outStream << "Positions:" << endl;
    for (int i = 0; i < cluster.size(); i++)
        outStream << cluster[i].properties.position << endl;
    outStream << "Velocities:" << endl;
    for (int i = 0; i < cluster.size(); i++)
        outStream << cluster[i].properties.velocity << endl;
    outStream << "Acceleration:" << endl;
    for (int i = 0; i < cluster.size(); i++)
        outStream << cluster[i].properties.acceleration << endl;
    outStream << "Potential Energies:" << endl;
    for (int i = 0; i < cluster.size(); i++)
        outStream << cluster[i].properties.potentialEnergy << endl;
    outStream << "Kinetic Energies:" << endl;
    for (int i = 0; i < cluster.size(); i++)
        outStream << cluster[i].properties.kineticEnergy;
    return outStream;
}

bool operator==(const Cluster &lhs, const Cluster &rhs)
{
    if (lhs.size() != rhs.size())
        return false;
    if (lhs.m_potentialEnergyTotal != rhs.m_potentialEnergyTotal)
        return false;
    if (lhs.m_kineticEnergyTotal != rhs.m_kineticEnergyTotal)
        return false;
    if (lhs.m_positionCentreMass != rhs.m_positionCentreMass)
        return false;
    if (lhs.m_velocityCentreMass != rhs.m_velocityCentreMass)
        return false;
    for (int i = 0; i < lhs.size(); i++)
    {
        if (lhs[i].properties.position != rhs[i].properties.position)
            return false;
        else if (lhs[i].properties.velocity != rhs[i].properties.velocity)
            return false;
        else if (lhs[i].properties.acceleration != rhs[i].properties.acceleration)
            return false;
        else if (lhs[i].properties.kineticEnergy != rhs[i].properties.kineticEnergy)
            return false;
        else if (lhs[i].properties.potentialEnergy != rhs[i].properties.potentialEnergy)
            return false;
    }
    return true;
}

---