Index: src/Fragmentation/Exporters/ExportGraph.cpp
===================================================================
--- src/Fragmentation/Exporters/ExportGraph.cpp	(revision 9fd44f7f53b271add9c0ff9fe035e67c63361f09)
+++ src/Fragmentation/Exporters/ExportGraph.cpp	(revision 06804b923c9b82a6dd4f7d932a21781d6e4eae6e)
@@ -137,129 +137,129 @@
 }
 
-/** Internal helper to create from each keyset a molecule
- *
- */
-void ExportGraph::prepareMolecule()
-{
-  size_t count = 0;
-  for(Graph::const_iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
-    KeySet test = (*runner).first;
-    LOG(2, "DEBUG: Fragment No." << (*runner).second.first << " with TEFactor "
-        << (*runner).second.second << ".");
-    BondFragments.insert(StoreFragmentFromKeySet(test, World::getInstance().getConfig()));
-    ++count;
-  }
-  LOG(1, "INFO: " << count << "/" << BondFragments.ListOfMolecules.size()
-      << " fragments generated from the keysets.");
-}
-
-/** Stores a fragment from \a KeySet into \a molecule.
- * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
- * molecule and adds missing hydrogen where bonds were cut.
- * \param &Leaflet pointer to KeySet structure
- * \param IsAngstroem whether we have Ansgtroem or bohrradius
- * \return pointer to constructed molecule
- */
-molecule * ExportGraph::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
-{
-  Info info(__func__);
-  ListOfLocalAtoms_t SonList;
-  molecule *Leaf = World::getInstance().createMolecule();
-
-  StoreFragmentFromKeySet_Init(Leaf, Leaflet, SonList);
-  // create the bonds between all: Make it an induced subgraph and add hydrogen
-//  LOG(2, "Creating bonds from father graph (i.e. induced subgraph creation).");
-  CreateInducedSubgraphOfFragment(Leaf, SonList, IsAngstroem);
-
-  //Leaflet->Leaf->ScanForPeriodicCorrection(out);
-  return Leaf;
-}
-
-/** Initializes some value for putting fragment of \a *mol into \a *Leaf.
- * \param *Leaf fragment molecule
- * \param &Leaflet pointer to KeySet structure
- * \param SonList calloc'd list which atom of \a *Leaf is a son of which atom in \a *mol
- * \return number of atoms in fragment
- */
-int ExportGraph::StoreFragmentFromKeySet_Init(molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList)
-{
-  atom *FatherOfRunner = NULL;
-
-  // first create the minimal set of atoms from the KeySet
-  World &world = World::getInstance();
-  int size = 0;
-  for(KeySet::const_iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
-    FatherOfRunner = world.getAtom(AtomById(*runner));  // find the id
-    SonList.insert( std::make_pair(FatherOfRunner->getNr(), Leaf->AddCopyAtom(FatherOfRunner) ) );
-    size++;
-  }
-  return size;
-}
-
-/** Creates an induced subgraph out of a fragmental key set, adding bonds and hydrogens (if treated specially).
- * \param *Leaf fragment molecule
- * \param IsAngstroem whether we have Ansgtroem or bohrradius
- * \param SonList list which atom of \a *Leaf is another atom's son
- */
-void ExportGraph::CreateInducedSubgraphOfFragment(molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem)
-{
-  bool LonelyFlag = false;
-  atom *OtherFather = NULL;
-  atom *FatherOfRunner = NULL;
-
-  // we increment the iter just before skipping the hydrogen
-  // as we use AddBond, we cannot have a const_iterator here
-  for (molecule::iterator iter = Leaf->begin(); iter != Leaf->end();) {
-    LonelyFlag = true;
-    FatherOfRunner = (*iter)->father;
-    ASSERT(FatherOfRunner,"Atom without father found");
-    if (SonList.find(FatherOfRunner->getNr()) != SonList.end())  {  // check if this, our father, is present in list
-      // create all bonds
-      const BondList& ListOfBonds = FatherOfRunner->getListOfBonds();
-      for (BondList::const_iterator BondRunner = ListOfBonds.begin();
-          BondRunner != ListOfBonds.end();
-          ++BondRunner) {
-        OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner);
-        if (SonList.find(OtherFather->getNr()) != SonList.end()) {
-//          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
-//              << " is bound to " << *OtherFather << ", whose son is "
-//              << *SonList[OtherFather->getNr()] << ".");
-          if (OtherFather->getNr() > FatherOfRunner->getNr()) { // add bond (Nr check is for adding only one of both variants: ab, ba)
-            std::stringstream output;
-//            output << "ACCEPT: Adding Bond: "
-            output << Leaf->AddBond((*iter), SonList[OtherFather->getNr()], (*BondRunner)->getDegree());
-//            LOG(3, output.str());
-            //NumBonds[(*iter)->getNr()]++;
-          } else {
-//            LOG(3, "REJECY: Not adding bond, labels in wrong order.");
-          }
-          LonelyFlag = false;
-        } else {
-//          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
-//              << " is bound to " << *OtherFather << ", who has no son in this fragment molecule.");
-          if (saturation == DoSaturate) {
-//          LOG(3, "ACCEPT: Adding Hydrogen to " << (*iter)->Name << " and a bond in between.");
-            if (!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem))
-              exit(1);
-          } else if ((treatment == ExcludeHydrogen) && (OtherFather->getElementNo() == (atomicNumber_t)1)) {
-            // just copy the atom if it's a hydrogen
-            atom * const OtherWalker = Leaf->AddCopyAtom(OtherFather);
-            Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
-          }
-          //NumBonds[(*iter)->getNr()] += Binder->getDegree();
-        }
-      }
-    } else {
-      ELOG(1, "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->getNr()] << "!");
-    }
-    if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) {
-      LOG(0, **iter << "has got bonds only to hydrogens!");
-    }
-    ++iter;
-    if (saturation == DoSaturate) {
-      while ((iter != Leaf->end()) && ((*iter)->getType()->getAtomicNumber() == 1)){ // skip added hydrogen
-        iter++;
-      }
-    }
-  }
-}
+///** Internal helper to create from each keyset a molecule
+// *
+// */
+//void ExportGraph::prepareMolecule()
+//{
+//  size_t count = 0;
+//  for(Graph::const_iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
+//    KeySet test = (*runner).first;
+//    LOG(2, "DEBUG: Fragment No." << (*runner).second.first << " with TEFactor "
+//        << (*runner).second.second << ".");
+//    BondFragments.insert(StoreFragmentFromKeySet(test, World::getInstance().getConfig()));
+//    ++count;
+//  }
+//  LOG(1, "INFO: " << count << "/" << BondFragments.ListOfMolecules.size()
+//      << " fragments generated from the keysets.");
+//}
+//
+///** Stores a fragment from \a KeySet into \a molecule.
+// * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
+// * molecule and adds missing hydrogen where bonds were cut.
+// * \param &Leaflet pointer to KeySet structure
+// * \param IsAngstroem whether we have Ansgtroem or bohrradius
+// * \return pointer to constructed molecule
+// */
+//molecule * ExportGraph::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
+//{
+//  Info info(__func__);
+//  ListOfLocalAtoms_t SonList;
+//  molecule *Leaf = World::getInstance().createMolecule();
+//
+//  StoreFragmentFromKeySet_Init(Leaf, Leaflet, SonList);
+//  // create the bonds between all: Make it an induced subgraph and add hydrogen
+////  LOG(2, "Creating bonds from father graph (i.e. induced subgraph creation).");
+//  CreateInducedSubgraphOfFragment(Leaf, SonList, IsAngstroem);
+//
+//  //Leaflet->Leaf->ScanForPeriodicCorrection(out);
+//  return Leaf;
+//}
+//
+///** Initializes some value for putting fragment of \a *mol into \a *Leaf.
+// * \param *Leaf fragment molecule
+// * \param &Leaflet pointer to KeySet structure
+// * \param SonList calloc'd list which atom of \a *Leaf is a son of which atom in \a *mol
+// * \return number of atoms in fragment
+// */
+//int ExportGraph::StoreFragmentFromKeySet_Init(molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList)
+//{
+//  atom *FatherOfRunner = NULL;
+//
+//  // first create the minimal set of atoms from the KeySet
+//  World &world = World::getInstance();
+//  int size = 0;
+//  for(KeySet::const_iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
+//    FatherOfRunner = world.getAtom(AtomById(*runner));  // find the id
+//    SonList.insert( std::make_pair(FatherOfRunner->getNr(), Leaf->AddCopyAtom(FatherOfRunner) ) );
+//    size++;
+//  }
+//  return size;
+//}
+//
+///** Creates an induced subgraph out of a fragmental key set, adding bonds and hydrogens (if treated specially).
+// * \param *Leaf fragment molecule
+// * \param IsAngstroem whether we have Ansgtroem or bohrradius
+// * \param SonList list which atom of \a *Leaf is another atom's son
+// */
+//void ExportGraph::CreateInducedSubgraphOfFragment(molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem)
+//{
+//  bool LonelyFlag = false;
+//  atom *OtherFather = NULL;
+//  atom *FatherOfRunner = NULL;
+//
+//  // we increment the iter just before skipping the hydrogen
+//  // as we use AddBond, we cannot have a const_iterator here
+//  for (molecule::iterator iter = Leaf->begin(); iter != Leaf->end();) {
+//    LonelyFlag = true;
+//    FatherOfRunner = (*iter)->father;
+//    ASSERT(FatherOfRunner,"Atom without father found");
+//    if (SonList.find(FatherOfRunner->getNr()) != SonList.end())  {  // check if this, our father, is present in list
+//      // create all bonds
+//      const BondList& ListOfBonds = FatherOfRunner->getListOfBonds();
+//      for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+//          BondRunner != ListOfBonds.end();
+//          ++BondRunner) {
+//        OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner);
+//        if (SonList.find(OtherFather->getNr()) != SonList.end()) {
+////          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
+////              << " is bound to " << *OtherFather << ", whose son is "
+////              << *SonList[OtherFather->getNr()] << ".");
+//          if (OtherFather->getNr() > FatherOfRunner->getNr()) { // add bond (Nr check is for adding only one of both variants: ab, ba)
+//            std::stringstream output;
+////            output << "ACCEPT: Adding Bond: "
+//            output << Leaf->AddBond((*iter), SonList[OtherFather->getNr()], (*BondRunner)->getDegree());
+////            LOG(3, output.str());
+//            //NumBonds[(*iter)->getNr()]++;
+//          } else {
+////            LOG(3, "REJECY: Not adding bond, labels in wrong order.");
+//          }
+//          LonelyFlag = false;
+//        } else {
+////          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
+////              << " is bound to " << *OtherFather << ", who has no son in this fragment molecule.");
+//          if (saturation == DoSaturate) {
+////          LOG(3, "ACCEPT: Adding Hydrogen to " << (*iter)->Name << " and a bond in between.");
+//            if (!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem))
+//              exit(1);
+//          } else if ((treatment == ExcludeHydrogen) && (OtherFather->getElementNo() == (atomicNumber_t)1)) {
+//            // just copy the atom if it's a hydrogen
+//            atom * const OtherWalker = Leaf->AddCopyAtom(OtherFather);
+//            Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
+//          }
+//          //NumBonds[(*iter)->getNr()] += Binder->getDegree();
+//        }
+//      }
+//    } else {
+//      ELOG(1, "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->getNr()] << "!");
+//    }
+//    if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) {
+//      LOG(0, **iter << "has got bonds only to hydrogens!");
+//    }
+//    ++iter;
+//    if (saturation == DoSaturate) {
+//      while ((iter != Leaf->end()) && ((*iter)->getType()->getAtomicNumber() == 1)){ // skip added hydrogen
+//        iter++;
+//      }
+//    }
+//  }
+//}
Index: src/Fragmentation/Exporters/ExportGraph.hpp
===================================================================
--- src/Fragmentation/Exporters/ExportGraph.hpp	(revision 9fd44f7f53b271add9c0ff9fe035e67c63361f09)
+++ src/Fragmentation/Exporters/ExportGraph.hpp	(revision 06804b923c9b82a6dd4f7d932a21781d6e4eae6e)
@@ -77,8 +77,8 @@
   void releaseFragment(SaturatedFragment_ptr &_ptr);
 
