source: src/atom.cpp@ b8c9d8

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

#include "atom.hpp"
#include "bond.hpp"
#include "config.hpp"
#include "element.hpp"
#include "lists.hpp"
#include "memoryallocator.hpp"
#include "parser.hpp"
#include "vector.hpp"

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


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

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


/** Destructor of class atom.
 */
atom::~atom()
{
  unlink(this);
};


/** 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()
{
  atom *walker = this;
  do {
    if (walker == walker->father) // top most father is the one that points on itself
      break;
    walker = walker->father;
  } while (walker != NULL);
  return walker;
};

/** 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::IsInParallelepiped(const Vector offset, const double *parallelepiped) const
{
  return (node->IsInParallelepiped(offset, parallelepiped));
};

/** 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.x[0] << "\t" << x.x[1] << "\t" << x.x[2];
    *out << "\t" << FixedIon;
    if (v.Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << v.x[0] << "\t" << v.x[1] << "\t" << v.x[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(ofstream * 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.x[0] << "\t" << x.x[1] << "\t" << x.x[2];
    *out << "\t" << FixedIon;
    if (v.Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << v.x[0] << "\t" << v.x[1] << "\t" << v.x[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.x[0] << "\t" << x.x[1] << "\t" << x.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).x[0] << "\t" << Trajectory.R.at(step).x[1] << "\t" << Trajectory.R.at(step).x[2];
    *out << "\t" << FixedIon;
    if (Trajectory.U.at(step).Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << Trajectory.U.at(step).x[0] << "\t" << Trajectory.U.at(step).x[1] << "\t" << Trajectory.U.at(step).x[2] << "\t";
    if (Trajectory.F.at(step).Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << Trajectory.F.at(step).x[0] << "\t" << Trajectory.F.at(step).x[1] << "\t" << Trajectory.F.at(step).x[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).x[0] << "\t";
    *out << Trajectory.R.at(step).x[1] << "\t";
    *out << Trajectory.R.at(step).x[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(ofstream * const out, const Vector *center, int *AtomNo = NULL) const
{
  *out << "\t\t" << type->symbol << " [ " << x.x[0]-center->x[0] << "\t" << x.x[1]-center->x[1] << "\t" << x.x[2]-center->x[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)
    Free(&ComponentNr);
  ComponentNr = Malloc<int>(ListOfBonds.size()+1, "atom::InitComponentNumbers: *ComponentNr");
  for (int i=ListOfBonds.size()+1;i--;)
    ComponentNr[i] = -1;
};


bool operator < (atom &a, atom &b)
{
  return a.Compare(b);
};