source: src/molecule.cpp@ e24c78

/** \file molecules.cpp
 *
 * Functions for the class molecule.
 *
 */

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

#include "Helpers/MemDebug.hpp"

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

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

#include "World.hpp"
#include "atom.hpp"
#include "bond.hpp"
#include "config.hpp"
#include "element.hpp"
#include "graph.hpp"
#include "helpers.hpp"
#include "leastsquaremin.hpp"
#include "linkedcell.hpp"
#include "lists.hpp"
#include "log.hpp"
#include "molecule.hpp"

#include "periodentafel.hpp"
#include "stackclass.hpp"
#include "tesselation.hpp"
#include "vector.hpp"
#include "Matrix.hpp"
#include "World.hpp"
#include "Box.hpp"
#include "Plane.hpp"
#include "Exceptions/LinearDependenceException.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(const periodentafel * const teil) :
  Observable("molecule"),
  elemente(teil),  MDSteps(0),  BondCount(0), NoNonHydrogen(0), NoNonBonds(0),
  NoCyclicBonds(0), BondDistance(0.),  ActiveFlag(false), IndexNr(-1),
  AtomCount(this,boost::bind(&molecule::doCountAtoms,this),"AtomCount"), last_atom(0),  InternalPointer(atoms.begin())
{

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

molecule *NewMolecule(){
  return new molecule(World::getInstance().getPeriode());
}

/** 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(){
  return std::string(name);
}

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

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

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

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

const std::string molecule::getFormula(){
  return formula.toString();
}

std::string molecule::calcFormula(){
  std::map<atomicNumber_t,unsigned int> counts;
  stringstream sstr;
  periodentafel *periode = World::getInstance().getPeriode();
  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    counts[(*iter)->type->getNumber()]++;
  }
  std::map<atomicNumber_t,unsigned int>::reverse_iterator iter;
  for(iter = counts.rbegin(); iter != counts.rend(); ++iter) {
    atomicNumber_t Z = (*iter).first;
    sstr << periode->FindElement(Z)->symbol << (*iter).second;
  }
  return sstr.str();
}

unsigned int molecule::getElementCount(){
  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::iterator molecule::begin(){
  return molecule::iterator(atoms.begin(),this);
}

molecule::const_iterator molecule::begin() const{
  return atoms.begin();
}

molecule::iterator molecule::end(){
  return molecule::iterator(atoms.end(),this);
}

molecule::const_iterator molecule::end() const{
  return atoms.end();
}

bool molecule::empty() const
{
  return (begin() == end());
}

size_t molecule::size() const
{
  size_t counter = 0;
  for (molecule::const_iterator iter = begin(); iter != end (); ++iter)
    counter++;
  return counter;
}

molecule::const_iterator molecule::erase( const_iterator loc )
{
  molecule::const_iterator iter = loc;
  iter--;
  atom* atom = *loc;
  atomIds.erase( atom->getId() );
  atoms.remove( atom );
  atom->removeFromMolecule();
  return iter;
}

molecule::const_iterator molecule::erase( atom * key )
{
  molecule::const_iterator iter = find(key);
  if (iter != end()){
    atomIds.erase( key->getId() );
    atoms.remove( key );
    key->removeFromMolecule();
  }
  return iter;
}

molecule::const_iterator molecule::find ( atom * key ) const
{
  molecule::const_iterator iter;
  for (molecule::const_iterator Runner = begin(); Runner != end(); ++Runner) {
    if (*Runner == key)
      return molecule::const_iterator(Runner);
  }
  return molecule::const_iterator(atoms.end());
}

pair<molecule::iterator,bool> molecule::insert ( atom * const key )
{
  pair<atomIdSet::iterator,bool> res = atomIds.insert(key->getId());
  if (res.second) { // push atom if went well
    atoms.push_back(key);
    return pair<iterator,bool>(molecule::iterator(--end()),res.second);
  } else {
    return pair<iterator,bool>(molecule::iterator(end()),res.second);
  }
}

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

/** 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) {
    pointer->sort = &pointer->nr;
    if (pointer->type != NULL) {
      formula += pointer->type;
      if (pointer->type->Z != 1)
        NoNonHydrogen++;
      if(pointer->getName() == "Unknown"){
        stringstream sstr;
        sstr << pointer->type->symbol << pointer->nr+1;
        pointer->setName(sstr.str());
      }
    }
    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->nr = last_atom++;  // increase number within molecule
    insert(walker);
    if ((pointer->type != NULL) && (pointer->type->Z != 1))
      NoNonHydrogen++;
    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)
{
  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 Matrix &matrix =  World::getInstance().getDomain().getM();
  bond *Binder = NULL;

//  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
  // create vector in direction of bond
  InBondvector = TopReplacement->x - TopOrigin->x;
  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!
  if (bondlength > BondDistance) {
//    Log() << Verbose(4) << "InBondvector is: ";
//    InBondvector.Output(out);
//    Log() << Verbose(0) << endl;
    Orthovector1.Zero();
    for (int i=NDIM;i--;) {
      l = TopReplacement->x[i] - TopOrigin->x[i];
      if (fabs(l) > BondDistance) { // is component greater than bond distance
        Orthovector1[i] = (l < 0) ? -1. : +1.;
      } // (signs are correct, was tested!)
    }
    Orthovector1 *= matrix;
    InBondvector -= Orthovector1; // subtract just the additional translation
    bondlength = InBondvector.Norm();
//    Log() << Verbose(4) << "Corrected InBondvector is now: ";
//    InBondvector.Output(out);
//    Log() << Verbose(0) << endl;
  } // periodic correction finished

  InBondvector.Normalize();
  // get typical bond length and store as scale factor for later
  ASSERT(TopOrigin->type != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
  BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
  if (BondRescale == -1) {
    DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!" << endl);
    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->type = elemente->FindElement(1);  // element is Hydrogen
      FirstOtherAtom->v = TopReplacement->v; // copy velocity
      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
      if (TopReplacement->type->Z == 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->x = TopOrigin->x; // set coordination to origin ...
      FirstOtherAtom->x += InBondvector;  // ... and add distance vector to replacement atom
      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
//      Log() << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
//      FirstOtherAtom->x.Output(out);
//      Log() << Verbose(0) << endl;
      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = TreeEdge;
      break;
    case 2:
      // determine two other bonds (warning if there are more than two other) plus valence sanity check
      for (BondList::const_iterator Runner = TopOrigin->ListOfBonds.begin(); Runner != TopOrigin->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 {
            DoeLog(2) && (eLog()<< Verbose(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() << Verbose(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." << endl;

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

      // create the two Hydrogens ...
      FirstOtherAtom = World::getInstance().createAtom();
      SecondOtherAtom = World::getInstance().createAtom();
      FirstOtherAtom->type = elemente->FindElement(1);
      SecondOtherAtom->type = elemente->FindElement(1);
      FirstOtherAtom->v = TopReplacement->v; // copy velocity
      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
      SecondOtherAtom->v = TopReplacement->v; // copy velocity
      SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
      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->type->HBondAngle[1];
      if (bondangle == -1) {
        DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!" << endl);
        return false;
        bondangle = 0;
      }
      bondangle *= M_PI/180./2.;
//      Log() << Verbose(3) << "ReScaleCheck: InBondvector ";
//      InBondvector.Output(out);
//      Log() << Verbose(0) << endl;
//      Log() << Verbose(3) << "ReScaleCheck: Orthovector ";
//      Orthovector1.Output(out);
//      Log() << Verbose(0) << endl;
//      Log() << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl;
      FirstOtherAtom->x.Zero();
      SecondOtherAtom->x.Zero();
      for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
        FirstOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle));
        SecondOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle));
      }
      FirstOtherAtom->x *= BondRescale;  // rescale by correct BondDistance
      SecondOtherAtom->x *= BondRescale;
      //Log() << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
      for(int i=NDIM;i--;) { // and make relative to origin atom
        FirstOtherAtom->x[i] += TopOrigin->x[i];
        SecondOtherAtom->x[i] += TopOrigin->x[i];
      }
      // ... and add to molecule
      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
      AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
//      Log() << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
//      FirstOtherAtom->x.Output(out);
//      Log() << Verbose(0) << endl;
//      Log() << Verbose(4) << "Added " << *SecondOtherAtom << " at: ";
//      SecondOtherAtom->x.Output(out);
//      Log() << Verbose(0) << endl;
      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = TreeEdge;
      Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = 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->type = elemente->FindElement(1);
      SecondOtherAtom->type = elemente->FindElement(1);
      ThirdOtherAtom->type = elemente->FindElement(1);
      FirstOtherAtom->v = TopReplacement->v; // copy velocity
      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
      SecondOtherAtom->v = TopReplacement->v; // copy velocity
      SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
      ThirdOtherAtom->v = TopReplacement->v; // copy velocity
      ThirdOtherAtom->FixedIon = TopReplacement->FixedIon;
      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() << Verbose(3) << "Orthovector1: ";
//      Orthovector1.Output(out);
//      Log() << Verbose(0) << endl;
      try{
        Orthovector2 = Plane(InBondvector, Orthovector1,0).getNormal();
      }
      catch(LinearDependenceException &excp) {
        Log() << Verbose(0) << excp;
        AllWentWell = false;
      }
//      Log() << Verbose(3) << "Orthovector2: ";
//      Orthovector2.Output(out);
//      Log() << Verbose(0) << endl;

      // create correct coordination for the three atoms
      alpha = (TopOrigin->type->HBondAngle[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() << Verbose(3) << "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", " << endl;
//      Log() << Verbose(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) << endl;
      factors[0] = d;
      factors[1] = f;
      factors[2] = 0.;
      FirstOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
      factors[1] = -0.5*f;
      factors[2] = g;
      SecondOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
      factors[2] = -g;
      ThirdOtherAtom->x.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->x += TopOrigin->x;
      SecondOtherAtom->x += TopOrigin->x;
      ThirdOtherAtom->x += TopOrigin->x;

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

//  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
  return AllWentWell;
};

/** Adds given atom \a *pointer from molecule list.
 * Increases molecule::last_atom and gives last number to added atom.
 * \param filename name and path of xyz file
 * \return true - succeeded, false - file not found
 */
bool molecule::AddXYZFile(string filename)
{

  istringstream *input = NULL;
  int NumberOfAtoms = 0; // atom number in xyz read
  int i, j; // loop variables
  atom *Walker = NULL;  // pointer to added atom
  char shorthand[3];  // shorthand for atom name
  ifstream xyzfile;   // xyz file
  string line;    // currently parsed line
  double x[3];    // atom coordinates

  xyzfile.open(filename.c_str());
  if (!xyzfile)
    return false;

  OBSERVE;
  getline(xyzfile,line,'\n'); // Read numer of atoms in file
  input = new istringstream(line);
  *input >> NumberOfAtoms;
  DoLog(0) && (Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl);
  getline(xyzfile,line,'\n'); // Read comment
  DoLog(1) && (Log() << Verbose(1) << "Comment: " << line << endl);

  if (MDSteps == 0) // no atoms yet present
    MDSteps++;
  for(i=0;i<NumberOfAtoms;i++){
    Walker = World::getInstance().createAtom();
    getline(xyzfile,line,'\n');
    istringstream *item = new istringstream(line);
    //istringstream input(line);
    //Log() << Verbose(1) << "Reading: " << line << endl;
    *item >> shorthand;
    *item >> x[0];
    *item >> x[1];
    *item >> x[2];
    Walker->type = elemente->FindElement(shorthand);
    if (Walker->type == NULL) {
      DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.");
      Walker->type = elemente->FindElement(1);
    }
    if (Walker->Trajectory.R.size() <= (unsigned int)MDSteps) {
      Walker->Trajectory.R.resize(MDSteps+10);
      Walker->Trajectory.U.resize(MDSteps+10);
      Walker->Trajectory.F.resize(MDSteps+10);
    }
    for(j=NDIM;j--;) {
      Walker->x[j] = x[j];
      Walker->Trajectory.R.at(MDSteps-1)[j] = x[j];
      Walker->Trajectory.U.at(MDSteps-1)[j] = 0;
      Walker->Trajectory.F.at(MDSteps-1)[j] = 0;
    }
    AddAtom(Walker);  // add to molecule
    delete(item);
  }
  xyzfile.close();
  delete(input);
  return true;
};

/** Creates a copy of this molecule.
 * \return copy of molecule
 */
molecule *molecule::CopyMolecule()
{
  molecule *copy = World::getInstance().createMolecule();
  atom *LeftAtom = NULL, *RightAtom = NULL;

  // copy all atoms
  ActOnCopyWithEachAtom ( &molecule::AddCopyAtom, copy );

  // copy all bonds
  bond *Binder = NULL;
  bond *NewBond = NULL;
  for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
    for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); !(*AtomRunner)->ListOfBonds.empty(); BondRunner = (*AtomRunner)->ListOfBonds.begin())
      if ((*BondRunner)->leftatom == *AtomRunner) {
        Binder = (*BondRunner);

        // get the pendant atoms of current bond in the copy molecule
        copy->ActOnAllAtoms( &atom::EqualsFather, (const atom *)Binder->leftatom, (const atom **)&LeftAtom );
        copy->ActOnAllAtoms( &atom::EqualsFather, (const atom *)Binder->rightatom, (const atom **)&RightAtom );

        NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
        NewBond->Cyclic = Binder->Cyclic;
        if (Binder->Cyclic)
          copy->NoCyclicBonds++;
        NewBond->Type = Binder->Type;
      }
  // correct fathers
  ActOnAllAtoms( &atom::CorrectFather );

  // copy values
  if (hasBondStructure()) {  // if adjaceny list is present
    copy->BondDistance = BondDistance;
  }

  return copy;
};