-  void prepareMolecule();
-  molecule * StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem);
-  int StoreFragmentFromKeySet_Init(molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList);
-  void CreateInducedSubgraphOfFragment(molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem);
+//  void prepareMolecule();
+//  molecule * StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem);
+//  int StoreFragmentFromKeySet_Init(molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList);
+//  void CreateInducedSubgraphOfFragment(molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem);
 
 protected:
Index: c/Graph/BreadthFirstSearchAdd.cpp
===================================================================
--- src/Graph/BreadthFirstSearchAdd.cpp	(revision 9fd44f7f53b271add9c0ff9fe035e67c63361f09)
+++ 	(revision )
@@ -1,190 +1,0 @@
-/*
- * Project: MoleCuilder
- * Description: creates and alters molecular systems
- * Copyright (C)  2010-2012 University of Bonn. 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/>.
- */
-
-/*
- * BreadthFirstSearchAdd.cpp
- *
- *  Created on: Feb 16, 2011
- *      Author: heber
- */
-
-// include config.h
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "CodePatterns/MemDebug.hpp"
-
-#include "BreadthFirstSearchAdd.hpp"
-
-#include <sstream>
-
-#include "Atom/atom.hpp"
-#include "Bond/bond.hpp"
-#include "CodePatterns/Assert.hpp"
-#include "CodePatterns/Info.hpp"
-#include "CodePatterns/Log.hpp"
-#include "CodePatterns/Verbose.hpp"
-#include "molecule.hpp"
-#include "World.hpp"
-
-
-BreadthFirstSearchAdd::BreadthFirstSearchAdd(atom *&_Root, int _BondOrder, bool _IsAngstroem, const enum HydrogenSaturation _saturation) :
-  BondOrder(_BondOrder),
-  Root(_Root),
-  saturation(_saturation),
-  IsAngstroem(_IsAngstroem)
-{
-  BFSStack.push_front(Root);
-}
-
-
-BreadthFirstSearchAdd::~BreadthFirstSearchAdd()
-{}
-
-void BreadthFirstSearchAdd::Init(atom *&_Root, int _BondOrder)
-{
-  BondOrder = _BondOrder;
-  Root = _Root;
-  BFSStack.clear();
-  BFSStack.push_front(Root);
-}
-
-void BreadthFirstSearchAdd::InitNode(atomId_t atom_id)
-{
-  PredecessorList[atom_id] = NULL;
-  ShortestPathList[atom_id] = -1;
-  if (AddedAtomList.find(atom_id) != AddedAtomList.end()) // mark already present atoms (i.e. Root and maybe others) as visited
-    ColorList[atom_id] = GraphEdge::lightgray;
-  else
-    ColorList[atom_id] = GraphEdge::white;
-}
-
-
-void BreadthFirstSearchAdd::UnvisitedNode(molecule *Mol, atom *&Walker, atom *&OtherAtom, bond::ptr &Binder, bond::ptr &Bond)
-{
-  if (Binder != Bond) // let other atom GraphEdge::white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already GraphEdge::black, thus no problem)
-    ColorList[OtherAtom->getNr()] = GraphEdge::lightgray;
-  PredecessorList[OtherAtom->getNr()] = Walker; // Walker is the predecessor
-  ShortestPathList[OtherAtom->getNr()] = ShortestPathList[Walker->getNr()] + 1;
-  LOG(2, "Coloring OtherAtom " << OtherAtom->getName() << " " << GraphEdge::getColorName(ColorList[OtherAtom->getNr()]) << ", its predecessor is " << Walker->getName() << " and its Shortest Path is " << ShortestPathList[OtherAtom->getNr()] << " egde(s) long.");
-  if ((((ShortestPathList[OtherAtom->getNr()] < BondOrder) && (Binder != Bond)))) { // Check for maximum distance
-    std::stringstream output;
-    if (AddedAtomList[OtherAtom->getNr()] == NULL) { // add if it's not been so far
-      AddedAtomList[OtherAtom->getNr()] = Mol->AddCopyAtom(OtherAtom);
-      output << "Added OtherAtom " << OtherAtom->getName();
-      AddedBondList[Binder] = Mol->CopyBond(AddedAtomList[Walker->getNr()], AddedAtomList[OtherAtom->getNr()], Binder);
-      output << " and bond " << *(AddedBondList[Binder]) << ", ";
-    } else { // this code should actually never come into play (all GraphEdge::white atoms are not yet present in BondMolecule, that's why they are GraphEdge::white in the first place)
-      output << "Not adding OtherAtom " << OtherAtom->getName();
-      if (AddedBondList[Binder] == NULL) {
-        AddedBondList[Binder] = Mol->CopyBond(AddedAtomList[Walker->getNr()], AddedAtomList[OtherAtom->getNr()], Binder);
-        output << ", added Bond " << *(AddedBondList[Binder]);
-      } else
-        output << ", not added Bond ";
-    }
-    output << ", putting OtherAtom into queue.";
-    LOG(0, output.str());
-    BFSStack.push_front(OtherAtom);
-  } else { // out of bond order, then replace
-    if ((AddedAtomList[OtherAtom->getNr()] == NULL) && (Binder->Cyclic))
-      ColorList[OtherAtom->getNr()] = GraphEdge::white; // unmark if it has not been queued/added, to make it available via its other bonds (cyclic)
-    {
-      std::stringstream output;
-      if (Binder == Bond)
-        output << "Not Queueing, is the Root bond";
-      else if (ShortestPathList[OtherAtom->getNr()] >= BondOrder)
-        output << "Not Queueing, is out of Bond Count of " << BondOrder;
-      if (!Binder->Cyclic)
-        output << ", is not part of a cyclic bond, saturating bond with Hydrogen.";
-      LOG(3, output.str());
-    }
-    if (AddedBondList[Binder] == NULL) {
-      if ((AddedAtomList[OtherAtom->getNr()] != NULL)) { // .. whether we add or saturate
-        AddedBondList[Binder] = Mol->CopyBond(AddedAtomList[Walker->getNr()], AddedAtomList[OtherAtom->getNr()], Binder);
-      } else {
-        if (saturation == DoSaturate)
-          if (!Mol->AddHydrogenReplacementAtom(Binder, AddedAtomList[Walker->getNr()], Walker, OtherAtom, IsAngstroem))
-            exit(1);
-      }
-    }
-  }
-}
-
-
-void BreadthFirstSearchAdd::VisitedNode(molecule *Mol, atom *&Walker, atom *&OtherAtom, bond::ptr &Binder, bond::ptr &Bond)
-{
-  LOG(3, "Not Adding, has already been visited.");
-  // This has to be a cyclic bond, check whether it's present ...
-  if (AddedBondList[Binder] == NULL) {
-    if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->getNr()] + 1) < BondOrder))) {
-      AddedBondList[Binder] = Mol->CopyBond(AddedAtomList[Walker->getNr()], AddedAtomList[OtherAtom->getNr()], Binder);
-    } else { // if it's root bond it has to broken (otherwise we would not create the fragments)
-      if (saturation == DoSaturate)
-        if(!Mol->AddHydrogenReplacementAtom(Binder, AddedAtomList[Walker->getNr()], Walker, OtherAtom, IsAngstroem))
-          exit(1);
-    }
-  }
-}
-
-
-void BreadthFirstSearchAdd::operator()(molecule *Mol, atom *_Root, bond::ptr Bond, int _BondOrder)
-{
-  Info FunctionInfo("BreadthFirstSearchAdd");
-  atom *Walker = NULL, *OtherAtom = NULL;
-  bond::ptr Binder;
-
-  // add Root if not done yet
-  if (AddedAtomList[_Root->getNr()] == NULL) // add Root if not yet present
-    AddedAtomList[_Root->getNr()] = Mol->AddCopyAtom(_Root);
-
-  Init(_Root, _BondOrder);
-
-  // and go on ... Queue always contains all GraphEdge::lightgray vertices
-  while (!BFSStack.empty()) {
-    // we have to pop the oldest atom from stack. This keeps the atoms on the stack always of the same ShortestPath distance.
-    // e.g. if current atom is 2, push to end of stack are of length 3, but first all of length 2 would be popped. They again
-    // append length of 3 (their neighbours). Thus on stack we have always atoms of a certain length n at bottom of stack and
-    // followed by n+1 till top of stack.
-    Walker = BFSStack.front(); // pop oldest added
-    BFSStack.pop_front();
-    const BondList& ListOfBonds = Walker->getListOfBonds();
-    LOG(1, "Current Walker is: " << Walker->getName() << ", and has " << ListOfBonds.size() << " bonds.");
-    for (BondList::const_iterator Runner = ListOfBonds.begin();
-        Runner != ListOfBonds.end();
-        ++Runner) {
-      if ((*Runner) != NULL) { // don't look at bond equal NULL
-        Binder = (*Runner);
-        OtherAtom = (*Runner)->GetOtherAtom(Walker);
-        LOG(2, "Current OtherAtom is: " << OtherAtom->getName() << " for bond " << *(*Runner) << ".");
-        if (ColorList[OtherAtom->getNr()] == GraphEdge::white) {
-          UnvisitedNode(Mol, Walker, OtherAtom, Binder, Bond);
-        } else {
-          VisitedNode(Mol, Walker, OtherAtom, Binder, Bond);
-        }
-      }
-    }
-    ColorList[Walker->getNr()] = GraphEdge::black;
-    LOG(1, "Coloring Walker " << Walker->getName() << " GraphEdge::black.");
-  }
-}
-
Index: c/Graph/BreadthFirstSearchAdd.hpp
===================================================================
--- src/Graph/BreadthFirstSearchAdd.hpp	(revision 9fd44f7f53b271add9c0ff9fe035e67c63361f09)
+++ 	(revision )
@@ -1,82 +1,0 @@
-/*
- * BreadthFirstSearchAdd.hpp
- *
- *  Created on: Feb 16, 2011
- *      Author: heber
- */
-
-#ifndef BREADTHFIRSTSEARCHADD_HPP_
-#define BREADTHFIRSTSEARCHADD_HPP_
-
-// include config.h
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <deque>
-#include <map>
-
-#include "Bond/bond.hpp"
-#include "Bond/GraphEdge.hpp"
-#include "Fragmentation/HydrogenSaturation_enum.hpp"
-#include "Helpers/defs.hpp"
-#include "types.hpp"
-
-class atom;
-class molecule;
-
-class BreadthFirstSearchAdd
-{
-public:
-  BreadthFirstSearchAdd(atom *&_Root, int _BondOrder, bool _IsAngstroem, const enum HydrogenSaturation _saturation);
-  ~BreadthFirstSearchAdd();
-
-  /** Reinitializes the structure for a new BFS run.
-   *
-   * @param _Root Root atom to start BFS from
-   * @param _BondOrder horizon in bond counts when to stop BFS
-   */
-  void Init(atom *&_Root, int _BondOrder);
-
-  /** Adds atoms up to \a BondCount distance from \a *Root and notes them down in \a AddedAtomList.
-   * Gray vertices are always enqueued in an std::deque<atom *> FIFO queue, the rest is usual BFS with adding vertices found was
-   * white and putting into queue.
-   * \param *Mol Molecule class to add atoms to
-   * \param *Root root vertex for BFS
-   * \param *Bond bond not to look beyond
-   * \param BondOrder maximum distance for vertices to add
-   * \param IsAngstroem lengths are in angstroem or bohrradii
-   */
-  void operator()(molecule *Mol, atom *_Root, bond::ptr Bond, int _BondOrder);
-
-private:
-  void UnvisitedNode(molecule *Mol, atom *&Walker, atom *&OtherAtom, bond::ptr &Binder, bond::ptr &Bond);
-  void VisitedNode(molecule *Mol, atom *&Walker, atom *&OtherAtom, bond::ptr &Binder, bond::ptr &Bond);
-
-  /** initialise vertex as white with no predecessor, empty queue, color Root lightgray.
-   *
-   */
-  void InitNode(atomId_t atom_id);
-
-  std::map<atomId_t, atom *> PredecessorList;
-  std::map<atomId_t, int > ShortestPathList;
-  std::map<atomId_t, enum GraphEdge::Shading> ColorList;
-  std::deque<atom *> BFSStack;
-  int BondOrder;
-  atom *Root;
-
-  bool BackStepping;
-  int CurrentGraphNr;
-  int ComponentNr;
-  std::map<atomId_t, atom *> AddedAtomList;
-  std::map<bond::ptr , bond::ptr > AddedBondList; //!< maps from father bond to son
-
-  //!> whether to treat hydrogen special or not
-  const enum HydrogenSaturation saturation;
-
-  //!> is angstroem our unit of length
-  bool IsAngstroem;
-};
-
-
-#endif /* BREADTHFIRSTSEARCHADD_HPP_ */
Index: src/Graph/Makefile.am
===================================================================
--- src/Graph/Makefile.am	(revision 9fd44f7f53b271add9c0ff9fe035e67c63361f09)
+++ src/Graph/Makefile.am	(revision 06804b923c9b82a6dd4f7d932a21781d6e4eae6e)
@@ -4,5 +4,4 @@
 GRAPHSOURCE = \
 	Graph/BondGraph.cpp \
