source: src/molecule.cpp@ e38750

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

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

#include "CodePatterns/MemDebug.hpp"

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

#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 "CodePatterns/Observer/Observable.hpp"
#include "CodePatterns/Observer/Notification.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"

using namespace boost::assign;

// static entities
static Observable::channels_t getBoundingBoxChannels()
{
  Observable::channels_t channels;
  channels += molecule::AtomInserted, molecule::AtomRemoved, molecule::AtomMoved;
  return channels;
}

/************************************* 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),
  NoNonBonds(0),
  NoCyclicBonds(0),
  ActiveFlag(false),
  IndexNr(-1),
  NoNonHydrogen(this,boost::bind(&molecule::doCountNoNonHydrogen,this),"NoNonHydrogen"),
  BondCount(this,boost::bind(&molecule::doCountBonds,this),"BondCount"),
  atomIdPool(1, 20, 100),
  BoundingBoxSweepingAxis(std::vector<AtomDistanceMap_t>(NDIM)),
  _lastchangedatomid(-1),
  last_atom(0)
{
  // add specific channels
  Channels *OurChannel = new Channels;
  Observable::insertNotificationChannel( std::make_pair( static_cast<Observable *>(this), OurChannel) );
  for (size_t type = 0; type < (size_t)NotificationType_MAX; ++type)
    OurChannel->addChannel(type);

  // cannot initialize in initializer body as then channels have not been setup yet
  BoundingBox.reset(
      new Cacheable<BoundingBoxInfo>(
          this, boost::bind(&molecule::updateBoundingBox, boost::cref(this)), "molecule_BoundingBox", getBoundingBoxChannels()));

  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()
{
  // inform all UI elements about imminent removal before anything is lost
  {
    OBSERVE;
    NOTIFY(AboutToBeRemoved);
  }
  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 atomIds.size();
}

size_t molecule::getNoNonHydrogen() const{
  return *NoNonHydrogen;
}

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

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

void molecule::InsertLocalToGlobalId(atom * const pointer)
{
#ifndef NDEBUG
  std::pair< LocalToGlobalId_t::iterator, bool > inserter =
#endif
      LocalToGlobalId.insert( std::make_pair(pointer->getNr(), pointer) );
  ASSERT( inserter.second,
      "molecule::AddAtom() - local number "+toString(pointer->getNr())+" appears twice.");
}

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

bool molecule::changeAtomId(int oldId, int newId)
{
  OBSERVE;
  if ((!atomIds.contains( oldId )) || (atomIds.contains( newId )))
    return false;
  atomIds.erase( oldId );
  atomIds.insert( newId );
  // also update BoundingBoxSweepingAxis
  for (int i=0;i<NDIM;++i) {
    AtomDistanceMap_t::left_iterator iter = BoundingBoxSweepingAxis[i].left.find(oldId);
    ASSERT(iter != BoundingBoxSweepingAxis[i].left.end(),
        "molecule::changeAtomId() - could not find atom "+toString(oldId)
        +" in BoundingBoxSweepingAxis.");
    const double component = iter->second;
    BoundingBoxSweepingAxis[i].left.erase(iter);
    BoundingBoxSweepingAxis[i].left.insert( std::make_pair(newId, component) );
  }
  return true;
}

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)){
    OBSERVE;
    NOTIFY(IndexChanged);
    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);
  {
    _lastchangedatomid = _atom->getId();
    NOTIFY(AtomRemoved);
  }
  atomIds.erase( _atom->getId() );
  {
    BoundingBoxInfo oldinfo = updateBoundingBox();
    for (int i=0;i<NDIM;++i)
      BoundingBoxSweepingAxis[i].left.erase( _atom->getId() );
    BoundingBoxInfo newinfo = updateBoundingBox();
    if (oldinfo != newinfo)
      NOTIFY(BoundingBoxChanged);
  }
  {
    NOTIFY(AtomNrChanged);
    atomIdPool.releaseId(_atom->getNr());
    LocalToGlobalId.erase(_atom->getNr());
    _atom->setNr(-1);
  }
  NOTIFY(FormulaChanged);
  formula-=_atom->getType();
  _atom->removeFromMolecule();
  return iter;
}

molecule::const_iterator molecule::erase( atom * key )
{
  const_iterator iter = const_cast<const molecule &>(*this).find(key);
  if (iter != const_cast<const molecule &>(*this).end())
    return erase(iter);
  else
    return iter;
}

pair<molecule::iterator,bool> molecule::insert ( atom * const key )
{
  OBSERVE;
  NOTIFY(AtomInserted);
  _lastchangedatomid = key->getId();
  std::pair<iterator,bool> res = atomIds.insert(key->getId());
  if (res.second) { // push atom if went well
    {
      BoundingBoxInfo oldinfo = updateBoundingBox();
      for (int i=0;i<NDIM;++i)
        BoundingBoxSweepingAxis[i].left.insert( std::make_pair(key->getId(), key->getPosition()[i]));
      BoundingBoxInfo newinfo = updateBoundingBox();
      if (oldinfo != newinfo)
        NOTIFY(BoundingBoxChanged);
    }
    NOTIFY(AtomNrChanged);
    key->setNr(atomIdPool.getNextId());
    InsertLocalToGlobalId(key);
    setAtomName(key);
    NOTIFY(FormulaChanged);
    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());
}

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

World::ConstAtomComposite molecule::getAtomSet() const
{
  World::ConstAtomComposite vector_of_atoms;
  for (molecule::const_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)
{
  if (pointer != NULL) {
    // molecule::insert() is called by setMolecule()
    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;
  if (pointer != NULL) {
    atom *walker = pointer->clone();
    AddAtom(walker);
    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::ptr 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
//  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::ptr FirstBond;
//  bond::ptr SecondBond; // 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::ptr Binder;
//
//  // 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->getDegree()-1);
//  if (BondRescale == -1) {
//    ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->getDegree() << "!");
//    return false;
//    BondRescale = bondlength;
//  } else {
//    if (!IsAngstroem)
//      BondRescale /= (1.*AtomicLengthToAngstroem);
//  }
//
//  // discern single, double and triple bonds
//  switch(TopBond->getDegree()) {
//    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->getDegree() << "!");
//        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)
{
  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 = const_cast<const molecule &>(*this).begin();
      AtomRunner != const_cast<const molecule &>(*this).end();
      ++AtomRunner) {
    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
        BondRunner != ListOfBonds.end();
        ++BondRunner)
      if ((*BondRunner)->leftatom == *AtomRunner) {
        bond::ptr 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::ptr const NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->getDegree());
        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 removeAtomsinMolecule(molecule *&_mol)
{
  // copy list of atoms from molecule as it will be changed
  std::vector<atom *> atoms;
  atoms.resize(_mol->getAtomCount(), NULL);
  std::copy(_mol->begin(), _mol->end(), atoms.begin());
  // remove each atom from world
  for(std::vector<atom *>::iterator AtomRunner = atoms.begin();
      AtomRunner != atoms.end(); ++AtomRunner)
    World::getInstance().destroyAtom(*AtomRunner);
  // make sure that pointer os not usable
  _mol = NULL;
};


/**
 * 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) {
  molecule *copy = World::getInstance().createMolecule();

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

  // copy all bonds
  for(molecule::const_iterator AtomRunner = const_cast<const molecule &>(*this).begin();
      AtomRunner != const_cast<const molecule &>(*this).end();
      ++AtomRunner) {
    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
        BondRunner != ListOfBonds.end();
        ++BondRunner)
      if ((*BondRunner)->leftatom == *AtomRunner) {
        bond::ptr 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::ptr const NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->getDegree());
          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::ptr molecule::AddBond(atom *atom1, atom *atom2, int degree)
{
  bond::ptr Binder;

  // 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.reset(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;
};

/** Set molecule::name from the basename without suffix in the given \a *filename.
 * \param *filename filename
 */
