source: src/molecule.cpp@ a275b3

/*
 * 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 molecules.cpp
 *
 * Functions for the class molecule.
 *
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include <cstring>
#include <boost/bind.hpp>
#include <boost/foreach.hpp>

#include <gsl/gsl_inline.h>
#include <gsl/gsl_heapsort.h>

#include "molecule.hpp"

#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "Box.hpp"
#include "CodePatterns/enumeration.hpp"
#include "CodePatterns/Log.hpp"
#include "config.hpp"
#include "Descriptors/AtomIdDescriptor.hpp"
#include "Element/element.hpp"
#include "Graph/BondGraph.hpp"
#include "LinearAlgebra/Exceptions.hpp"
#include "LinearAlgebra/leastsquaremin.hpp"
#include "LinearAlgebra/Plane.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinkedCell/linkedcell.hpp"
#include "IdPool_impl.hpp"
#include "Shapes/BaseShapes.hpp"
#include "Tesselation/tesselation.hpp"
#include "World.hpp"
#include "WorldTime.hpp"


/************************************* Functions for class molecule *********************************/

/** Constructor of class molecule.
 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
 */
molecule::molecule() :
  Observable("molecule"),
  MDSteps(0),
  NoNonHydrogen(0),
  NoNonBonds(0),
  NoCyclicBonds(0),
  ActiveFlag(false),
  IndexNr(-1),
  AtomCount(this,boost::bind(&molecule::doCountAtoms,this),"AtomCount"),
  BondCount(this,boost::bind(&molecule::doCountBonds,this),"BondCount"),
  atomIdPool(1, 20, 100),
  last_atom(0)
{

  strcpy(name,World::getInstance().getDefaultName().c_str());
};

molecule *NewMolecule(){
  return new molecule();
}

/** Destructor of class molecule.
 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
 */
molecule::~molecule()
{
  CleanupMolecule();
};


void DeleteMolecule(molecule *mol){
  delete mol;
}

// getter and setter
const std::string molecule::getName() const{
  return std::string(name);
}

int molecule::getAtomCount() const{
  return *AtomCount;
}

int molecule::getBondCount() const{
  return *BondCount;
}

void molecule::setName(const std::string _name){
  OBSERVE;
  cout << "Set name of molecule " << getId() << " to " << _name << endl;
  strncpy(name,_name.c_str(),MAXSTRINGSIZE);
}

bool molecule::changeAtomNr(int oldNr, int newNr, atom* target){
  OBSERVE;
  if(atomIdPool.reserveId(newNr)){
    if (oldNr != -1)  // -1 is reserved and indicates no number
      atomIdPool.releaseId(oldNr);
    ASSERT (target,
        "molecule::changeAtomNr() - given target is NULL, cannot set Nr or name.");
    target->setNr(newNr);
    setAtomName(target);
    return true;
  } else{
    return false;
  }
}

bool molecule::changeId(moleculeId_t newId){
  // first we move ourselves in the world
  // the world lets us know if that succeeded
  if(World::getInstance().changeMoleculeId(id,newId,this)){
    id = newId;
    return true;
  }
  else{
    return false;
  }
}


moleculeId_t molecule::getId() const {
  return id;
}

void molecule::setId(moleculeId_t _id){
  id =_id;
}

const Formula &molecule::getFormula() const {
  return formula;
}

unsigned int molecule::getElementCount() const{
  return formula.getElementCount();
}

bool molecule::hasElement(const element *element) const{
  return formula.hasElement(element);
}

bool molecule::hasElement(atomicNumber_t Z) const{
  return formula.hasElement(Z);
}

bool molecule::hasElement(const string &shorthand) const{
  return formula.hasElement(shorthand);
}

/************************** Access to the List of Atoms ****************/

molecule::const_iterator molecule::erase( const_iterator loc )
{
  OBSERVE;
  const_iterator iter = loc;
  ++iter;
  atom * const _atom = const_cast<atom *>(*loc);
  atomIds.erase( _atom->getId() );
  formula-=_atom->getType();
  _atom->removeFromMolecule();
  return iter;
}

molecule::const_iterator molecule::erase( atom * key )
{
  OBSERVE;
  const_iterator iter = find(key);
  if (iter != end()){
    ++iter;
    atomIds.erase( key->getId() );
    atomIdPool.releaseId(key->getNr());
    key->setNr(-1);
    formula-=key->getType();
    key->removeFromMolecule();
  }
  return iter;
}