/**
 * 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();

  BOOST_FOREACH(atom *iter,atoms){
    if(iter->IsInShape(region)){
      copy->AddCopyAtom(iter);
    }
  }

  //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, "First atom in bond-creation was an invalid pointer");
  ASSERT(atom2, "Second atom in bond-creation was an invalid pointer");
  ASSERT(FindAtom(atom1->nr),"First atom in bond-creation was not part of molecule");
  ASSERT(FindAtom(atom2->nr),"Second atom in bond-creation was not parto of molecule");

  Binder = new bond(atom1, atom2, degree, BondCount++);
  atom1->RegisterBond(Binder);
  atom2->RegisterBond(Binder);
  if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
    NoNonBonds++;

  return Binder;
};

/** Remove bond from bond chain list and from the both atom::ListOfBonds.
 * \todo Function not implemented yet
 * \param *pointer bond pointer
 * \return true - bound found and removed, false - bond not found/removed
 */
bool molecule::RemoveBond(bond *pointer)
{
  //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
  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)
{
  //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
  BondList::const_iterator ForeRunner;
  while (!BondPartner->ListOfBonds.empty()) {
    ForeRunner = BondPartner->ListOfBonds.begin();
    RemoveBond(*ForeRunner);
  }
  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)
{
  Matrix 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 deletes it.
 * \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;
  formula-=pointer->type;
  RemoveBonds(pointer);
  erase(pointer);
  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;
  formula-=pointer->type;
  erase(pointer);
  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())
      erase(iter);
  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::const_iterator iter = begin();
  for (; iter != end(); ++iter)
    if ((*iter)->nr == Nr)
      break;
  if (iter != end()) {
    //Log() << Verbose(0) << "Found Atom Nr. " << walker->nr << endl;
    return (*iter);
  } else {
    DoLog(0) && (Log() << Verbose(0) << "Atom not found in list." << endl);
    return NULL;
  }
};

