/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2010 University of Bonn. All rights reserved. * Please see the LICENSE file or "Copyright notice" in builder.cpp for details. */ /** \file atom.cpp * * Function implementations for the class atom. * */ // include config.h #ifdef HAVE_CONFIG_H #include #endif #include "CodePatterns/MemDebug.hpp" #include "atom.hpp" #include "Bond/bond.hpp" #include "CodePatterns/Log.hpp" #include "config.hpp" #include "Element/element.hpp" #include "LinearAlgebra/Vector.hpp" #include "World.hpp" #include "molecule.hpp" #include "Shapes/Shape.hpp" #include #include /************************************* Functions for class atom *************************************/ atom::atom() : father(this), sort(&Nr), mol(0) {}; atom::atom(atom *pointer) : ParticleInfo(pointer), father(pointer), sort(&Nr), mol(0) { setType(pointer->getType()); // copy element of atom AtomicPosition = pointer->AtomicPosition; // copy trajectory of coordination AtomicVelocity = pointer->AtomicVelocity; // copy trajectory of velocity AtomicForce = pointer->AtomicForce; setFixedIon(pointer->getFixedIon()); }; atom *atom::clone(){ atom *res = new atom(this); World::getInstance().registerAtom(res); return res; } /** Destructor of class atom. */ atom::~atom() { removeFromMolecule(); }; void atom::UpdateSteps() { LOG(4,"atom::UpdateSteps() called."); // append to position, velocity and force vector AtomInfo::AppendTrajectoryStep(); // append to ListOfBonds vector BondedParticleInfo::AppendTrajectoryStep(); } atom *atom::GetTrueFather() { if(father == this){ // top most father is the one that points on itself return this; } else if(!father) { return 0; } else { return father->GetTrueFather(); } }; /** Sets father to itself or its father in case of copying a molecule. */ void atom::CorrectFather() { if (father->father != father) // same atom in copy's father points to itself // father = this; // set father to itself (copy of a whole molecule) // else father = father->father; // set father to original's father }; void atom::EqualsFather ( const atom *ptr, const atom **res ) const { if ( ptr == father ) *res = this; }; bool atom::isFather(const atom *ptr){ return ptr==father; } bool atom::IsInShape(const Shape& shape) const { return shape.isInside(getPosition()); }; bool atom::OutputIndexed(ofstream * const out, const int ElementNo, const int AtomNo, const char *comment) const { if (out != NULL) { *out << "Ion_Type" << ElementNo << "_" << AtomNo << "\t" << fixed << setprecision(9) << showpoint; *out << at(0) << "\t" << at(1) << "\t" << at(2); *out << "\t" << (int)(getFixedIon()); if (getAtomicVelocity().Norm() > MYEPSILON) *out << "\t" << scientific << setprecision(6) << getAtomicVelocity()[0] << "\t" << getAtomicVelocity()[1] << "\t" << getAtomicVelocity()[2] << "\t"; if (comment != NULL) *out << " # " << comment << endl; else *out << " # molecule nr " << getNr() << endl; return true; } else return false; }; bool atom::OutputArrayIndexed(ostream * const out,const enumeration &elementLookup, int *AtomNo, const char *comment) const { AtomNo[getType()->getAtomicNumber()]++; // increment number if (out != NULL) { const element *elemental = getType(); ASSERT(elementLookup.there.find(elemental)!=elementLookup.there.end(),"Type of this atom was not in the formula upon enumeration"); *out << "Ion_Type" << elementLookup.there.find(elemental)->second << "_" << AtomNo[elemental->getAtomicNumber()] << "\t" << fixed << setprecision(9) << showpoint; *out << at(0) << "\t" << at(1) << "\t" << at(2); *out << "\t" << getFixedIon(); if (getAtomicVelocity().Norm() > MYEPSILON) *out << "\t" << scientific << setprecision(6) << getAtomicVelocity()[0] << "\t" << getAtomicVelocity()[1] << "\t" << getAtomicVelocity()[2] << "\t"; if (comment != NULL) *out << " # " << comment << endl; else *out << " # molecule nr " << getNr() << endl; return true; } else return false; }; bool atom::OutputXYZLine(ofstream *out) const { if (out != NULL) { *out << getType()->getSymbol() << "\t" << at(0) << "\t" << at(1) << "\t" << at(2) << "\t" << endl; return true; } else return false; }; bool atom::OutputTrajectory(ofstream * const out, const enumeration &elementLookup, int *AtomNo, const int step) const { AtomNo[getType()->getAtomicNumber()]++; if (out != NULL) { const element *elemental = getType(); ASSERT(elementLookup.there.find(elemental)!=elementLookup.there.end(),"Type of this atom was not in the formula upon enumeration"); *out << "Ion_Type" << elementLookup.there.find(elemental)->second << "_" << AtomNo[getType()->getAtomicNumber()] << "\t" << fixed << setprecision(9) << showpoint; *out << getPositionAtStep(step)[0] << "\t" << getPositionAtStep(step)[1] << "\t" << getPositionAtStep(step)[2]; *out << "\t" << (int)(getFixedIon()); if (getAtomicVelocityAtStep(step).Norm() > MYEPSILON) *out << "\t" << scientific << setprecision(6) << getAtomicVelocityAtStep(step)[0] << "\t" << getAtomicVelocityAtStep(step)[1] << "\t" << getAtomicVelocityAtStep(step)[2] << "\t"; if (getAtomicForceAtStep(step).Norm() > MYEPSILON) *out << "\t" << scientific << setprecision(6) << getAtomicForceAtStep(step)[0] << "\t" << getAtomicForceAtStep(step)[1] << "\t" << getAtomicForceAtStep(step)[2] << "\t"; *out << "\t# Number in molecule " << getNr() << endl; return true; } else return false; }; bool atom::OutputTrajectoryXYZ(ofstream * const out, const int step) const { if (out != NULL) { *out << getType()->getSymbol() << "\t"; *out << getPositionAtStep(step)[0] << "\t"; *out << getPositionAtStep(step)[1] << "\t"; *out << getPositionAtStep(step)[2] << endl; return true; } else return false; }; void atom::OutputMPQCLine(ostream * const out, const Vector *center) const { Vector recentered(getPosition()); recentered -= *center; *out << "\t\t" << getType()->getSymbol() << " [ " << recentered[0] << "\t" << recentered[1] << "\t" << recentered[2] << " ]" << endl; }; void atom::OutputPsi3Line(ostream * const out, const Vector *center) const { Vector recentered(getPosition()); recentered -= *center; *out << "\t( " << getType()->getSymbol() << "\t" << recentered[0] << "\t" << recentered[1] << "\t" << recentered[2] << " )" << endl; }; bool atom::Compare(const atom &ptr) const { if (getNr() < ptr.getNr()) return true; else return false; }; double atom::DistanceSquaredToVector(const Vector &origin) const { return DistanceSquared(origin); }; double atom::DistanceToVector(const Vector &origin) const { return distance(origin); }; void atom::InitComponentNr() { if (ComponentNr != NULL) delete[](ComponentNr); const BondList& ListOfBonds = getListOfBonds(); ComponentNr = new int[ListOfBonds.size()+1]; for (int i=ListOfBonds.size()+1;i--;) ComponentNr[i] = -1; }; void atom::resetGraphNr(){ GraphNr=-1; } std::ostream & atom::operator << (std::ostream &ost) const { ParticleInfo::operator<<(ost); ost << "," << getPosition(); return ost; } std::ostream & operator << (std::ostream &ost, const atom &a) { a.ParticleInfo::operator<<(ost); ost << "," << a.getPosition(); return ost; } bool operator < (atom &a, atom &b) { return a.Compare(b); }; World *atom::getWorld(){ return world; } void atom::setWorld(World* _world){ world = _world; } bool atom::changeId(atomId_t newId){ // first we move ourselves in the world // the world lets us know if that succeeded if(world->changeAtomId(id,newId,this)){ id = newId; return true; } else{ return false; } } void atom::setId(atomId_t _id) { id=_id; } atomId_t atom::getId() const { return id; } void atom::setMolecule(molecule *_mol){ // take this atom from the old molecule removeFromMolecule(); mol = _mol; if(!mol->containsAtom(this)){ mol->insert(this); } } void atom::unsetMolecule() { // take this atom from the old molecule ASSERT(!mol->containsAtom(this), "atom::unsetMolecule() - old molecule "+toString(mol)+" still contains us!"); mol = NULL; } molecule* atom::getMolecule() const { return mol; } void atom::removeFromMolecule(){ if(mol){ if(mol->containsAtom(this)){ mol->erase(this); } mol=0; } } int atom::getNr() const{ return ParticleInfo::getNr(); } atom* NewAtom(atomId_t _id){ atom * res =new atom(); res->setId(_id); return res; } void DeleteAtom(atom* atom){ delete atom; } bool compareAtomElements(atom* atom1,atom* atom2){ return atom1->getType()->getNumber() < atom2->getType()->getNumber(); }