pair<molecule::iterator,bool> molecule::insert ( atom * const key )
{
  OBSERVE;
  std::pair<iterator,bool> res = atomIds.insert(key->getId());
  if (res.second) { // push atom if went well
    key->setNr(atomIdPool.getNextId());
    setAtomName(key);
    formula+=key->getType();
    return res;
  } else {
    return pair<iterator,bool>(end(),res.second);
  }
}

void molecule::setAtomName(atom *_atom) const
{
  std::stringstream sstr;
  sstr << _atom->getType()->getSymbol() << _atom->getNr();
  _atom->setName(sstr.str());
}


bool molecule::containsAtom(atom* key){
  return (find(key) != end());
}

World::AtomComposite molecule::getAtomSet() const
{
  World::AtomComposite vector_of_atoms;
  for (molecule::iterator iter = begin(); iter != end(); ++iter)
    vector_of_atoms.push_back(*iter);
  return vector_of_atoms;
}

/** Adds given atom \a *pointer from molecule list.
 * Increases molecule::last_atom and gives last number to added atom and names it according to its element::abbrev and molecule::AtomCount
 * \param *pointer allocated and set atom
 * \return true - succeeded, false - atom not found in list
 */
bool molecule::AddAtom(atom *pointer)
{
  OBSERVE;
  if (pointer != NULL) {
    if (pointer->getType()->getAtomicNumber() != 1)
      NoNonHydrogen++;
    insert(pointer);
    pointer->setMolecule(this);
  }
  return true;
};

/** Adds a copy of the given atom \a *pointer from molecule list.
 * Increases molecule::last_atom and gives last number to added atom.
 * \param *pointer allocated and set atom
 * \return pointer to the newly added atom
 */
atom * molecule::AddCopyAtom(atom *pointer)
{
  atom *retval = NULL;
  OBSERVE;
  if (pointer != NULL) {
    atom *walker = pointer->clone();
    walker->setName(pointer->getName());
    walker->setNr(last_atom++);  // increase number within molecule
    insert(walker);
    if ((pointer->getType() != NULL) && (pointer->getType()->getAtomicNumber() != 1))
      NoNonHydrogen++;
    walker->setMolecule(this);
    retval=walker;
  }
  return retval;
};

/** Adds a Hydrogen atom in replacement for the given atom \a *partner in bond with a *origin.
 * Here, we have to distinguish between single, double or triple bonds as stated by \a BondDegree, that each demand
 * a different scheme when adding \a *replacement atom for the given one.
 * -# Single Bond: Simply add new atom with bond distance rescaled to typical hydrogen one
 * -# Double Bond: Here, we need the **BondList of the \a *origin atom, by scanning for the other bonds instead of
 *    *Bond, we use the through these connected atoms to determine the plane they lie in, vector::MakeNormalvector().
 *    The orthonormal vector to this plane along with the vector in *Bond direction determines the plane the two
 *    replacing hydrogens shall lie in. Now, all remains to do is take the usual hydrogen double bond angle for the
 *    element of *origin and form the sin/cos admixture of both plane vectors for the new coordinates of the two
 *    hydrogens forming this angle with *origin.
 * -# Triple Bond: The idea is to set up a tetraoid (C1-H1-H2-H3) (however the lengths \f$b\f$ of the sides of the base
 *    triangle formed by the to be added hydrogens are not equal to the typical bond distance \f$l\f$ but have to be
 *    determined from the typical angle \f$\alpha\f$ for a hydrogen triple connected to the element of *origin):
 *    We have the height \f$d\f$ as the vector in *Bond direction (from triangle C1-H1-H2).
 *    \f[ h = l \cdot \cos{\left (\frac{\alpha}{2} \right )} \qquad b = 2l \cdot \sin{\left (\frac{\alpha}{2} \right)} \quad \rightarrow \quad d = l \cdot \sqrt{\cos^2{\left (\frac{\alpha}{2} \right)}-\frac{1}{3}\cdot\sin^2{\left (\frac{\alpha}{2}\right )}}
 *    \f]
 *    vector::GetNormalvector() creates one orthonormal vector from this *Bond vector and vector::MakeNormalvector creates
 *    the third one from the former two vectors. The latter ones form the plane of the base triangle mentioned above.
 *    The lengths for these are \f$f\f$ and \f$g\f$ (from triangle H1-H2-(center of H1-H2-H3)) with knowledge that
 *    the median lines in an isosceles triangle meet in the center point with a ratio 2:1.
 *    \f[ f = \frac{b}{\sqrt{3}} \qquad g = \frac{b}{2}
 *    \f]
 *    as the coordination of all three atoms in the coordinate system of these three vectors:
 *    \f$\pmatrix{d & f & 0}\f$, \f$\pmatrix{d & -0.5 \cdot f & g}\f$ and \f$\pmatrix{d & -0.5 \cdot f & -g}\f$.
 *
 * \param *out output stream for debugging
 * \param *Bond pointer to bond between \a *origin and \a *replacement
 * \param *TopOrigin son of \a *origin of upper level molecule (the atom added to this molecule as a copy of \a *origin)
 * \param *origin pointer to atom which acts as the origin for scaling the added hydrogen to correct bond length
 * \param *replacement pointer to the atom which shall be copied as a hydrogen atom in this molecule
 * \param isAngstroem whether the coordination of the given atoms is in AtomicLength (false) or Angstrom(true)
 * \return number of atoms added, if < bond::BondDegree then something went wrong
 * \todo double and triple bonds splitting (always use the tetraeder angle!)
 */