/** Asks for atom number, and checks whether in list.
 * \param *text question before entering
 */
atom * molecule::AskAtom(string text)
{
  int No;
  atom *ion = NULL;
  do {
    //Log() << Verbose(0) << "============Atom list==========================" << endl;
    //mol->Output((ofstream *)&cout);
    //Log() << Verbose(0) << "===============================================" << endl;
    DoLog(0) && (Log() << Verbose(0) << 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 Matrix &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(ofstream * const output)
{
  int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS];

  for (int i=0;i<MAX_ELEMENTS;++i) {
    AtomNo[i] = 0;
    ElementNo[i] = 0;
  }
  if (output == NULL) {
    return false;
  } else {
    *output << "#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)" << endl;
    SetIndexedArrayForEachAtomTo ( ElementNo, &element::Z, &AbsoluteValue, 1);
    int current=1;
    for (int i=0;i<MAX_ELEMENTS;++i) {
      if (ElementNo[i] == 1)
        ElementNo[i] = current++;
    }
    ActOnAllAtoms( &atom::OutputArrayIndexed, (ostream * const) output, (const int *)ElementNo, (int *)AtomNo, (const char *) NULL );
    return true;
  }
};

/** Prints molecule with all atomic trajectory positions to *out.
 * \param *out output stream
 */
bool molecule::OutputTrajectories(ofstream * const output)
{
  int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS];

  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;
      }
      for (int i=0;i<MAX_ELEMENTS;++i) {
        AtomNo[i] = 0;
        ElementNo[i] = 0;
      }
      SetIndexedArrayForEachAtomTo ( ElementNo, &element::Z, &AbsoluteValue, 1);
      int current=1;
      for (int i=0;i<MAX_ELEMENTS;++i) {
        if (ElementNo[i] == 1)
          ElementNo[i] = current++;
      }
      ActOnAllAtoms( &atom::OutputTrajectory, output, (const int *)ElementNo, AtomNo, (const int)step );
    }
    return true;
  }
};