-	Graph/BreadthFirstSearchAdd.cpp \
 	Graph/BuildInducedSubgraph.cpp \
 	Graph/AdjacencyList.cpp \
@@ -13,5 +12,4 @@
 GRAPHHEADER = \
 	Graph/BondGraph.hpp \
-	Graph/BreadthFirstSearchAdd.hpp \
 	Graph/BuildInducedSubgraph.hpp \
 	Graph/AdjacencyList.hpp \
Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision 9fd44f7f53b271add9c0ff9fe035e67c63361f09)
+++ src/molecule.cpp	(revision 06804b923c9b82a6dd4f7d932a21781d6e4eae6e)
@@ -339,251 +339,251 @@
  * \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;
-};
+//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.
Index: src/molecule.hpp
===================================================================
--- src/molecule.hpp	(revision 9fd44f7f53b271add9c0ff9fe035e67c63361f09)
+++ src/molecule.hpp	(revision 06804b923c9b82a6dd4f7d932a21781d6e4eae6e)
@@ -261,5 +261,5 @@
   /// Add/remove atoms to/from molecule.
   atom * AddCopyAtom(atom *pointer);
-  bool AddHydrogenReplacementAtom(bond::ptr Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem);
+//  bool AddHydrogenReplacementAtom(bond::ptr Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem);
   bond::ptr AddBond(atom *first, atom *second, int degree = 1);
   bool hasBondStructure() const;