bool molecule::AddHydrogenReplacementAtom(bond *TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem)
{
//  Info info(__func__);
  bool AllWentWell = true;    // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
  OBSERVE;
  double bondlength;  // bond length of the bond to be replaced/cut
  double bondangle;  // bond angle of the bond to be replaced/cut
  double BondRescale;   // rescale value for the hydrogen bond length
  bond *FirstBond = NULL, *SecondBond = NULL; // Other bonds in double bond case to determine "other" plane
  atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
  double b,l,d,f,g, alpha, factors[NDIM];    // hold temporary values in triple bond case for coordination determination
  Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
  Vector InBondvector;    // vector in direction of *Bond
  const RealSpaceMatrix &matrix =  World::getInstance().getDomain().getM();
  bond *Binder = NULL;

  // create vector in direction of bond
  InBondvector = TopReplacement->getPosition() - TopOrigin->getPosition();
  bondlength = InBondvector.Norm();

   // is greater than typical bond distance? Then we have to correct periodically
   // the problem is not the H being out of the box, but InBondvector have the wrong direction
   // due to TopReplacement or Origin being on the wrong side!
  const BondGraph * const BG = World::getInstance().getBondGraph();
  const range<double> MinMaxBondDistance(
      BG->getMinMaxDistance(TopOrigin,TopReplacement));
  if (!MinMaxBondDistance.isInRange(bondlength)) {
//    LOG(4, "InBondvector is: " << InBondvector << ".");
    Orthovector1.Zero();
    for (int i=NDIM;i--;) {
      l = TopReplacement->at(i) - TopOrigin->at(i);
      if (fabs(l) > MinMaxBondDistance.last) { // is component greater than bond distance (check against min not useful here)
        Orthovector1[i] = (l < 0) ? -1. : +1.;
      } // (signs are correct, was tested!)
    }
    Orthovector1 *= matrix;
    InBondvector -= Orthovector1; // subtract just the additional translation
    bondlength = InBondvector.Norm();
//    LOG(4, "INFO: Corrected InBondvector is now: " << InBondvector << ".");
  } // periodic correction finished

  InBondvector.Normalize();
  // get typical bond length and store as scale factor for later
  ASSERT(TopOrigin->getType() != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
  BondRescale = TopOrigin->getType()->getHBondDistance(TopBond->BondDegree-1);
  if (BondRescale == -1) {
    ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!");
    return false;
    BondRescale = bondlength;
  } else {
    if (!IsAngstroem)
      BondRescale /= (1.*AtomicLengthToAngstroem);
  }

  // discern single, double and triple bonds
  switch(TopBond->BondDegree) {
    case 1:
      FirstOtherAtom = World::getInstance().createAtom();    // new atom
      FirstOtherAtom->setType(1);  // element is Hydrogen
      FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
      if (TopReplacement->getType()->getAtomicNumber() == 1) { // neither rescale nor replace if it's already hydrogen
        FirstOtherAtom->father = TopReplacement;
        BondRescale = bondlength;
      } else {
        FirstOtherAtom->father = NULL;  // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
      }
      InBondvector *= BondRescale;   // rescale the distance vector to Hydrogen bond length
      FirstOtherAtom->setPosition(TopOrigin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
//      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = GraphEdge::TreeEdge;
      break;
    case 2:
      {
        // determine two other bonds (warning if there are more than two other) plus valence sanity check
        const BondList& ListOfBonds = TopOrigin->getListOfBonds();
        for (BondList::const_iterator Runner = ListOfBonds.begin();
            Runner != ListOfBonds.end();
            ++Runner) {
          if ((*Runner) != TopBond) {
            if (FirstBond == NULL) {
              FirstBond = (*Runner);
              FirstOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
            } else if (SecondBond == NULL) {
              SecondBond = (*Runner);
              SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
            } else {
              ELOG(2, "Detected more than four bonds for atom " << TopOrigin->getName());
            }
          }
        }
      }
      if (SecondOtherAtom == NULL) {  // then we have an atom with valence four, but only 3 bonds: one to replace and one which is TopBond (third is FirstBond)
        SecondBond = TopBond;
        SecondOtherAtom = TopReplacement;
      }
      if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
//        LOG(3, "Regarding the double bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") to be constructed: Taking " << FirstOtherAtom->Name << " and " << SecondOtherAtom->Name << " along with " << TopOrigin->Name << " to determine orthogonal plane.");

        // determine the plane of these two with the *origin
        try {
          Orthovector1 = Plane(TopOrigin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
        }
        catch(LinearDependenceException &excp){
          LOG(0, boost::diagnostic_information(excp));
          // TODO: figure out what to do with the Orthovector in this case
          AllWentWell = false;
        }
      } else {
        Orthovector1.GetOneNormalVector(InBondvector);
      }
      //LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
      // orthogonal vector and bond vector between origin and replacement form the new plane
      Orthovector1.MakeNormalTo(InBondvector);
      Orthovector1.Normalize();
      //LOG(3, "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << ".");

      // create the two Hydrogens ...
      FirstOtherAtom = World::getInstance().createAtom();
      SecondOtherAtom = World::getInstance().createAtom();
      FirstOtherAtom->setType(1);
      SecondOtherAtom->setType(1);
      FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
      SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
      SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
      FirstOtherAtom->father = NULL;  // we are just an added hydrogen with no father
      SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
      bondangle = TopOrigin->getType()->getHBondAngle(1);
      if (bondangle == -1) {
        ELOG(1, "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!");
        return false;
        bondangle = 0;
      }
      bondangle *= M_PI/180./2.;
//      LOG(3, "INFO: ReScaleCheck: InBondvector " << InBondvector << ", " << Orthovector1 << ".");
//      LOG(3, "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle));
      FirstOtherAtom->Zero();
      SecondOtherAtom->Zero();
      for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
        FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
        SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
      }
      FirstOtherAtom->Scale(BondRescale);  // rescale by correct BondDistance
      SecondOtherAtom->Scale(BondRescale);
      //LOG(3, "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << ".");
      *FirstOtherAtom += TopOrigin->getPosition();
      *SecondOtherAtom += TopOrigin->getPosition();
      // ... and add to molecule
      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
      AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
//      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
//      LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = GraphEdge::TreeEdge;
      Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = GraphEdge::TreeEdge;
      break;
    case 3:
      // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
      FirstOtherAtom = World::getInstance().createAtom();
      SecondOtherAtom = World::getInstance().createAtom();
      ThirdOtherAtom = World::getInstance().createAtom();
      FirstOtherAtom->setType(1);
      SecondOtherAtom->setType(1);
      ThirdOtherAtom->setType(1);
      FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
      SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
      SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
      ThirdOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
      ThirdOtherAtom->setFixedIon(TopReplacement->getFixedIon());
      FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
      SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
      ThirdOtherAtom->father = NULL;  //  we are just an added hydrogen with no father

      // we need to vectors orthonormal the InBondvector
      AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(InBondvector);
//      LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
      try{
        Orthovector2 = Plane(InBondvector, Orthovector1,0).getNormal();
      }
      catch(LinearDependenceException &excp) {
        LOG(0, boost::diagnostic_information(excp));
        AllWentWell = false;
      }
//      LOG(3, "INFO: Orthovector2: " << Orthovector2 << ".")

      // create correct coordination for the three atoms
      alpha = (TopOrigin->getType()->getHBondAngle(2))/180.*M_PI/2.;  // retrieve triple bond angle from database
      l = BondRescale;        // desired bond length
      b = 2.*l*sin(alpha);    // base length of isosceles triangle
      d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.);   // length for InBondvector
      f = b/sqrt(3.);   // length for Orthvector1
      g = b/2.;         // length for Orthvector2
//      LOG(3, "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", ");
//      LOG(3, "The three Bond lengths: " << sqrt(d*d+f*f) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << ", "  << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g));
      factors[0] = d;
      factors[1] = f;
      factors[2] = 0.;
      FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
      factors[1] = -0.5*f;
      factors[2] = g;
      SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
      factors[2] = -g;
      ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);

      // rescale each to correct BondDistance
//      FirstOtherAtom->x.Scale(&BondRescale);
//      SecondOtherAtom->x.Scale(&BondRescale);
//      ThirdOtherAtom->x.Scale(&BondRescale);

      // and relative to *origin atom
      *FirstOtherAtom += TopOrigin->getPosition();
      *SecondOtherAtom += TopOrigin->getPosition();
      *ThirdOtherAtom += TopOrigin->getPosition();

      // ... and add to molecule
      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
      AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
      AllWentWell = AllWentWell && AddAtom(ThirdOtherAtom);
//      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
//      LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
//      LOG(4, "INFO: Added " << *ThirdOtherAtom << " at: " << ThirdOtherAtom->x << ".");
      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = GraphEdge::TreeEdge;
      Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = GraphEdge::TreeEdge;
      Binder = AddBond(BottomOrigin, ThirdOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = GraphEdge::TreeEdge;
      break;
    default:
      ELOG(1, "BondDegree does not state single, double or triple bond!");
      AllWentWell = false;
      break;
  }

  return AllWentWell;
};