void molecule::SetNameFromFilename(const char *filename)
{
  OBSERVE;
  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';
};

/** 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
{
  LocalToGlobalId_t::const_iterator iter = LocalToGlobalId.find(Nr);
  if (iter != LocalToGlobalId.end()) {
    //LOG(0, "Found Atom Nr. " << walker->getNr());
    return iter->second;
  } else {
    ELOG(1, "Atom with Nr " << Nr << " 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) const
{
  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;
  }
};

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

/** Brings molecule::AtomCount and atom::*Name up-to-date.
 * \param *out output stream for debugging
 */
size_t molecule::doCountNoNonHydrogen() const
{
  int temp = 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
      ++temp;
  return temp;
};

/** 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(const molecule * const 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 = const_cast<const molecule &>(*this).begin();
        iter != const_cast<const molecule &>(*this).end();
        ++iter) {
      if ((*iter)->getFather() == 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)->getFather() << " with the other one " << (*runner)->getFather() << ".");
          if ((*iter)->getFather() == (*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 double scale) const
{
  // create Sphere around every atom
  if (empty())
    return Nowhere();
  const_iterator iter = begin();
  const Vector center = (*iter)->getPosition();
  const double vdWRadius = (*iter)->getElement().getVanDerWaalsRadius();
  Shape BoundingShape = Sphere(center, vdWRadius*scale);
  for(++iter; iter != end(); ++iter) {
    const Vector center = (*iter)->getPosition();
    const double vdWRadius = (*iter)->getElement().getVanDerWaalsRadius();
    if (vdWRadius*scale != 0.)
      BoundingShape = Sphere(center, vdWRadius*scale) || BoundingShape;
  }
  return BoundingShape;
}

Shape molecule::getBoundingSphere(const double boundary) const
{
  // get center and radius
  Vector center;
  double radius = 0.;
  {
    center.Zero();
    for(const_iterator iter = begin(); iter != end(); ++iter)
      center += (*iter)->getPosition();
    if (begin() != end())
      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) + boundary + 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 << ".");

  // TODO: When we do not use a Sphere here anymore, then FillRegularGridAction will
  // will not work as it expects a sphere due to possible random rotations.
  Shape BoundingShape(Sphere(center, radius));
  LOG(1, "INFO: Created sphere at " << BoundingShape.getCenter() << " and radius "
      << BoundingShape.getRadius() << ".");
  return BoundingShape;
}

molecule::BoundingBoxInfo molecule::updateBoundingBox() const
{
  BoundingBoxInfo info;
  Vector min = zeroVec;
  Vector max = zeroVec;
  for (int i=0;i<NDIM;++i) {
    if (!BoundingBoxSweepingAxis[i].right.empty()) {
      min[i] = BoundingBoxSweepingAxis[i].right.begin()->first;
      max[i] = BoundingBoxSweepingAxis[i].right.rbegin()->first;
    }
  }
  info.radius = (.5*(max-min)).Norm();
  info.position = .5*(max+min);
  return info;
}

molecule::BoundingBoxInfo molecule::getBoundingBox() const
{
  return **BoundingBox;
}

void molecule::update(Observable *publisher)
{
  ASSERT(0, "molecule::update() - did not sign on for any general updates.");
}

void molecule::recieveNotification(Observable *publisher, Notification_ptr notification)
{
  const atom * const _atom = dynamic_cast<atom *>(publisher);
  if ((_atom != NULL) && containsAtom(_atom)) {
#ifdef LOG_OBSERVER
    observerLog().addMessage() << "++ Update of Observer "<< observerLog().getName(static_cast<Observer *>(this))
          << " received notification from atom " << _atom->getId() << " for channel "
          << notification->getChannelNo() << ".";
#endif
    switch (notification->getChannelNo()) {
      case AtomObservable::PositionChanged:
      {
        // emit others about one of our atoms moved
        _lastchangedatomid = _atom->getId();
        // update entry in map
        BoundingBoxInfo oldinfo = updateBoundingBox();
        for (int i=0;i<NDIM;++i) {
          AtomDistanceMap_t::left_iterator iter = BoundingBoxSweepingAxis[i].left.find(_atom->getId());
          ASSERT(iter != BoundingBoxSweepingAxis[i].left.end(),
              "molecule::recieveNotification() - could not find atom "+toString(_atom->getId())
              +" in BoundingBoxSweepingAxis.");
          BoundingBoxSweepingAxis[i].left.erase(iter);
          BoundingBoxSweepingAxis[i].left.insert(
              std::make_pair(_atom->getId(), _atom->getPosition()[i]) );
        }
        BoundingBoxInfo newinfo = updateBoundingBox();
        OBSERVE;
        NOTIFY(AtomMoved);
        if (oldinfo != newinfo)
          NOTIFY(BoundingBoxChanged);
        break;
      }
      case AtomObservable::ElementChanged:
      {
        // emit others about one of our atoms moved
        _lastchangedatomid = _atom->getId();
        OBSERVE;
        NOTIFY(FormulaChanged);
        resetFormula();
        break;
      }
      default:
        ASSERT( 0, "molecule::recieveNotification() - we did not sign up for channel "
            +toString(notification->getChannelNo()));
        break;
    }
  }
}

void molecule::subjectKilled(Observable *publisher)
{
  // do nothing, atom does it all
}

void molecule::resetFormula()
{
  // clear
  formula.clear();

  for (molecule::const_iterator iter = const_cast<const molecule *>(this)->begin();
      iter != const_cast<const molecule *>(this)->end(); ++iter)
    formula+=(*iter)->getType();
}

// construct idpool
CONSTRUCT_IDPOOL(atomId_t, continuousId)