source: src/atom.cpp@ 081e5d

/*
 * 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 <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include "atom.hpp"
#include "bond.hpp"
#include "config.hpp"
#include "element.hpp"
#include "lists.hpp"
#include "parser.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "World.hpp"
#include "molecule.hpp"
#include "Shapes/Shape.hpp"

#include <iomanip>
#include <iostream>

/************************************* Functions for class atom *************************************/


/** Constructor of class atom.
 */
atom::atom() :
  father(this),
  sort(&nr),
  mol(0)
{};

/** Constructor of class atom.
 */
atom::atom(atom *pointer) :
    ParticleInfo(pointer),
    father(pointer),
    sort(&nr),
    mol(0)
{
  setType(pointer->getType());  // copy element of atom
  setPosition(pointer->getPosition()); // copy coordination
  AtomicVelocity = pointer->AtomicVelocity; // copy velocity
  FixedIon = pointer->FixedIon;
};

atom *atom::clone(){
  atom *res = new atom(this);
  World::getInstance().registerAtom(res);
  return res;
}


/** Destructor of class atom.
 */
atom::~atom()
{
  removeFromMolecule();
  for(BondList::iterator iter=ListOfBonds.begin(); iter!=ListOfBonds.end();){
    // deleting the bond will invalidate the iterator !!!
    bond *bond =*(iter++);
    delete(bond);
  }
};


/** Climbs up the father list until NULL, last is returned.
 * \return true father, i.e. whose father points to itself, NULL if it could not be found or has none (added hydrogen)
 */
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

};

/** Check whether father is equal to given atom.
 * \param *ptr atom to compare father to
 * \param **res return value (only set if atom::father is equal to \a *ptr)
 */
void atom::EqualsFather ( const atom *ptr, const atom **res ) const
{
  if ( ptr == father )
    *res = this;
};

bool atom::isFather(const atom *ptr){
  return ptr==father;
}

/** Checks whether atom is within the given box.
 * \param offset offset to box origin
 * \param *parallelepiped box matrix
 * \return true - is inside, false - is not
 */
bool atom::IsInShape(const Shape& shape) const
{
  return shape.isInside(getPosition());
};

/** Counts the number of bonds weighted by bond::BondDegree.
 * \param bonds times bond::BondDegree
 */
int BondedParticle::CountBonds() const
{
  int NoBonds = 0;
  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
    NoBonds += (*Runner)->BondDegree;
  return NoBonds;
};

/** Output of a single atom with given numbering.
 * \param ElementNo cardinal number of the element
 * \param AtomNo cardinal number among these atoms of the same element
 * \param *out stream to output to
 * \param *comment commentary after '#' sign
  * \return true - \a *out present, false - \a *out is NULL
 */
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" << FixedIon;
    if (AtomicVelocity.Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << AtomicVelocity[0] << "\t" << AtomicVelocity[1] << "\t" << AtomicVelocity[2] << "\t";
    if (comment != NULL)
      *out << " # " << comment << endl;
    else
      *out << " # molecule nr " << nr << endl;
    return true;
  } else
    return false;
};

/** Output of a single atom with numbering from array according to atom::type.
 * \param *ElementNo cardinal number of the element
 * \param *AtomNo cardinal number among these atoms of the same element
 * \param *out stream to output to
 * \param *comment commentary after '#' sign
  * \return true - \a *out present, false - \a *out is NULL
 */
bool atom::OutputArrayIndexed(ostream * const out,const enumeration<const element*> &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" << FixedIon;
    if (AtomicVelocity.Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << AtomicVelocity[0] << "\t" << AtomicVelocity[1] << "\t" << AtomicVelocity[2] << "\t";
    if (comment != NULL)
      *out << " # " << comment << endl;
    else
      *out << " # molecule nr " << nr << endl;
    return true;
  } else
    return false;
};

/** Output of a single atom as one lin in xyz file.
 * \param *out stream to output to
  * \return true - \a *out present, false - \a *out is NULL
 */
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;
};

/** Output of a single atom as one lin in xyz file.
 * \param *out stream to output to
 * \param *ElementNo array with ion type number in the config file this atom's element shall have
 * \param *AtomNo array with atom number in the config file this atom shall have, is increase by one automatically
 * \param step Trajectory time step to output
  * \return true - \a *out present, false - \a *out is NULL
 */
bool atom::OutputTrajectory(ofstream * const out, const enumeration<const element*> &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 << Trajectory.R.at(step)[0] << "\t" << Trajectory.R.at(step)[1] << "\t" << Trajectory.R.at(step)[2];
    *out << "\t" << FixedIon;
    if (Trajectory.U.at(step).Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << Trajectory.U.at(step)[0] << "\t" << Trajectory.U.at(step)[1] << "\t" << Trajectory.U.at(step)[2] << "\t";
    if (Trajectory.F.at(step).Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << Trajectory.F.at(step)[0] << "\t" << Trajectory.F.at(step)[1] << "\t" << Trajectory.F.at(step)[2] << "\t";
    *out << "\t# Number in molecule " << nr << endl;
    return true;
  } else
    return false;
};

/** Output of a single atom as one lin in xyz file.
 * \param *out stream to output to
 * \param step Trajectory time step to output
 * \return true - \a *out present, false - \a *out is NULL
 */
bool atom::OutputTrajectoryXYZ(ofstream * const out, const int step) const
{
  if (out != NULL) {
    *out << getType()->getSymbol() << "\t";
    *out << Trajectory.R.at(step)[0] << "\t";
    *out << Trajectory.R.at(step)[1] << "\t";
    *out << Trajectory.R.at(step)[2] << endl;
    return true;
  } else
    return false;
};

/** Outputs the MPQC configuration line for this atom.
 * \param *out output stream
 * \param *center center of molecule subtracted from position
 * \param *AtomNo pointer to atom counter that is increased by one
 */
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;
};

/** Compares the indices of \a this atom with a given \a ptr.
 * \param ptr atom to compare index against
 * \return true - this one's is smaller, false - not
 */
bool atom::Compare(const atom &ptr) const
{
  if (nr < ptr.nr)
    return true;
  else
    return false;
};

/** Returns squared distance to a given vector.
 * \param origin vector to calculate distance to
 * \return distance squared
 */
double atom::DistanceSquaredToVector(const Vector &origin) const
{
  return DistanceSquared(origin);
};

/** Returns distance to a given vector.
 * \param origin vector to calculate distance to
 * \return distance
 */
double atom::DistanceToVector(const Vector &origin) const
{
  return distance(origin);
};

/** Initialises the component number array.
 * Size is set to atom::ListOfBonds.size()+1 (last is th encode end by -1)
 */
void atom::InitComponentNr()
{
  if (ComponentNr != NULL)
    delete[](ComponentNr);
  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;
}

/** Makes the atom be contained in the new molecule \a *_mol.
 * Uses atom::removeFromMolecule() to delist from old molecule.
 * \param *_mol pointer to new molecule
 */
void atom::setMolecule(molecule *_mol){
  // take this atom from the old molecule
  removeFromMolecule();
  mol = _mol;
  if(!mol->containsAtom(this)){
    mol->insert(this);
  }
}

/** Returns pointer to the molecule which atom belongs to.
 * \return containing molecule
 */
molecule* atom::getMolecule() const {
  return mol;
}

/** Erases the atom in atom::mol's list of atoms and sets it to zero.
 */
void atom::removeFromMolecule(){
  if(mol){
    if(mol->containsAtom(this)){
      mol->erase(this);
    }
    mol=0;
  }
}

double atom::getMass() const{
  return getType()->getMass();
}

int atom::getNr() const{
  return nr;
}

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();
}