/** Creates a copy of this molecule.
 * \param offset translation Vector for the new molecule relative to old one
 * \return copy of molecule
 */
molecule *molecule::CopyMolecule(const Vector &offset) const
{
  molecule *copy = World::getInstance().createMolecule();

  // copy all atoms
  std::map< const atom *, atom *> FatherFinder;
  for (iterator iter = begin(); iter != end(); ++iter) {
    atom * const copy_atom = copy->AddCopyAtom(*iter);
    copy_atom->setPosition(copy_atom->getPosition() + offset);
    FatherFinder.insert( std::make_pair( *iter, copy_atom ) );
  }

  // copy all bonds
  for(const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
        BondRunner != ListOfBonds.end();
        ++BondRunner)
      if ((*BondRunner)->leftatom == *AtomRunner) {
        bond *Binder = (*BondRunner);
        // get the pendant atoms of current bond in the copy molecule
        ASSERT(FatherFinder.count(Binder->leftatom),
            "molecule::CopyMolecule() - No copy of original left atom "
            +toString(Binder->leftatom)+" for bond copy found");
        ASSERT(FatherFinder.count(Binder->rightatom),
            "molecule::CopyMolecule() - No copy of original right atom "
            +toString(Binder->rightatom)+"  for bond copy found");
        atom * const LeftAtom = FatherFinder[Binder->leftatom];
        atom * const RightAtom = FatherFinder[Binder->rightatom];

        bond * const NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
        NewBond->Cyclic = Binder->Cyclic;
        if (Binder->Cyclic)
          copy->NoCyclicBonds++;
        NewBond->Type = Binder->Type;
      }
  }
  // correct fathers
  //for_each(begin(),end(),mem_fun(&atom::CorrectFather));

  return copy;
};


