source: src/atom.cpp@ ce3d2b

/** \file atom.cpp
 *
 * Function implementations for the class atom.
 *
 */

#include "Helpers/MemDebug.hpp"

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

#include <iomanip>

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


/** Constructor of class atom.
 */
atom::atom() :
  father(this), sort(&nr), mol(0)
{
  node = &x;  // TesselPoint::x can only be referenced from here
};

/** Constructor of class atom.
 */
atom::atom(atom *pointer) :
    ParticleInfo(pointer),father(pointer), sort(&nr)
{
  type = pointer->type;  // copy element of atom
  x = pointer->x; // copy coordination
  v = pointer->v; // copy velocity
  FixedIon = pointer->FixedIon;
  node = &x;
  mol = 0;
};

atom *atom::clone(){
  atom *res = new atom(this);
  res->father = this;
  res->sort = &res->nr;
  res->type = type;
  res->x = this->x;
  res->v = this->v;
  res->FixedIon = FixedIon;
  res->node = &x;
  res->mol = 0;
  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;
};

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

/** 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 << x[0] << "\t" << x[1] << "\t" << x[2];
    *out << "\t" << FixedIon;
    if (v.Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << v[0] << "\t" << v[1] << "\t" << v[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 int *ElementNo, int *AtomNo, const char *comment) const
{
  AtomNo[type->Z]++;  // increment number
  if (out != NULL) {
    *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
    *out << x[0] << "\t" << x[1] << "\t" << x[2];
    *out << "\t" << FixedIon;
    if (v.Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << v[0] << "\t" << v[1] << "\t" << v[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 << type->symbol << "\t" << x[0] << "\t" << x[1] << "\t" << x[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 int *ElementNo, int *AtomNo, const int step) const
{
  AtomNo[type->Z]++;
  if (out != NULL) {
    *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\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 << type->symbol << "\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, int *AtomNo = NULL) const
{
  *out << "\t\t" << type->symbol << " [ " << x[0]-center->at(0) << "\t" << x[1]-center->at(1) << "\t" << x[2]-center->at(2) << " ]" << endl;
  if (AtomNo != NULL)
    *AtomNo++;
};

/** 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 origin.DistanceSquared(x);
};

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

/** 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;
};


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

void atom::removeFromMolecule(){
  if(mol){
    if(mol->containsAtom(this)){
      mol->erase(this);
    }
    mol=0;
  }
}


atom* NewAtom(atomId_t _id){
  atom * res =new atom();
  res->setId(_id);
  return res;
}

void DeleteAtom(atom* atom){
  delete atom;
}