source: src/atom.cpp@ ccd9f5

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

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

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

/** Constructor of class atom.
 */
atom::atom()
{
  father = this;  // generally, father is itself
  previous = NULL;
  next = NULL;
  Ancestor = NULL;
  type = NULL;
  sort = NULL;
  FixedIon = 0;
  GraphNr = -1;
  ComponentNr = NULL;
  IsCyclic = false;
  SeparationVertex = false;
  LowpointNr = -1;
  AdaptiveOrder = 0;
  MaxOrder = false;
  // set LCNode::Vector to our Vector
  node = &x;
};

/** Constructor of class atom.
 */
atom::atom(atom *pointer)
{
  Name = NULL;
  previous = NULL;
  next = NULL;
  father = pointer;  // generally, father is itself
  Ancestor = NULL;
  GraphNr = -1;
  ComponentNr = NULL;
  IsCyclic = false;
  SeparationVertex = false;
  LowpointNr = -1;
  AdaptiveOrder = 0;
  MaxOrder = false;
  type = pointer->type;  // copy element of atom
  x.CopyVector(&pointer->x); // copy coordination
  v.CopyVector(&pointer->v); // copy velocity
  FixedIon = pointer->FixedIon;
  nr = -1;
  sort = &nr;
  node = &x;
}


/** Destructor of class atom.
 */
atom::~atom()
{
  Free<int>(&ComponentNr, "atom::~atom: *ComponentNr");
  Free<char>(&Name, "atom::~atom: *Name");
  Trajectory.R.clear();
  Trajectory.U.clear();
  Trajectory.F.clear();
};


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

/** Output of a single atom.
 * \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::Output(ofstream *out, int ElementNo, 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;
};
bool atom::Output(ofstream *out, int *ElementNo, int *AtomNo, const char *comment)
{
  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 *out, int *ElementNo, int *AtomNo, 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 *out, 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;
};

/** Prints all bonds of this atom from given global lists.
 * \param *out stream to output to
 * \param *NumberOfBondsPerAtom array with number of bonds per atomic index
 * \param ***ListOfBondsPerAtom array per atomic index of array with pointer to bond
 * \return true - \a *out present, false - \a *out is NULL
 */
bool atom::OutputBondOfAtom(ofstream *out, int *NumberOfBondsPerAtom, bond ***ListOfBondsPerAtom) const
{
  if (out != NULL) {
#ifdef ADDHYDROGEN
    if (type->Z != 1) {   // regard only non-hydrogen
#endif
      *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << NumberOfBondsPerAtom[nr] << " bonds: ";
      int TotalDegree = 0;
      for (int j=0;j<NumberOfBondsPerAtom[nr];j++) {
        *out << *ListOfBondsPerAtom[nr][j] << "\t";
        TotalDegree += ListOfBondsPerAtom[nr][j]->BondDegree;
      }
      *out << " -- TotalDegree: " << TotalDegree << endl;
#ifdef ADDHYDROGEN
    }
#endif
    return true;
  } else
    return false;
};

ostream & operator << (ostream &ost, const atom &a)
{
  ost << "[" << a.Name << "|" << &a << "]";
  return ost;
};

ostream & atom::operator << (ostream &ost)
{
  ost << "[" << Name << "|" << this << "]";
  return ost;
};

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

/** Adds kinetic energy of this atom to given temperature value.
 * \param *temperature add on this value
 * \param step given step of trajectory to add
 */
void atom::AddKineticToTemperature(double *temperature, int step) const
{
  for (int i=NDIM;i--;)
    *temperature += type->mass * Trajectory.U.at(step).x[i]* Trajectory.U.at(step).x[i];
};

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

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

/** Evaluates some constraint potential if atom moves from \a startstep at once to \endstep in trajectory.
 * \param startstep trajectory begins at
 * \param endstep trajectory ends at
 * \param **PermutationMap if atom switches places with some other atom, there is no translation but a permutaton noted here (not in the trajectories of each).
 * \param *Force Force matrix to store result in
 */
void atom::EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
{
  double constant = 10.;
  atom *Sprinter = PermutationMap[nr];
  // set forces
  for (int i=NDIM;i++;)
    Force->Matrix[0][nr][5+i] += 2.*constant*sqrt(Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
};