/** Destroys all atoms inside this molecule.
 */
void molecule::removeAtomsinMolecule()
{
  // remove each atom from world
  for(iterator AtomRunner = begin(); !empty(); AtomRunner = begin())
    World::getInstance().destroyAtom(*AtomRunner);
};


/**
 * Copies all atoms of a molecule which are within the defined parallelepiped.
 *
 * @param offest for the origin of the parallelepiped
 * @param three vectors forming the matrix that defines the shape of the parallelpiped
 */
molecule* molecule::CopyMoleculeFromSubRegion(const Shape &region) const {
  molecule *copy = World::getInstance().createMolecule();

  // copy all atoms
  std::map< const atom *, atom *> FatherFinder;
  for (iterator iter = begin(); iter != end(); ++iter) {
    if((*iter)->IsInShape(region)){
      atom * const copy_atom = copy->AddCopyAtom(*iter);
      FatherFinder.insert( std::make_pair( *iter, copy_atom ) );
    }
  }

  // copy all bonds
  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
        BondRunner != ListOfBonds.end();
        ++BondRunner)
      if ((*BondRunner)->leftatom == *AtomRunner) {
        bond *Binder = (*BondRunner);
        if ((FatherFinder.count(Binder->leftatom))
            && (FatherFinder.count(Binder->rightatom))) {
          // if copy present, then it must be from subregion
          atom * const LeftAtom = FatherFinder[Binder->leftatom];
          atom * const RightAtom = FatherFinder[Binder->rightatom];

          bond * const NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
          NewBond->Cyclic = Binder->Cyclic;
          if (Binder->Cyclic)
            copy->NoCyclicBonds++;
          NewBond->Type = Binder->Type;
        }
      }
  }
  // correct fathers
  //for_each(begin(),end(),mem_fun(&atom::CorrectFather));

  //TODO: copy->BuildInducedSubgraph(this);

  return copy;
}