/** Outputs contents of each atom::ListOfBonds.
 * \param *out output stream
 */
void molecule::OutputListOfBonds() const
{
  DoLog(2) && (Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl);
  ActOnAllAtoms (&atom::OutputBondOfAtom );
  DoLog(0) && (Log() << Verbose(0) << endl);
};

/** 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);
      ActOnAllAtoms( &atom::OutputTrajectoryXYZ, 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);
    ActOnAllAtoms( &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();
  int i = 0;
  NoNonHydrogen = 0;
  for (molecule::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    (*iter)->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
    if ((*iter)->type->Z != 1) // count non-hydrogen atoms whilst at it
      NoNonHydrogen++;
    stringstream sstr;
    sstr << (*iter)->type->symbol << (*iter)->nr+1;
    (*iter)->setName(sstr.str());
    DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->getName() << "." << endl);
    i++;
  }
  return res;
};

/** Determines whether two molecules actually contain the same atoms and coordination.
 * \param *out output stream for debugging
 * \param *OtherMolecule the molecule to compare this one to
 * \param threshold upper limit of difference when comparing the coordination.
 * \return NULL - not equal, otherwise an allocated (molecule::AtomCount) permutation map of the atom numbers (which corresponds to which)
 */
int * molecule::IsEqualToWithinThreshold(molecule *OtherMolecule, double threshold)
{
  int flag;
  double *Distances = NULL, *OtherDistances = NULL;
  Vector CenterOfGravity, OtherCenterOfGravity;
  size_t *PermMap = NULL, *OtherPermMap = NULL;
  int *PermutationMap = NULL;
  bool result = true; // status of comparison

  DoLog(3) && (Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl);
  /// first count both their atoms and elements and update lists thereby ...
  //Log() << Verbose(0) << "Counting atoms, updating list" << endl;

  /// ... and compare:
  /// -# AtomCount
  if (result) {
    if (getAtomCount() != OtherMolecule->getAtomCount()) {
      DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl);
      result = false;
    } else Log() << Verbose(4) << "AtomCounts match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;
  }
  /// -# Formula
  if (result) {
    if (formula != OtherMolecule->formula) {
      DoLog(4) && (Log() << Verbose(4) << "Formulas don't match: " << formula << " == " << OtherMolecule->formula << endl);
      result = false;
    } else Log() << Verbose(4) << "Formulas match: " << formula << " == " << OtherMolecule->formula << endl;
  }
  /// then determine and compare center of gravity for each molecule ...
  if (result) {
    DoLog(5) && (Log() << Verbose(5) << "Calculating Centers of Gravity" << endl);
    DeterminePeriodicCenter(CenterOfGravity);
    OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
    DoLog(5) && (Log() << Verbose(5) << "Center of Gravity: " << CenterOfGravity << endl);
    DoLog(5) && (Log() << Verbose(5) << "Other Center of Gravity: " << OtherCenterOfGravity << endl);
    if (CenterOfGravity.DistanceSquared(OtherCenterOfGravity) > threshold*threshold) {
      DoLog(4) && (Log() << Verbose(4) << "Centers of gravity don't match." << endl);
      result = false;
    }
  }

  /// ... then make a list with the euclidian distance to this center for each atom of both molecules
  if (result) {
    DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);
    Distances = new double[getAtomCount()];
    OtherDistances = new double[getAtomCount()];
    SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
    SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
    for(int i=0;i<getAtomCount();i++) {
      Distances[i] = 0.;
      OtherDistances[i] = 0.;
    }

    /// ... sort each list (using heapsort (o(N log N)) from GSL)
    DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl);
    PermMap = new size_t[getAtomCount()];
    OtherPermMap = new size_t[getAtomCount()];
    for(int i=0;i<getAtomCount();i++) {
      PermMap[i] = 0;
      OtherPermMap[i] = 0;
    }
    gsl_heapsort_index (PermMap, Distances, getAtomCount(), sizeof(double), CompareDoubles);
    gsl_heapsort_index (OtherPermMap, OtherDistances, getAtomCount(), sizeof(double), CompareDoubles);
    PermutationMap = new int[getAtomCount()];
    for(int i=0;i<getAtomCount();i++)
      PermutationMap[i] = 0;
    DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl);
    for(int i=getAtomCount();i--;)
      PermutationMap[PermMap[i]] = (int) OtherPermMap[i];

    /// ... and compare them step by step, whether the difference is individually(!) below \a threshold for all
    DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl);
    flag = 0;
    for (int i=0;i<getAtomCount();i++) {
      DoLog(5) && (Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl);
      if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold)
        flag = 1;
    }

    // free memory
    delete[](PermMap);
    delete[](OtherPermMap);
    delete[](Distances);
    delete[](OtherDistances);
    if (flag) { // if not equal
      delete[](PermutationMap);
      result = false;
    }
  }
  /// return pointer to map if all distances were below \a threshold
  DoLog(3) && (Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl);
  if (result) {
    DoLog(3) && (Log() << Verbose(3) << "Result: Equal." << endl);
    return PermutationMap;
  } else {
    DoLog(3) && (Log() << Verbose(3) << "Result: Not equal." << endl);
    return NULL;
  }
};

