source: src/Atom/atom.cpp@ d649b7

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * Copyright (C)  2013 Frederik Heber. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/** \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 "AtomObserver.hpp"
#include "Bond/bond.hpp"
#include "CodePatterns/Log.hpp"
#include "config.hpp"
#include "Element/element.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "World.hpp"
#include "WorldTime.hpp"
#include "molecule.hpp"
#include "Shapes/Shape.hpp"

#include <iomanip>
#include <iostream>

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


atom::atom() :
  father(this),
  sort(&Nr),
  mol(0)
{
  // sign on to global atom change tracker
  AtomObserver::getInstance().AtomInserted(this);
}

atom::atom(atom *pointer) :
    ParticleInfo(*pointer),
    AtomInfo(*pointer),
    father(pointer),
    sort(&Nr),
    mol(0)
{
  AtomicPosition = pointer->AtomicPosition; // copy trajectory of coordination
  AtomicVelocity = pointer->AtomicVelocity; // copy trajectory of velocity
  AtomicForce = pointer->AtomicForce;
  // sign on to global atom change tracker
  AtomObserver::getInstance().AtomInserted(this);
};

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


/** Destructor of class atom.
 */
atom::~atom()
{
  removeFromMolecule();
  // sign off from global atom change tracker
  AtomObserver::getInstance().AtomRemoved(this);
}


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::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::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)){
    OBSERVE;
    id = newId;
    NOTIFY(IndexChanged);
    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;
  }
}

bool atom::changeNr(const int newNr)
{
  if ((mol) && (mol->changeAtomNr(getNr(),newNr,this))) {
    return true;
  } else{
    return false;
  }
}

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

atom* NewAtom(atomId_t _id){
  atom * res = new atom();
  // extent trajectory to current time step
  const size_t CurrentTime = WorldTime::getTime();
  for (size_t step = res->getTrajectorySize(); step <= CurrentTime; ++step)
    res->UpdateSteps();
  res->setId(_id);
  return res;
}

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

bool compareAtomElements(atom* atom1,atom* atom2){
  return atom1->getType()->getAtomicNumber() < atom2->getType()->getAtomicNumber();
}