/** Adds a bond to a the molecule specified by two atoms, \a *first and \a *second.
 * Also updates molecule::BondCount and molecule::NoNonBonds.
 * \param *first first atom in bond
 * \param *second atom in bond
 * \return pointer to bond or NULL on failure
 */
bond * molecule::AddBond(atom *atom1, atom *atom2, int degree)
{
  OBSERVE;
  bond *Binder = NULL;

  // some checks to make sure we are able to create the bond
  ASSERT(atom1,
      "molecule::AddBond() - First atom "+toString(atom1)
      +" is not a invalid pointer");
  ASSERT(atom2,
      "molecule::AddBond() - Second atom "+toString(atom2)
      +" is not a invalid pointer");
  ASSERT(isInMolecule(atom1),
      "molecule::AddBond() - First atom "+toString(atom1)
      +" is not part of molecule");
  ASSERT(isInMolecule(atom2),
      "molecule::AddBond() - Second atom "+toString(atom2)
      +" is not part of molecule");

  Binder = new bond(atom1, atom2, degree);
  atom1->RegisterBond(WorldTime::getTime(), Binder);
  atom2->RegisterBond(WorldTime::getTime(), Binder);
  if ((atom1->getType() != NULL)
      && (atom1->getType()->getAtomicNumber() != 1)
      && (atom2->getType() != NULL)
      && (atom2->getType()->getAtomicNumber() != 1))
    NoNonBonds++;

  return Binder;
};

/** Remove bond from bond chain list and from the both atom::ListOfBonds.
 * Bond::~Bond takes care of bond removal
 * \param *pointer bond pointer
 * \return true - bound found and removed, false - bond not found/removed
 */
bool molecule::RemoveBond(bond *pointer)
{
  //ELOG(1, "molecule::RemoveBond: Function not implemented yet.");
  delete(pointer);
  return true;
};

/** Remove every bond from bond chain list that atom \a *BondPartner is a constituent of.
 * \todo Function not implemented yet
 * \param *BondPartner atom to be removed
 * \return true - bounds found and removed, false - bonds not found/removed
 */
bool molecule::RemoveBonds(atom *BondPartner)
{
  //ELOG(1, "molecule::RemoveBond: Function not implemented yet.");
  BondPartner->removeAllBonds();
  return false;
};

/** Set molecule::name from the basename without suffix in the given \a *filename.
 * \param *filename filename
 */
void molecule::SetNameFromFilename(const char *filename)
{
  int length = 0;
  const char *molname = strrchr(filename, '/');
  if (molname != NULL)
    molname += sizeof(char);  // search for filename without dirs
  else
    molname = filename; // contains no slashes
  const char *endname = strchr(molname, '.');
  if ((endname == NULL) || (endname < molname))
    length = strlen(molname);
  else
    length = strlen(molname) - strlen(endname);
  cout << "Set name of molecule " << getId() << " to " << molname << endl;
  strncpy(name, molname, length);
  name[length]='\0';
};

/** Sets the molecule::cell_size to the components of \a *dim (rectangular box)
 * \param *dim vector class
 */
void molecule::SetBoxDimension(Vector *dim)
{
  RealSpaceMatrix domain;
  for(int i =0; i<NDIM;++i)
    domain.at(i,i) = dim->at(i);
  World::getInstance().setDomain(domain);
};

/** Removes atom from molecule list and removes all of its bonds.
 * \param *pointer atom to be removed
 * \return true - succeeded, false - atom not found in list
 */
bool molecule::RemoveAtom(atom *pointer)
{
  ASSERT(pointer, "Null pointer passed to molecule::RemoveAtom().");
  OBSERVE;
  RemoveBonds(pointer);
  pointer->removeFromMolecule();
  return true;
};

/** Removes atom from molecule list, but does not delete it.
 * \param *pointer atom to be removed
 * \return true - succeeded, false - atom not found in list
 */
bool molecule::UnlinkAtom(atom *pointer)
{
  if (pointer == NULL)
    return false;
  pointer->removeFromMolecule();
  return true;
};

/** Removes every atom from molecule list.
 * \return true - succeeded, false - atom not found in list
 */