/** 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)
{
  DoLog(3) && (Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl);
  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;
    DoLog(4) && (Log() << Verbose(4) << "Map is trivial." << endl);
  } else {
    DoLog(4) && (Log() << Verbose(4) << "Map is ");
    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
      if ((*iter)->father == NULL) {
        AtomicMap[(*iter)->nr] = -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() << Verbose(4) << "Comparing father " << (*iter)->father << " with the other one " << (*runner)->father << "." << endl;
          if ((*iter)->father == (*runner))
            AtomicMap[(*iter)->nr] = (*runner)->nr;
        }
      }
      DoLog(0) && (Log() << Verbose(0) << AtomicMap[(*iter)->nr] << "\t");
    }
    DoLog(0) && (Log() << Verbose(0) << endl);
  }
  DoLog(3) && (Log() << Verbose(3) << "End of GetFatherAtomicMap." << endl);
  return AtomicMap;
};

/** Stores the temperature evaluated from velocities in molecule::Trajectories.
 * We simply use the formula equivaleting temperature and kinetic energy:
 * \f$k_B T = \sum_i m_i v_i^2\f$
 * \param *output output stream of temperature file
 * \param startstep first MD step in molecule::Trajectories
 * \param endstep last plus one MD step in molecule::Trajectories
 * \return file written (true), failure on writing file (false)
 */
bool molecule::OutputTemperatureFromTrajectories(ofstream * const output, int startstep, int endstep)
{
  double temperature;
  // test stream
  if (output == NULL)
    return false;
  else
    *output << "# Step Temperature [K] Temperature [a.u.]" << endl;
  for (int step=startstep;step < endstep; step++) { // loop over all time steps
    temperature = 0.;
    ActOnAllAtoms( &TrajectoryParticle::AddKineticToTemperature, &temperature, step);
    *output << step << "\t" << temperature*AtomicEnergyToKelvin << "\t" << temperature << endl;
  }
  return true;
};

void molecule::SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::*index) const
{
  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    array[((*iter)->*index)] = (*iter);
  }
};

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