source: src/Atom/atom.cpp@ ac5c36

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 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/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 <iomanip>
#include <iostream>

/************************************* 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()
{
  const atom *father = const_cast<const atom *>(this)->GetTrueFather();
  return const_cast<atom *>(father);
}

const atom *atom::GetTrueFather() const
{
  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<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" << 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<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 << 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) && (!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();
}