bool molecule::CleanupMolecule()
{
  for (molecule::iterator iter = begin(); !empty(); iter = begin())
    (*iter)->removeFromMolecule();
  return empty();
};

/** Finds an atom specified by its continuous number.
 * \param Nr number of atom withim molecule
 * \return pointer to atom or NULL
 */
atom * molecule::FindAtom(int Nr)  const
{
  molecule::iterator iter = begin();
  for (; iter != end(); ++iter)
  if ((*iter)->getNr() == Nr)
    break;
  if (iter != end()) {
    //LOG(0, "Found Atom Nr. " << walker->getNr());
    return (*iter);
  } else {
    ELOG(1, "Atom not found in molecule " << getName() << "'s list.");
    return NULL;
  }
}

/** Checks whether the given atom is a member of this molecule.
 *
 *  We make use here of molecule::atomIds to get a result on
 *
 * @param _atom atom to check
 * @return true - is member, false - is not
 */
bool molecule::isInMolecule(const atom * const _atom)
{
  ASSERT(_atom->getMolecule() == this,
      "molecule::isInMolecule() - atom is not designated to be in molecule '"
      +toString(this->getName())+"'.");
  molecule::const_iterator iter = atomIds.find(_atom->getId());
  return (iter != atomIds.end());
}

/** Asks for atom number, and checks whether in list.
 * \param *text question before entering
 */
atom * molecule::AskAtom(std::string text)
{
  int No;
  atom *ion = NULL;
  do {
    //std::cout << "============Atom list==========================" << std::endl;
    //mol->Output((ofstream *)&cout);
    //std::cout << "===============================================" << std::endl;
    std::cout << text;
    cin >> No;
    ion = this->FindAtom(No);
  } while (ion == NULL);
  return ion;
};

/** Checks if given coordinates are within cell volume.
 * \param *x array of coordinates
 * \return true - is within, false - out of cell
 */
bool molecule::CheckBounds(const Vector *x) const
{
  const RealSpaceMatrix &domain = World::getInstance().getDomain().getM();
  bool result = true;
  for (int i=0;i<NDIM;i++) {
    result = result && ((x->at(i) >= 0) && (x->at(i) < domain.at(i,i)));
  }
  //return result;
  return true; /// probably not gonna use the check no more
};

/** Prints molecule to *out.
 * \param *out output stream
 */
bool molecule::Output(ostream * const output) const
{
  if (output == NULL) {
    return false;
  } else {
    int AtomNo[MAX_ELEMENTS];
    memset(AtomNo,0,(MAX_ELEMENTS-1)*sizeof(*AtomNo));
    enumeration<const element*> elementLookup = formula.enumerateElements();
    *output << "#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)" << endl;
    for_each(begin(),end(),boost::bind(&atom::OutputArrayIndexed,_1,output,elementLookup,AtomNo,(const char*)0));
    return true;
  }
};

/** Prints molecule with all atomic trajectory positions to *out.
 * \param *out output stream
 */
bool molecule::OutputTrajectories(ofstream * const output) const
{
  if (output == NULL) {
    return false;
  } else {
    for (int step = 0; step < MDSteps; step++) {
      if (step == 0) {
        *output << "#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)" << endl;
      } else {
        *output << "# ====== MD step " << step << " =========" << endl;
      }
      int AtomNo[MAX_ELEMENTS];
      memset(AtomNo,0,(MAX_ELEMENTS-1)*sizeof(*AtomNo));
      enumeration<const element*> elementLookup = formula.enumerateElements();
      for_each(begin(),end(),boost::bind(&atom::OutputTrajectory,_1,output,elementLookup, AtomNo, (const int)step));
    }
    return true;
  }
};

/** Outputs contents of each atom::ListOfBonds.
 * \param *out output stream
 */
void molecule::OutputListOfBonds() const
{
  std::stringstream output;
  LOG(2, "From Contents of ListOfBonds, all atoms:");
  for (molecule::const_iterator iter = begin();
      iter != end();
      ++iter) {
    (*iter)->OutputBondOfAtom(output);
    output << std::endl << "\t\t";
  }
  LOG(2, output.str());
}

/** Output of element before the actual coordination list.
 * \param *out stream pointer
 */
bool molecule::Checkout(ofstream * const output)  const
{
  return formula.checkOut(output);
};

/** Prints molecule with all its trajectories to *out as xyz file.
 * \param *out output stream
 */
bool molecule::OutputTrajectoriesXYZ(ofstream * const output)
{
  time_t now;

  if (output != NULL) {
    now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
    for (int step=0;step<MDSteps;step++) {
      *output << getAtomCount() << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
      for_each(begin(),end(),boost::bind(&atom::OutputTrajectoryXYZ,_1,output,step));
    }
    return true;
  } else
    return false;
};

/** Prints molecule to *out as xyz file.
* \param *out output stream
 */
bool molecule::OutputXYZ(ofstream * const output) const
{
  time_t now;

  if (output != NULL) {
    now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
    *output << getAtomCount() << "\n\tCreated by molecuilder on " << ctime(&now);
    for_each(begin(),end(),bind2nd(mem_fun(&atom::OutputXYZLine),output));
    return true;
  } else
    return false;
};

/** Brings molecule::AtomCount and atom::*Name up-to-date.
 * \param *out output stream for debugging
 */
int molecule::doCountAtoms()
{
  int res = size();
  NoNonHydrogen = 0;
  // go through atoms and look for new ones
  for (molecule::const_iterator iter = begin(); iter != end(); ++iter)
    if ((*iter)->getType()->getAtomicNumber() != 1) // count non-hydrogen atoms whilst at it
      NoNonHydrogen++;
  return res;
};

/** Counts the number of present bonds.
 * \return number of bonds
 */
int molecule::doCountBonds() const
{
  unsigned int counter = 0;
  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
        BondRunner != ListOfBonds.end();
        ++BondRunner)
      if ((*BondRunner)->leftatom == *AtomRunner)
        counter++;
  }
  return counter;
}


/** Returns an index map for two father-son-molecules.
 * The map tells which atom in this molecule corresponds to which one in the other molecul with their fathers.
 * \param *out output stream for debugging
 * \param *OtherMolecule corresponding molecule with fathers
 * \return allocated map of size molecule::AtomCount with map
 * \todo make this with a good sort O(n), not O(n^2)
 */
int * molecule::GetFatherSonAtomicMap(molecule *OtherMolecule)
{
  LOG(3, "Begin of GetFatherAtomicMap.");
  int *AtomicMap = new int[getAtomCount()];
  for (int i=getAtomCount();i--;)
    AtomicMap[i] = -1;
  if (OtherMolecule == this) {  // same molecule
    for (int i=getAtomCount();i--;) // no need as -1 means already that there is trivial correspondence
      AtomicMap[i] = i;
    LOG(4, "Map is trivial.");
  } else {
    std::stringstream output;
    output << "Map is ";
    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
      if ((*iter)->father == NULL) {
        AtomicMap[(*iter)->getNr()] = -2;
      } else {
        for (molecule::const_iterator runner = OtherMolecule->begin(); runner != OtherMolecule->end(); ++runner) {
      //for (int i=0;i<AtomCount;i++) { // search atom
        //for (int j=0;j<OtherMolecule->getAtomCount();j++) {
          //LOG(4, "Comparing father " << (*iter)->father << " with the other one " << (*runner)->father << ".");
          if ((*iter)->father == (*runner))
            AtomicMap[(*iter)->getNr()] = (*runner)->getNr();
        }
      }
      output << AtomicMap[(*iter)->getNr()] << "\t";
    }
    LOG(4, output.str());
  }
  LOG(3, "End of GetFatherAtomicMap.");
  return AtomicMap;
};


void molecule::flipActiveFlag(){
  ActiveFlag = !ActiveFlag;
}

Shape molecule::getBoundingShape() const
{
  // get center and radius
  Vector center;
  double radius = 0.;
  {
    center.Zero();
    for(const_iterator iter = begin(); iter != end(); ++iter)
      center += (*iter)->getPosition();
    center *= 1./(double)size();
    for(const_iterator iter = begin(); iter != end(); ++iter) {
      const Vector &position = (*iter)->getPosition();
      const double temp_distance = position.DistanceSquared(center);
      if (temp_distance > radius)
        radius = temp_distance;
    }
  }
  // convert radius to true value and add some small boundary
  radius = sqrt(radius) + 1e+6*std::numeric_limits<double>::epsilon();
  LOG(1, "INFO: The " << size() << " atoms of the molecule are contained in a sphere at "
      << center << " with radius " << radius << ".");

  Shape BoundingShape(Sphere(center, radius));
  LOG(1, "INFO: Created sphere at " << BoundingShape.getCenter() << " and radius "
      << BoundingShape.getRadius() << ".");
  return BoundingShape;
}

// construct idpool
CONSTRUCT_IDPOOL(atomId_t, continuousId)