Changes in src/molecule_graph.cpp [f66195:99593f]

File:

• src/molecule_graph.cpp (modified) (27 diffs)

src/molecule_graph.cpp

@@ -8 +8 @@
 #include "atom.hpp"
 #include "bond.hpp"
+#include "bondgraph.hpp"
 #include "config.hpp"
 #include "element.hpp"
 #include "helpers.hpp"
+#include "linkedcell.hpp"
 #include "lists.hpp"
 #include "memoryallocator.hpp"
 #include "molecule.hpp"
 
+struct BFSAccounting
+{
+  atom **PredecessorList;
+  int *ShortestPathList;
+  enum Shading *ColorList;
+  class StackClass<atom *> *BFSStack;
+  class StackClass<atom *> *TouchedStack;
+  int AtomCount;
+  int BondOrder;
+  atom *Root;
+  bool BackStepping;
+  int CurrentGraphNr;
+  int ComponentNr;
+};
+
+/** Accounting data for Depth First Search.
+ */
+struct DFSAccounting
+{
+  class StackClass<atom *> *AtomStack;
+  class StackClass<bond *> *BackEdgeStack;
+  int CurrentGraphNr;
+  int ComponentNumber;
+  atom *Root;
+  bool BackStepping;
+};
+
 /************************************* Functions for class molecule *********************************/
-
 
 /** Creates an adjacency list of the molecule.
  * We obtain an outside file with the indices of atoms which are bondmembers.
  */
-void molecule::CreateAdjacencyList2(ofstream *out, ifstream *input)
+void molecule::CreateAdjacencyListFromDbondFile(ofstream *out, ifstream *input)
 {
 
   // 1 We will parse bonds out of the dbond file created by tremolo.
-      int atom1, atom2, temp;
-      atom *Walker, *OtherWalker;
-
-          if (!input)
-          {
-            cout << Verbose(1) << "Opening silica failed \n";
-          };
-
-      *input >> ws >> atom1;
-      *input >> ws >> atom2;
-          cout << Verbose(1) << "Scanning file\n";
-          while (!input->eof()) // Check whether we read everything already
-          {
-        *input >> ws >> atom1;
-        *input >> ws >> atom2;
-            if(atom2<atom1) //Sort indices of atoms in order
-            {
-              temp=atom1;
-              atom1=atom2;
-              atom2=temp;
-            };
-
-            Walker=start;
-            while(Walker-> nr != atom1) // Find atom corresponding to first index
-            {
-              Walker = Walker->next;
-            };
-            OtherWalker = Walker->next;
-            while(OtherWalker->nr != atom2) // Find atom corresponding to second index
-            {
-              OtherWalker= OtherWalker->next;
-            };
-            AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
-
-          }
-
-          CreateListOfBondsPerAtom(out);
-
-};
-
+  int atom1, atom2;
+  atom *Walker, *OtherWalker;
+
+  if (!input) {
+    cout << Verbose(1) << "Opening silica failed \n";
+  };
+
+  *input >> ws >> atom1;
+  *input >> ws >> atom2;
+  cout << Verbose(1) << "Scanning file\n";
+  while (!input->eof()) // Check whether we read everything already
+  {
+    *input >> ws >> atom1;
+    *input >> ws >> atom2;
+
+    if (atom2 < atom1) //Sort indices of atoms in order
+      flip(atom1, atom2);
+    Walker = FindAtom(atom1);
+    OtherWalker = FindAtom(atom2);
+    AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
+  }
+}
+;
 
 /** Creates an adjacency list of the molecule.
@@ -77 +88 @@
  *  -# put each atom into its corresponding cell
  *  -# go through every cell, check the atoms therein against all possible bond partners in the 27 adjacent cells, add bond if true
- *  -# create the list of bonds via CreateListOfBondsPerAtom()
  *  -# correct the bond degree iteratively (single->double->triple bond)
  *  -# finally print the bond list to \a *out if desired
@@ -83 +93 @@
  * \param bonddistance length of linked cells (i.e. maximum minimal length checked)
  * \param IsAngstroem whether coordinate system is gauged to Angstroem or Bohr radii
- */
-void molecule::CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem)
-{
-
-  atom *Walker = NULL, *OtherWalker = NULL, *Candidate = NULL;
-  int No, NoBonds, CandidateBondNo;
-  int NumberCells, divisor[NDIM], n[NDIM], N[NDIM], index, Index, j;
-  molecule **CellList;
-  double distance, MinDistance, MaxDistance;
-  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
-  Vector x;
-  int FalseBondDegree = 0;
+ * \param *minmaxdistance function to give upper and lower bound on whether particle is bonded to some other
+ * \param *BG BondGraph with the member function above or NULL, if just standard covalent should be used.
+ */
+void molecule::CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem, void (BondGraph::*minmaxdistance)(BondedParticle * const , BondedParticle * const , double &, double &, bool), BondGraph *BG)
+{
+  atom *Walker = NULL;
+  atom *OtherWalker = NULL;
+  atom **AtomMap = NULL;
+  int n[NDIM];
+  double MinDistance, MaxDistance;
+  LinkedCell *LC = NULL;
+  bool free_BG = false;
+
+  if (BG == NULL) {
+    BG = new BondGraph(IsAngstroem);
+    free_BG = true;
+  }
 
   BondDistance = bonddistance; // * ((IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem);
   *out << Verbose(0) << "Begin of CreateAdjacencyList." << endl;
   // remove every bond from the list
-  if ((first->next != last) && (last->previous != first)) {  // there are bonds present
-    cleanup(first,last);
-  }
+  bond *Binder = NULL;
+  while (last->previous != first) {
+    Binder = last->previous;
+    Binder->leftatom->UnregisterBond(Binder);
+    Binder->rightatom->UnregisterBond(Binder);
+    removewithoutcheck(Binder);
+  }
+  BondCount = 0;
 
   // count atoms in molecule = dimension of matrix (also give each unique name and continuous numbering)
   CountAtoms(out);
-  *out << Verbose(1) << "AtomCount " << AtomCount << "." << endl;
-
-  if (AtomCount != 0) {
-    // 1. find divisor for each axis, such that a sphere with radius of at least bonddistance can be placed into each cell
-    j=-1;
-    for (int i=0;i<NDIM;i++) {
-      j += i+1;
-      divisor[i] = (int)floor(cell_size[j]/bonddistance); // take smaller value such that size of linked cell is at least bonddistance
-      //*out << Verbose(1) << "divisor[" << i << "]  = " << divisor[i] << "." << endl;
-    }
-    // 2a. allocate memory for the cell list
-    NumberCells = divisor[0]*divisor[1]*divisor[2];
-    *out << Verbose(1) << "Allocating " << NumberCells << " cells." << endl;
-    CellList = Malloc<molecule*>(NumberCells, "molecule::CreateAdjacencyList - ** CellList");
-    for (int i=NumberCells;i--;)
-      CellList[i] = NULL;
-
-    // 2b. put all atoms into its corresponding list
+  *out << Verbose(1) << "AtomCount " << AtomCount << " and bonddistance is " << bonddistance << "." << endl;
+
+  if ((AtomCount > 1) && (bonddistance > 1.)) {
+    *out << Verbose(2) << "Creating Linked Cell structure ... " << endl;
+    LC = new LinkedCell(this, bonddistance);
+
+    // create a list to map Tesselpoint::nr to atom *
+    *out << Verbose(2) << "Creating TesselPoint to atom map ... " << endl;
+    AtomMap = Calloc<atom *> (AtomCount, "molecule::CreateAdjacencyList - **AtomCount");
     Walker = start;
-    while(Walker->next != end) {
+    while (Walker->next != end) {
       Walker = Walker->next;
-      //*out << Verbose(1) << "Current atom is " << *Walker << " with coordinates ";
-      //Walker->x.Output(out);
-      //*out << "." << endl;
-      // compute the cell by the atom's coordinates
-      j=-1;
-      for (int i=0;i<NDIM;i++) {
-        j += i+1;
-        x.CopyVector(&(Walker->x));
-        x.KeepPeriodic(out, matrix);
-        n[i] = (int)floor(x.x[i]/cell_size[j]*(double)divisor[i]);
-      }
-      index = n[2] + (n[1] + n[0] * divisor[1]) * divisor[2];
-      //*out << Verbose(1) << "Atom " << *Walker << " goes into cell number [" << n[0] << "," << n[1] << "," << n[2] << "] = " << index << "." << endl;
-      // add copy atom to this cell
-      if (CellList[index] == NULL)  // allocate molecule if not done
-        CellList[index] = new molecule(elemente);
-      OtherWalker = CellList[index]->AddCopyAtom(Walker); // add a copy of walker to this atom, father will be walker for later reference
-      //*out << Verbose(1) << "Copy Atom is " << *OtherWalker << "." << endl;
-    }
-    //for (int i=0;i<NumberCells;i++)
-      //*out << Verbose(1) << "Cell number " << i << ": " << CellList[i] << "." << endl;
-
+      AtomMap[Walker->nr] = Walker;
+    }
 
     // 3a. go through every cell
-    for (N[0]=divisor[0];N[0]--;)
-      for (N[1]=divisor[1];N[1]--;)
-        for (N[2]=divisor[2];N[2]--;) {
-          Index = N[2] + (N[1] + N[0] * divisor[1]) * divisor[2];
-          if (CellList[Index] != NULL) { // if there atoms in this cell
-            //*out << Verbose(1) << "Current cell is " << Index << "." << endl;
-            // 3b. for every atom therein
-            Walker = CellList[Index]->start;
-            while (Walker->next != CellList[Index]->end) {  // go through every atom
-              Walker = Walker->next;
+    *out << Verbose(2) << "Celling ... " << endl;
+    for (LC->n[0] = 0; LC->n[0] < LC->N[0]; LC->n[0]++)
+      for (LC->n[1] = 0; LC->n[1] < LC->N[1]; LC->n[1]++)
+        for (LC->n[2] = 0; LC->n[2] < LC->N[2]; LC->n[2]++) {
+          const LinkedNodes *List = LC->GetCurrentCell();
+          //*out << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
+          if (List != NULL) {
+            for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+              Walker = AtomMap[(*Runner)->nr];
               //*out << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
               // 3c. check for possible bond between each atom in this and every one in the 27 cells
-              for (n[0]=-1;n[0]<=1;n[0]++)
-                for (n[1]=-1;n[1]<=1;n[1]++)
-                  for (n[2]=-1;n[2]<=1;n[2]++) {
-                     // compute the index of this comparison cell and make it periodic
-                    index = ((N[2]+n[2]+divisor[2])%divisor[2]) + (((N[1]+n[1]+divisor[1])%divisor[1]) + ((N[0]+n[0]+divisor[0])%divisor[0]) * divisor[1]) * divisor[2];
-                    //*out << Verbose(1) << "Number of comparison cell is " << index << "." << endl;
-                    if (CellList[index] != NULL) {  // if there are any atoms in this cell
-                      OtherWalker = CellList[index]->start;
-                      while(OtherWalker->next != CellList[index]->end) {  // go through every atom in this cell
-                        OtherWalker = OtherWalker->next;
-                        //*out << Verbose(0) << "Current comparison atom is " << *OtherWalker << "." << endl;
-                        /// \todo periodic check is missing here!
-                        //*out << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
-                        MinDistance = OtherWalker->type->CovalentRadius + Walker->type->CovalentRadius;
-                        MinDistance *= (IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem;
-                        MaxDistance = MinDistance + BONDTHRESHOLD;
-                        MinDistance -= BONDTHRESHOLD;
-                        distance = OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size);
-                        if ((OtherWalker->father->nr > Walker->father->nr) && (distance <= MaxDistance*MaxDistance) && (distance >= MinDistance*MinDistance)) { // create bond if distance is smaller
-                          //*out << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
-                          AddBond(Walker->father, OtherWalker->father, 1);  // also increases molecule::BondCount
-                        } else {
-                          //*out << Verbose(1) << "Not Adding: Wrong label order or distance too great." << endl;
+              for (n[0] = -1; n[0] <= 1; n[0]++)
+                for (n[1] = -1; n[1] <= 1; n[1]++)
+                  for (n[2] = -1; n[2] <= 1; n[2]++) {
+                    const LinkedNodes *OtherList = LC->GetRelativeToCurrentCell(n);
+                    //*out << Verbose(2) << "Current relative cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
+                    if (OtherList != NULL) {
+                      for (LinkedNodes::const_iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
+                        if ((*OtherRunner)->nr > Walker->nr) {
+                          OtherWalker = AtomMap[(*OtherRunner)->nr];
+                          //*out << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
+                          (BG->*minmaxdistance)(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
+                          const double distance = OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size);
+                          const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
+                          if ((OtherWalker->father->nr > Walker->father->nr) && (status)) { // create bond if distance is smaller
+                            //*out << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
+                            AddBond(Walker->father, OtherWalker->father, 1); // also increases molecule::BondCount
+                          } else {
+                            //*out << Verbose(1) << "Not Adding: Wrong label order or distance too great." << endl;
+                          }
                         }
                       }
@@ -192 +178 @@
           }
         }
-
-
-
-    // 4. free the cell again
-    for (int i=NumberCells;i--;)
-      if (CellList[i] != NULL) {
-        delete(CellList[i]);
-      }
-    Free(&CellList);
-
-    // create the adjacency list per atom
-    CreateListOfBondsPerAtom(out);
-
-    // correct Bond degree of each bond by checking both bond partners for a mismatch between valence and current sum of bond degrees,
-    // iteratively increase the one first where the other bond partner has the fewest number of bonds (i.e. in general bonds oxygene
-    // preferred over carbon bonds). Beforehand, we had picked the first mismatching partner, which lead to oxygenes with single instead of
-    // double bonds as was expected.
-    if (BondCount != 0) {
-      NoCyclicBonds = 0;
-      *out << Verbose(1) << "Correcting Bond degree of each bond ... ";
-      do {
-        No = 0; // No acts as breakup flag (if 1 we still continue)
-        Walker = start;
-        while (Walker->next != end) { // go through every atom
-          Walker = Walker->next;
-          // count valence of first partner
-          NoBonds = 0;
-          for(j=0;j<NumberOfBondsPerAtom[Walker->nr];j++)
-            NoBonds += ListOfBondsPerAtom[Walker->nr][j]->BondDegree;
-          *out << Verbose(3) << "Walker " << *Walker << ": " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-          if ((int)(Walker->type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
-            Candidate = NULL;
-            CandidateBondNo = -1;
-            for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) { // go through each of its bond partners
-              OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
-              // count valence of second partner
-              NoBonds = 0;
-              for(j=0;j<NumberOfBondsPerAtom[OtherWalker->nr];j++)
-                NoBonds += ListOfBondsPerAtom[OtherWalker->nr][j]->BondDegree;
-              *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-              if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
-                if ((Candidate == NULL) || (NumberOfBondsPerAtom[Candidate->nr] > NumberOfBondsPerAtom[OtherWalker->nr])) { // pick the one with fewer number of bonds first
-                  Candidate = OtherWalker;
-                  CandidateBondNo = i;
-                  *out << Verbose(3) << "New candidate is " << *Candidate << "." << endl;
-                }
-              }
-            }
-            if ((Candidate != NULL) && (CandidateBondNo != -1)) {
-              ListOfBondsPerAtom[Walker->nr][CandidateBondNo]->BondDegree++;
-              *out << Verbose(2) << "Increased bond degree for bond " << *ListOfBondsPerAtom[Walker->nr][CandidateBondNo] << "." << endl;
-            } else
-              *out << Verbose(2) << "Could not find correct degree for atom " << *Walker << "." << endl;
-              FalseBondDegree++;
-          }
-        }
-      } while (No);
-    *out << " done." << endl;
-    } else
-      *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
-    *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
+    Free(&AtomMap);
+    delete (LC);
+    *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << BondDistance << "." << endl;
+
+    // correct bond degree by comparing valence and bond degree
+    *out << Verbose(2) << "Correcting bond degree ... " << endl;
+    CorrectBondDegree(out);
 
     // output bonds for debugging (if bond chain list was correctly installed)
-    *out << Verbose(1) << endl << "From contents of bond chain list:";
-    bond *Binder = first;
-    while(Binder->next != last) {
-      Binder = Binder->next;
-      *out << *Binder << "\t" << endl;
-    }
-    *out << endl;
+    ActOnAllAtoms(&atom::OutputBondOfAtom, out);
   } else
     *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
   *out << Verbose(0) << "End of CreateAdjacencyList." << endl;
-  Free(&matrix);
-
-};
+  if (free_BG)
+    delete(BG);
+}
+;
+
+/** Prints a list of all bonds to \a *out.
+ * \param output stream
+ */
+void molecule::OutputBondsList(ofstream *out) const
+{
+  *out << Verbose(1) << endl << "From contents of bond chain list:";
+  bond *Binder = first;
+  while (Binder->next != last) {
+    Binder = Binder->next;
+    *out << *Binder << "\t" << endl;
+  }
+  *out << endl;
+}
+;
+
+/** correct bond degree by comparing valence and bond degree.
+ * correct Bond degree of each bond by checking both bond partners for a mismatch between valence and current sum of bond degrees,
+ * iteratively increase the one first where the other bond partner has the fewest number of bonds (i.e. in general bonds oxygene
+ * preferred over carbon bonds). Beforehand, we had picked the first mismatching partner, which lead to oxygenes with single instead of
+ * double bonds as was expected.
+ * \param *out output stream for debugging
+ * \return number of bonds that could not be corrected
+ */
+int molecule::CorrectBondDegree(ofstream *out) const
+{
+  int No = 0, OldNo = -1;
+
+  if (BondCount != 0) {
+    *out << Verbose(1) << "Correcting Bond degree of each bond ... " << endl;
+    do {
+      OldNo = No;
+      No = SumPerAtom(&atom::CorrectBondDegree, out);
+    } while (OldNo != No);
+    *out << " done." << endl;
+  } else {
+    *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
+  }
+  *out << No << " bonds could not be corrected." << endl;
+
+  return (No);
+}
+;
 
 /** Counts all cyclic bonds and returns their number.
@@ -276 +246 @@
 int molecule::CountCyclicBonds(ofstream *out)
 {
-  int No = 0;
+  NoCyclicBonds = 0;
   int *MinimumRingSize = NULL;
   MoleculeLeafClass *Subgraphs = NULL;
@@ -286 +256 @@
     while (Subgraphs->next != NULL) {
       Subgraphs = Subgraphs->next;
-      delete(Subgraphs->previous);
-    }
-    delete(Subgraphs);
-    delete[](MinimumRingSize);
-  }
-  while(Binder->next != last) {
+      delete (Subgraphs->previous);
+    }
+    delete (Subgraphs);
+    delete[] (MinimumRingSize);
+  }
+  while (Binder->next != last) {
     Binder = Binder->next;
     if (Binder->Cyclic)
-      No++;
-  }
-  delete(BackEdgeStack);
-  return No;
-};
+      NoCyclicBonds++;
+  }
+  delete (BackEdgeStack);
+  return NoCyclicBonds;
+}
+;
+
 /** Returns Shading as a char string.
  * \param color the Shading
  * \return string of the flag
  */
-string molecule::GetColor(enum Shading color)
-{
-  switch(color) {
+string molecule::GetColor(enum Shading color) const
+{
+  switch (color) {
     case white:
       return "white";
@@ -322 +294 @@
       break;
   };
-};
+}
+;
+
+/** Sets atom::GraphNr and atom::LowpointNr to BFSAccounting::CurrentGraphNr.
+ * \param *out output stream for debugging
+ * \param *Walker current node
+ * \param &BFS structure with accounting data for BFS
+ */
+void DepthFirstSearchAnalysis_SetWalkersGraphNr(ofstream *out, atom *&Walker, struct DFSAccounting &DFS)
+{
+  if (!DFS.BackStepping) { // if we don't just return from (8)
+    Walker->GraphNr = DFS.CurrentGraphNr;
+    Walker->LowpointNr = DFS.CurrentGraphNr;
+    *out << Verbose(1) << "Setting Walker[" << Walker->Name << "]'s number to " << Walker->GraphNr << " with Lowpoint " << Walker->LowpointNr << "." << endl;
+    DFS.AtomStack->Push(Walker);
+    DFS.CurrentGraphNr++;
+  }
+}
+;
+
+/** During DFS goes along unvisited bond and touches other atom.
+ * Sets bond::type, if
+ *  -# BackEdge: set atom::LowpointNr and push on \a BackEdgeStack
+ *  -# TreeEgde: set atom::Ancestor and continue with Walker along this edge
+ * Continue until molecule::FindNextUnused() finds no more unused bonds.
+ * \param *out output stream for debugging
+ * \param *mol molecule with atoms and finding unused bonds
+ * \param *&Binder current edge
+ * \param &DFS DFS accounting data
+ */
+void DepthFirstSearchAnalysis_ProbeAlongUnusedBond(ofstream *out, const molecule * const mol, atom *&Walker, bond *&Binder, struct DFSAccounting &DFS)
+{
+  atom *OtherAtom = NULL;
+
+  do { // (3) if Walker has no unused egdes, go to (5)
+    DFS.BackStepping = false; // reset backstepping flag for (8)
+    if (Binder == NULL) // if we don't just return from (11), Binder is already set to next unused
+      Binder = mol->FindNextUnused(Walker);
+    if (Binder == NULL)
+      break;
+    *out << Verbose(2) << "Current Unused Bond is " << *Binder << "." << endl;
+    // (4) Mark Binder used, ...
+    Binder->MarkUsed(black);
+    OtherAtom = Binder->GetOtherAtom(Walker);
+    *out << Verbose(2) << "(4) OtherAtom is " << OtherAtom->Name << "." << endl;
+    if (OtherAtom->GraphNr != -1) {
+      // (4a) ... if "other" atom has been visited (GraphNr != 0), set lowpoint to minimum of both, go to (3)
+      Binder->Type = BackEdge;
+      DFS.BackEdgeStack->Push(Binder);
+      Walker->LowpointNr = (Walker->LowpointNr < OtherAtom->GraphNr) ? Walker->LowpointNr : OtherAtom->GraphNr;
+      *out << Verbose(3) << "(4a) Visited: Setting Lowpoint of Walker[" << Walker->Name << "] to " << Walker->LowpointNr << "." << endl;
+    } else {
+      // (4b) ... otherwise set OtherAtom as Ancestor of Walker and Walker as OtherAtom, go to (2)
+      Binder->Type = TreeEdge;
+      OtherAtom->Ancestor = Walker;
+      Walker = OtherAtom;
+      *out << Verbose(3) << "(4b) Not Visited: OtherAtom[" << OtherAtom->Name << "]'s Ancestor is now " << OtherAtom->Ancestor->Name << ", Walker is OtherAtom " << OtherAtom->Name << "." << endl;
+      break;
+    }
+    Binder = NULL;
+  } while (1); // (3)
+}
+;
+
+/** Checks whether we have a new component.
+ * if atom::LowpointNr of \a *&Walker is greater than atom::GraphNr of its atom::Ancestor, we have a new component.
+ * Meaning that if we touch upon a node who suddenly has a smaller atom::LowpointNr than its ancestor, then we
+ * have a found a new branch in the graph tree.
+ * \param *out output stream for debugging
+ * \param *mol molecule with atoms and finding unused bonds
+ * \param *&Walker current node
+ * \param &DFS DFS accounting data
+ */
+void DepthFirstSearchAnalysis_CheckForaNewComponent(ofstream *out, const molecule * const mol, atom *&Walker, struct DFSAccounting &DFS, MoleculeLeafClass *&LeafWalker)
+{
+  atom *OtherAtom = NULL;
+
+  // (5) if Ancestor of Walker is ...
+  *out << Verbose(1) << "(5) Number of Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "] is " << Walker->Ancestor->GraphNr << "." << endl;
+
+  if (Walker->Ancestor->GraphNr != DFS.Root->GraphNr) {
+    // (6)  (Ancestor of Walker is not Root)
+    if (Walker->LowpointNr < Walker->Ancestor->GraphNr) {
+      // (6a) set Ancestor's Lowpoint number to minimum of of its Ancestor and itself, go to Step(8)
+      Walker->Ancestor->LowpointNr = (Walker->Ancestor->LowpointNr < Walker->LowpointNr) ? Walker->Ancestor->LowpointNr : Walker->LowpointNr;
+      *out << Verbose(2) << "(6) Setting Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "]'s Lowpoint to " << Walker->Ancestor->LowpointNr << "." << endl;
+    } else {
+      // (7) (Ancestor of Walker is a separating vertex, remove all from stack till Walker (including), these and Ancestor form a component
+      Walker->Ancestor->SeparationVertex = true;
+      *out << Verbose(2) << "(7) Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "]'s is a separating vertex, creating component." << endl;
+      mol->SetNextComponentNumber(Walker->Ancestor, DFS.ComponentNumber);
+      *out << Verbose(3) << "(7) Walker[" << Walker->Name << "]'s Ancestor's Compont is " << DFS.ComponentNumber << "." << endl;
+      mol->SetNextComponentNumber(Walker, DFS.ComponentNumber);
+      *out << Verbose(3) << "(7) Walker[" << Walker->Name << "]'s Compont is " << DFS.ComponentNumber << "." << endl;
+      do {
+        OtherAtom = DFS.AtomStack->PopLast();
+        LeafWalker->Leaf->AddCopyAtom(OtherAtom);
+        mol->SetNextComponentNumber(OtherAtom, DFS.ComponentNumber);
+        *out << Verbose(3) << "(7) Other[" << OtherAtom->Name << "]'s Compont is " << DFS.ComponentNumber << "." << endl;
+      } while (OtherAtom != Walker);
+      DFS.ComponentNumber++;
+    }
+    // (8) Walker becomes its Ancestor, go to (3)
+    *out << Verbose(2) << "(8) Walker[" << Walker->Name << "] is now its Ancestor " << Walker->Ancestor->Name << ", backstepping. " << endl;
+    Walker = Walker->Ancestor;
+    DFS.BackStepping = true;
+  }
+}
+;
+
+/** Cleans the root stack when we have found a component.
+ * If we are not DFSAccounting::BackStepping, then we clear the root stack by putting everything into a
+ * component down till we meet DFSAccounting::Root.
+ * \param *out output stream for debugging
+ * \param *mol molecule with atoms and finding unused bonds
+ * \param *&Walker current node
+ * \param *&Binder current edge
+ * \param &DFS DFS accounting data
+ */
+void DepthFirstSearchAnalysis_CleanRootStackDownTillWalker(ofstream *out, const molecule * const mol, atom *&Walker, bond *&Binder, struct DFSAccounting &DFS, MoleculeLeafClass *&LeafWalker)
+{
+  atom *OtherAtom = NULL;
+
+  if (!DFS.BackStepping) { // coming from (8) want to go to (3)
+    // (9) remove all from stack till Walker (including), these and Root form a component
+    //DFS.AtomStack->Output(out);
+    mol->SetNextComponentNumber(DFS.Root, DFS.ComponentNumber);
+    *out << Verbose(3) << "(9) Root[" << DFS.Root->Name << "]'s Component is " << DFS.ComponentNumber << "." << endl;
+    mol->SetNextComponentNumber(Walker, DFS.ComponentNumber);
+    *out << Verbose(3) << "(9) Walker[" << Walker->Name << "]'s Component is " << DFS.ComponentNumber << "." << endl;
+    do {
+      OtherAtom = DFS.AtomStack->PopLast();
+      LeafWalker->Leaf->AddCopyAtom(OtherAtom);
+      mol->SetNextComponentNumber(OtherAtom, DFS.ComponentNumber);
+      *out << Verbose(3) << "(7) Other[" << OtherAtom->Name << "]'s Compont is " << DFS.ComponentNumber << "." << endl;
+    } while (OtherAtom != Walker);
+    DFS.ComponentNumber++;
+
+    // (11) Root is separation vertex,  set Walker to Root and go to (4)
+    Walker = DFS.Root;
+    Binder = mol->FindNextUnused(Walker);
+    *out << Verbose(1) << "(10) Walker is Root[" << DFS.Root->Name << "], next Unused Bond is " << Binder << "." << endl;
+    if (Binder != NULL) { // Root is separation vertex
+      *out << Verbose(1) << "(11) Root is a separation vertex." << endl;
+      Walker->SeparationVertex = true;
+    }
+  }
+}
+;
+
+/** Initializes DFSAccounting structure.
+ * \param *out output stream for debugging
+ * \param &DFS accounting structure to allocate
+ * \param *mol molecule with AtomCount, BondCount and all atoms
+ */
+void DepthFirstSearchAnalysis_Init(ofstream *out, struct DFSAccounting &DFS, const molecule * const mol)
+{
+  DFS.AtomStack = new StackClass<atom *> (mol->AtomCount);
+  DFS.CurrentGraphNr = 0;
+  DFS.ComponentNumber = 0;
+  DFS.BackStepping = false;
+  mol->ResetAllBondsToUnused();
+  mol->SetAtomValueToValue(-1, &atom::GraphNr);
+  mol->ActOnAllAtoms(&atom::InitComponentNr);
+  DFS.BackEdgeStack->ClearStack();
+}
+;
+
+/** Free's DFSAccounting structure.
+ * \param *out output stream for debugging
+ * \param &DFS accounting structure to free
+ */
+void DepthFirstSearchAnalysis_Finalize(ofstream *out, struct DFSAccounting &DFS)
+{
+  delete (DFS.AtomStack);
+  // delete (DFS.BackEdgeStack); // DON'T free, see DepthFirstSearchAnalysis(), is returned as allocated
+}
+;
 
 /** Performs a Depth-First search on this molecule.
@@ -332 +481 @@
  * \return list of each disconnected subgraph as an individual molecule class structure
  */
-MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(ofstream *out, class StackClass<bond *> *&BackEdgeStack)
-{
-  class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
+MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(ofstream *out, class StackClass<bond *> *&BackEdgeStack) const
+{
+  struct DFSAccounting DFS;
   BackEdgeStack = new StackClass<bond *> (BondCount);
+  DFS.BackEdgeStack = BackEdgeStack;
   MoleculeLeafClass *SubGraphs = new MoleculeLeafClass(NULL);
   MoleculeLeafClass *LeafWalker = SubGraphs;
-  int CurrentGraphNr = 0, OldGraphNr;
-  int ComponentNumber = 0;
-  atom *Walker = NULL, *OtherAtom = NULL, *Root = start->next;
+  int OldGraphNr = 0;
+  atom *Walker = NULL;
   bond *Binder = NULL;
-  bool BackStepping = false;
 
   *out << Verbose(0) << "Begin of DepthFirstSearchAnalysis" << endl;
-
-  ResetAllBondsToUnused();
-  ResetAllAtomNumbers();
-  InitComponentNumbers();
-  BackEdgeStack->ClearStack();
-  while (Root != end) { // if there any atoms at all
-    // (1) mark all edges unused, empty stack, set atom->GraphNr = 0 for all
-    AtomStack->ClearStack();
+  DepthFirstSearchAnalysis_Init(out, DFS, this);
+
+  DFS.Root = start->next;
+  while (DFS.Root != end) { // if there any atoms at all
+    // (1) mark all edges unused, empty stack, set atom->GraphNr = -1 for all
+    DFS.AtomStack->ClearStack();
 
     // put into new subgraph molecule and add this to list of subgraphs
     LeafWalker = new MoleculeLeafClass(LeafWalker);
     LeafWalker->Leaf = new molecule(elemente);
-    LeafWalker->Leaf->AddCopyAtom(Root);
-
-    OldGraphNr = CurrentGraphNr;
-    Walker = Root;
+    LeafWalker->Leaf->AddCopyAtom(DFS.Root);
+
+    OldGraphNr = DFS.CurrentGraphNr;
+    Walker = DFS.Root;
     do { // (10)
       do { // (2) set number and Lowpoint of Atom to i, increase i, push current atom
-        if (!BackStepping) { // if we don't just return from (8)
-          Walker->GraphNr = CurrentGraphNr;
-          Walker->LowpointNr = CurrentGraphNr;
-          *out << Verbose(1) << "Setting Walker[" << Walker->Name << "]'s number to " << Walker->GraphNr << " with Lowpoint " << Walker->LowpointNr << "." << endl;
-          AtomStack->Push(Walker);
-          CurrentGraphNr++;
-        }
-        do { // (3) if Walker has no unused egdes, go to (5)
-          BackStepping = false; // reset backstepping flag for (8)
-          if (Binder == NULL) // if we don't just return from (11), Binder is already set to next unused
-            Binder = FindNextUnused(Walker);
-          if (Binder == NULL)
-            break;
-          *out << Verbose(2) << "Current Unused Bond is " << *Binder << "." << endl;
-          // (4) Mark Binder used, ...
-          Binder->MarkUsed(black);
-          OtherAtom = Binder->GetOtherAtom(Walker);
-          *out << Verbose(2) << "(4) OtherAtom is " << OtherAtom->Name << "." << endl;
-          if (OtherAtom->GraphNr != -1) {
-            // (4a) ... if "other" atom has been visited (GraphNr != 0), set lowpoint to minimum of both, go to (3)
-            Binder->Type = BackEdge;
-            BackEdgeStack->Push(Binder);
-            Walker->LowpointNr = ( Walker->LowpointNr < OtherAtom->GraphNr ) ? Walker->LowpointNr : OtherAtom->GraphNr;
-            *out << Verbose(3) << "(4a) Visited: Setting Lowpoint of Walker[" << Walker->Name << "] to " << Walker->LowpointNr << "." << endl;
-          } else {
-            // (4b) ... otherwise set OtherAtom as Ancestor of Walker and Walker as OtherAtom, go to (2)
-            Binder->Type = TreeEdge;
-            OtherAtom->Ancestor = Walker;
-            Walker = OtherAtom;
-            *out << Verbose(3) << "(4b) Not Visited: OtherAtom[" << OtherAtom->Name << "]'s Ancestor is now " << OtherAtom->Ancestor->Name << ", Walker is OtherAtom " << OtherAtom->Name << "." << endl;
-            break;
-          }
-          Binder = NULL;
-        } while (1);  // (3)
+        DepthFirstSearchAnalysis_SetWalkersGraphNr(out, Walker, DFS);
+
+        DepthFirstSearchAnalysis_ProbeAlongUnusedBond(out, this, Walker, Binder, DFS);
+
         if (Binder == NULL) {
           *out << Verbose(2) << "No more Unused Bonds." << endl;
@@ -402 +518 @@
         } else
           Binder = NULL;
-      } while (1);  // (2)
+      } while (1); // (2)
 
       // if we came from backstepping, yet there were no more unused bonds, we end up here with no Ancestor, because Walker is Root! Then we are finished!
-      if ((Walker == Root) && (Binder == NULL))
+      if ((Walker == DFS.Root) && (Binder == NULL))
         break;
 
-      // (5) if Ancestor of Walker is ...
-      *out << Verbose(1) << "(5) Number of Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "] is " << Walker->Ancestor->GraphNr << "." << endl;
-      if (Walker->Ancestor->GraphNr != Root->GraphNr) {
-        // (6)  (Ancestor of Walker is not Root)
-        if (Walker->LowpointNr < Walker->Ancestor->GraphNr) {
-          // (6a) set Ancestor's Lowpoint number to minimum of of its Ancestor and itself, go to Step(8)
-          Walker->Ancestor->LowpointNr = (Walker->Ancestor->LowpointNr < Walker->LowpointNr) ? Walker->Ancestor->LowpointNr : Walker->LowpointNr;
-          *out << Verbose(2) << "(6) Setting Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "]'s Lowpoint to " << Walker->Ancestor->LowpointNr << "." << endl;
-        } else {
-          // (7) (Ancestor of Walker is a separating vertex, remove all from stack till Walker (including), these and Ancestor form a component
-          Walker->Ancestor->SeparationVertex = true;
-          *out << Verbose(2) << "(7) Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "]'s is a separating vertex, creating component." << endl;
-          SetNextComponentNumber(Walker->Ancestor, ComponentNumber);
-          *out << Verbose(3) << "(7) Walker[" << Walker->Name << "]'s Ancestor's Compont is " << ComponentNumber << "." << endl;
-          SetNextComponentNumber(Walker, ComponentNumber);
-          *out << Verbose(3) << "(7) Walker[" << Walker->Name << "]'s Compont is " << ComponentNumber << "." << endl;
-          do {
-            OtherAtom = AtomStack->PopLast();
-            LeafWalker->Leaf->AddCopyAtom(OtherAtom);
-            SetNextComponentNumber(OtherAtom, ComponentNumber);
-            *out << Verbose(3) << "(7) Other[" << OtherAtom->Name << "]'s Compont is " << ComponentNumber << "." << endl;
-          } while (OtherAtom != Walker);
-          ComponentNumber++;
-        }
-        // (8) Walker becomes its Ancestor, go to (3)
-        *out << Verbose(2) << "(8) Walker[" << Walker->Name << "] is now its Ancestor " << Walker->Ancestor->Name << ", backstepping. " << endl;
-        Walker = Walker->Ancestor;
-        BackStepping = true;
-      }
-      if (!BackStepping) {  // coming from (8) want to go to (3)
-        // (9) remove all from stack till Walker (including), these and Root form a component
-        AtomStack->Output(out);
-        SetNextComponentNumber(Root, ComponentNumber);
-        *out << Verbose(3) << "(9) Root[" << Root->Name << "]'s Component is " << ComponentNumber << "." << endl;
-        SetNextComponentNumber(Walker, ComponentNumber);
-        *out << Verbose(3) << "(9) Walker[" << Walker->Name << "]'s Component is " << ComponentNumber << "." << endl;
-        do {
-          OtherAtom = AtomStack->PopLast();
-          LeafWalker->Leaf->AddCopyAtom(OtherAtom);
-          SetNextComponentNumber(OtherAtom, ComponentNumber);
-          *out << Verbose(3) << "(7) Other[" << OtherAtom->Name << "]'s Compont is " << ComponentNumber << "." << endl;
-        } while (OtherAtom != Walker);
-        ComponentNumber++;
-
-        // (11) Root is separation vertex,  set Walker to Root and go to (4)
-        Walker = Root;
-        Binder = FindNextUnused(Walker);
-        *out << Verbose(1) << "(10) Walker is Root[" << Root->Name << "], next Unused Bond is " << Binder << "." << endl;
-        if (Binder != NULL) { // Root is separation vertex
-          *out << Verbose(1) << "(11) Root is a separation vertex." << endl;
-          Walker->SeparationVertex = true;
-        }
-      }
-    } while ((BackStepping) || (Binder != NULL)); // (10) halt only if Root has no unused edges
+      DepthFirstSearchAnalysis_CheckForaNewComponent(out, this, Walker, DFS, LeafWalker);
+
+      DepthFirstSearchAnalysis_CleanRootStackDownTillWalker(out, this, Walker, Binder, DFS, LeafWalker);
+
+    } while ((DFS.BackStepping) || (Binder != NULL)); // (10) halt only if Root has no unused edges
 
     // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph
-    *out << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << CurrentGraphNr << "." << endl;
+    *out << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << DFS.CurrentGraphNr << "." << endl;
     LeafWalker->Leaf->Output(out);
     *out << endl;
 
     // step on to next root
-    while ((Root != end) && (Root->GraphNr != -1)) {
+    while ((DFS.Root != end) && (DFS.Root->GraphNr != -1)) {
       //*out << Verbose(1) << "Current next subgraph root candidate is " << Root->Name << "." << endl;
-      if (Root->GraphNr != -1) // if already discovered, step on
-        Root = Root->next;
+      if (DFS.Root->GraphNr != -1) // if already discovered, step on
+        DFS.Root = DFS.Root->next;
     }
   }
   // set cyclic bond criterium on "same LP" basis
-  Binder = first;
-  while(Binder->next != last) {
+  CyclicBondAnalysis();
+
+  OutputGraphInfoPerAtom(out);
+
+  OutputGraphInfoPerBond(out);
+
+  // free all and exit
+  DepthFirstSearchAnalysis_Finalize(out, DFS);
+  *out << Verbose(0) << "End of DepthFirstSearchAnalysis" << endl;
+  return SubGraphs;
+}
+;
+
+/** Scans through all bonds and set bond::Cyclic to true where atom::LowpointNr of both ends is equal: LP criterion.
+ */
+void molecule::CyclicBondAnalysis() const
+{
+  NoCyclicBonds = 0;
+  bond *Binder = first;
+  while (Binder->next != last) {
     Binder = Binder->next;
     if (Binder->rightatom->LowpointNr == Binder->leftatom->LowpointNr) { // cyclic ??
@@ -484 +569 @@
     }
   }
-
-
+}
+;
+
+/** Output graph information per atom.
+ * \param *out output stream
+ */
+void molecule::OutputGraphInfoPerAtom(ofstream *out) const
+{
   *out << Verbose(1) << "Final graph info for each atom is:" << endl;
-  Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    *out << Verbose(2) << "Atom " << Walker->Name << " is " << ((Walker->SeparationVertex) ? "a" : "not a") << " separation vertex, components are ";
-    OutputComponentNumber(out, Walker);
-    *out << " with Lowpoint " << Walker->LowpointNr << " and Graph Nr. " << Walker->GraphNr << "." << endl;
-  }
-
+  ActOnAllAtoms(&atom::OutputGraphInfo, out);
+}
+;
+
+/** Output graph information per bond.
+ * \param *out output stream
+ */
+void molecule::OutputGraphInfoPerBond(ofstream *out) const
+{
   *out << Verbose(1) << "Final graph info for each bond is:" << endl;
-  Binder = first;
-  while(Binder->next != last) {
+  bond *Binder = first;
+  while (Binder->next != last) {
     Binder = Binder->next;
     *out << Verbose(2) << ((Binder->Type == TreeEdge) ? "TreeEdge " : "BackEdge ") << *Binder << ": <";
     *out << ((Binder->leftatom->SeparationVertex) ? "SP," : "") << "L" << Binder->leftatom->LowpointNr << " G" << Binder->leftatom->GraphNr << " Comp.";
-    OutputComponentNumber(out, Binder->leftatom);
+    Binder->leftatom->OutputComponentNumber(out);
     *out << " ===  ";
     *out << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp.";
-    OutputComponentNumber(out, Binder->rightatom);
+    Binder->rightatom->OutputComponentNumber(out);
     *out << ">." << endl;
     if (Binder->Cyclic) // cyclic ??
       *out << Verbose(3) << "Lowpoint at each side are equal: CYCLIC!" << endl;
   }
-
-  // free all and exit
-  delete(AtomStack);
-  *out << Verbose(0) << "End of DepthFirstSearchAnalysis" << endl;
-  return SubGraphs;
+}
+;
+
+/** Initialise each vertex as white with no predecessor, empty queue, color Root lightgray.
+ * \param *out output stream for debugging
+ * \param &BFS accounting structure
+ * \param AtomCount number of entries in the array to allocate
+ */
+void InitializeBFSAccounting(ofstream *out, struct BFSAccounting &BFS, int AtomCount)
+{
+  BFS.AtomCount = AtomCount;
+  BFS.PredecessorList = Calloc<atom*> (AtomCount, "molecule::BreadthFirstSearchAdd_Init: **PredecessorList");
+  BFS.ShortestPathList = Malloc<int> (AtomCount, "molecule::BreadthFirstSearchAdd_Init: *ShortestPathList");
+  BFS.ColorList = Calloc<enum Shading> (AtomCount, "molecule::BreadthFirstSearchAdd_Init: *ColorList");
+  BFS.BFSStack = new StackClass<atom *> (AtomCount);
+
+  for (int i = AtomCount; i--;)
+    BFS.ShortestPathList[i] = -1;
 };
+
+/** Free's accounting structure.
+ * \param *out output stream for debugging
+ * \param &BFS accounting structure
+ */
+void FinalizeBFSAccounting(ofstream *out, struct BFSAccounting &BFS)
+{
+  Free(&BFS.PredecessorList);
+  Free(&BFS.ShortestPathList);
+  Free(&BFS.ColorList);
+  delete (BFS.BFSStack);
+  BFS.AtomCount = 0;
+};
+
+/** Clean the accounting structure.
+ * \param *out output stream for debugging
+ * \param &BFS accounting structure
+ */
+void CleanBFSAccounting(ofstream *out, struct BFSAccounting &BFS)
+{
+  atom *Walker = NULL;
+  while (!BFS.TouchedStack->IsEmpty()) {
+    Walker = BFS.TouchedStack->PopFirst();
+    BFS.PredecessorList[Walker->nr] = NULL;
+    BFS.ShortestPathList[Walker->nr] = -1;
+    BFS.ColorList[Walker->nr] = white;
+  }
+};
+
+/** Resets shortest path list and BFSStack.
+ * \param *out output stream for debugging
+ * \param *&Walker current node, pushed onto BFSAccounting::BFSStack and BFSAccounting::TouchedStack
+ * \param &BFS accounting structure
+ */
+void ResetBFSAccounting(ofstream *out, atom *&Walker, struct BFSAccounting &BFS)
+{
+  BFS.ShortestPathList[Walker->nr] = 0;
+  BFS.BFSStack->ClearStack(); // start with empty BFS stack
+  BFS.BFSStack->Push(Walker);
+  BFS.TouchedStack->Push(Walker);
+};
+
+/** Performs a BFS from \a *Root, trying to find the same node and hence a cycle.
+ * \param *out output stream for debugging
+ * \param *&BackEdge the edge from root that we don't want to move along
+ * \param &BFS accounting structure
+ */
+void CyclicStructureAnalysis_CyclicBFSFromRootToRoot(ofstream *out, bond *&BackEdge, struct BFSAccounting &BFS)
+{
+  atom *Walker = NULL;
+  atom *OtherAtom = NULL;
+  do { // look for Root
+    Walker = BFS.BFSStack->PopFirst();
+    *out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *BFS.Root << "." << endl;
+    for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+      if ((*Runner) != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
+        OtherAtom = (*Runner)->GetOtherAtom(Walker);
+#ifdef ADDHYDROGEN
+        if (OtherAtom->type->Z != 1) {
+#endif
+        *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *(*Runner) << "." << endl;
+        if (BFS.ColorList[OtherAtom->nr] == white) {
+          BFS.TouchedStack->Push(OtherAtom);
+          BFS.ColorList[OtherAtom->nr] = lightgray;
+          BFS.PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
+          BFS.ShortestPathList[OtherAtom->nr] = BFS.ShortestPathList[Walker->nr] + 1;
+          *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << BFS.ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
+          //if (BFS.ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
+          *out << Verbose(3) << "Putting OtherAtom into queue." << endl;
+          BFS.BFSStack->Push(OtherAtom);
+          //}
+        } else {
+          *out << Verbose(3) << "Not Adding, has already been visited." << endl;
+        }
+        if (OtherAtom == BFS.Root)
+          break;
+#ifdef ADDHYDROGEN
+      } else {
+        *out << Verbose(2) << "Skipping hydrogen atom " << *OtherAtom << "." << endl;
+        BFS.ColorList[OtherAtom->nr] = black;
+      }
+#endif
+      } else {
+        *out << Verbose(2) << "Bond " << *(*Runner) << " not Visiting, is the back edge." << endl;
+      }
+    }
+    BFS.ColorList[Walker->nr] = black;
+    *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
+    if (OtherAtom == BFS.Root) { // if we have found the root, check whether this cycle wasn't already found beforehand
+      // step through predecessor list
+      while (OtherAtom != BackEdge->rightatom) {
+        if (!OtherAtom->GetTrueFather()->IsCyclic) // if one bond in the loop is not marked as cyclic, we haven't found this cycle yet
+          break;
+        else
+          OtherAtom = BFS.PredecessorList[OtherAtom->nr];
+      }
+      if (OtherAtom == BackEdge->rightatom) { // if each atom in found cycle is cyclic, loop's been found before already
+        *out << Verbose(3) << "This cycle was already found before, skipping and removing seeker from search." << endl;
+        do {
+          OtherAtom = BFS.TouchedStack->PopLast();
+          if (BFS.PredecessorList[OtherAtom->nr] == Walker) {
+            *out << Verbose(4) << "Removing " << *OtherAtom << " from lists and stacks." << endl;
+            BFS.PredecessorList[OtherAtom->nr] = NULL;
+            BFS.ShortestPathList[OtherAtom->nr] = -1;
+            BFS.ColorList[OtherAtom->nr] = white;
+            BFS.BFSStack->RemoveItem(OtherAtom);
+          }
+        } while ((!BFS.TouchedStack->IsEmpty()) && (BFS.PredecessorList[OtherAtom->nr] == NULL));
+        BFS.TouchedStack->Push(OtherAtom); // last was wrongly popped
+        OtherAtom = BackEdge->rightatom; // set to not Root
+      } else
+        OtherAtom = BFS.Root;
+    }
+  } while ((!BFS.BFSStack->IsEmpty()) && (OtherAtom != BFS.Root) && (OtherAtom != NULL)); // || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr])));
+};
+
+/** Climb back the BFSAccounting::PredecessorList and find cycle members.
+ * \param *out output stream for debugging
+ * \param *&OtherAtom
+ * \param *&BackEdge denotes the edge we did not want to travel along when doing CyclicBFSFromRootToRoot()
+ * \param &BFS accounting structure
+ * \param *&MinimumRingSize minimum distance from this node possible without encountering oneself, set on return for each atom
+ * \param &MinRingSize global minimum distance from one node without encountering oneself, set on return
+ */
+void CyclicStructureAnalysis_RetrieveCycleMembers(ofstream *out, atom *&OtherAtom, bond *&BackEdge, struct BFSAccounting &BFS, int *&MinimumRingSize, int &MinRingSize)
+{
+  atom *Walker = NULL;
+  int NumCycles = 0;
+  int RingSize = -1;
+
+  if (OtherAtom == BFS.Root) {
+    // now climb back the predecessor list and thus find the cycle members
+    NumCycles++;
+    RingSize = 1;
+    BFS.Root->GetTrueFather()->IsCyclic = true;
+    *out << Verbose(1) << "Found ring contains: ";
+    Walker = BFS.Root;
+    while (Walker != BackEdge->rightatom) {
+      *out << Walker->Name << " <-> ";
+      Walker = BFS.PredecessorList[Walker->nr];
+      Walker->GetTrueFather()->IsCyclic = true;
+      RingSize++;
+    }
+    *out << Walker->Name << "  with a length of " << RingSize << "." << endl << endl;
+    // walk through all and set MinimumRingSize
+    Walker = BFS.Root;
+    MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
+    while (Walker != BackEdge->rightatom) {
+      Walker = BFS.PredecessorList[Walker->nr];
+      if (RingSize < MinimumRingSize[Walker->GetTrueFather()->nr])
+        MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
+    }
+    if ((RingSize < MinRingSize) || (MinRingSize == -1))
+      MinRingSize = RingSize;
+  } else {
+    *out << Verbose(1) << "No ring containing " << *BFS.Root << " with length equal to or smaller than " << MinimumRingSize[Walker->GetTrueFather()->nr] << " found." << endl;
+  }
+};
+
+/** From a given node performs a BFS to touch the next cycle, for whose nodes \a *&MinimumRingSize is set and set it accordingly.
+ * \param *out output stream for debugging
+ * \param *&Root node to look for closest cycle from, i.e. \a *&MinimumRingSize is set for this node
+ * \param *&MinimumRingSize minimum distance from this node possible without encountering oneself, set on return for each atom
+ * \param AtomCount number of nodes in graph
+ */
+void CyclicStructureAnalysis_BFSToNextCycle(ofstream *out, atom *&Root, atom *&Walker, int *&MinimumRingSize, int AtomCount)
+{
+  struct BFSAccounting BFS;
+  atom *OtherAtom = Walker;
+
+  InitializeBFSAccounting(out, BFS, AtomCount);
+
+  ResetBFSAccounting(out, Walker, BFS);
+  while (OtherAtom != NULL) { // look for Root
+    Walker = BFS.BFSStack->PopFirst();
+    //*out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
+    for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+      // "removed (*Runner) != BackEdge) || " from next if, is u
+      if ((Walker->ListOfBonds.size() == 1)) { // only walk along DFS spanning tree (otherwise we always find SP of 1 being backedge Binder), but terminal hydrogens may be connected via backedge, hence extra check
+        OtherAtom = (*Runner)->GetOtherAtom(Walker);
+        //*out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
+        if (BFS.ColorList[OtherAtom->nr] == white) {
+          BFS.TouchedStack->Push(OtherAtom);
+          BFS.ColorList[OtherAtom->nr] = lightgray;
+          BFS.PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
+          BFS.ShortestPathList[OtherAtom->nr] = BFS.ShortestPathList[Walker->nr] + 1;
+          //*out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
+          if (OtherAtom->GetTrueFather()->IsCyclic) { // if the other atom is connected to a ring
+            MinimumRingSize[Root->GetTrueFather()->nr] = BFS.ShortestPathList[OtherAtom->nr] + MinimumRingSize[OtherAtom->GetTrueFather()->nr];
+            OtherAtom = NULL; //break;
+            break;
+          } else
+            BFS.BFSStack->Push(OtherAtom);
+        } else {
+          //*out << Verbose(3) << "Not Adding, has already been visited." << endl;
+        }
+      } else {
+        //*out << Verbose(3) << "Not Visiting, is a back edge." << endl;
+      }
+    }
+    BFS.ColorList[Walker->nr] = black;
+    //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
+  }
+  //CleanAccountingLists(TouchedStack, PredecessorList, ShortestPathList, ColorList);
+
+  FinalizeBFSAccounting(out, BFS);
+}
+;
+
+/** All nodes that are not in cycles get assigned a \a *&MinimumRingSizeby BFS to next cycle.
+ * \param *out output stream for debugging
+ * \param *&MinimumRingSize array with minimum distance without encountering onself for each atom
+ * \param &MinRingSize global minium distance
+ * \param &NumCyles number of cycles in graph
+ * \param *mol molecule with atoms
+ */
+void CyclicStructureAnalysis_AssignRingSizetoNonCycleMembers(ofstream *out, int *&MinimumRingSize, int &MinRingSize, int &NumCycles, const molecule * const mol)
+{
+  atom *Root = NULL;
+  atom *Walker = NULL;
+  if (MinRingSize != -1) { // if rings are present
+    // go over all atoms
+    Root = mol->start;
+    while (Root->next != mol->end) {
+      Root = Root->next;
+
+      if (MinimumRingSize[Root->GetTrueFather()->nr] == mol->AtomCount) { // check whether MinimumRingSize is set, if not BFS to next where it is
+        Walker = Root;
+
+        //*out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
+        CyclicStructureAnalysis_BFSToNextCycle(out, Root, Walker, MinimumRingSize, mol->AtomCount);
+
+      }
+      *out << Verbose(1) << "Minimum ring size of " << *Root << " is " << MinimumRingSize[Root->GetTrueFather()->nr] << "." << endl;
+    }
+    *out << Verbose(1) << "Minimum ring size is " << MinRingSize << ", over " << NumCycles << " cycles total." << endl;
+  } else
+    *out << Verbose(1) << "No rings were detected in the molecular structure." << endl;
+}
+;
 
 /** Analyses the cycles found and returns minimum of all cycle lengths.
@@ -526 +871 @@
  * \todo BFS from the not-same-LP to find back to starting point of tributary cycle over more than one bond
  */
-void molecule::CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *  BackEdgeStack, int *&MinimumRingSize)
-{
-  atom **PredecessorList = Malloc<atom*>(AtomCount, "molecule::CyclicStructureAnalysis: **PredecessorList");
-  int *ShortestPathList = Malloc<int>(AtomCount, "molecule::CyclicStructureAnalysis: *ShortestPathList");
-  enum Shading *ColorList = Malloc<enum Shading>(AtomCount, "molecule::CyclicStructureAnalysis: *ColorList");
-  class StackClass<atom *> *BFSStack = new StackClass<atom *> (AtomCount);   // will hold the current ring
-  class StackClass<atom *> *TouchedStack = new StackClass<atom *> (AtomCount);   // contains all "touched" atoms (that need to be reset after BFS loop)
-  atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL;
-  bond *Binder = NULL, *BackEdge = NULL;
-  int RingSize, NumCycles, MinRingSize = -1;
-
-  // initialise each vertex as white with no predecessor, empty queue, color Root lightgray
-  for (int i=AtomCount;i--;) {
-    PredecessorList[i] = NULL;
-    ShortestPathList[i] = -1;
-    ColorList[i] = white;
-  }
-
-  *out << Verbose(1) << "Back edge list - ";
-  BackEdgeStack->Output(out);
+void molecule::CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> * BackEdgeStack, int *&MinimumRingSize) const
+{
+  struct BFSAccounting BFS;
+  atom *Walker = NULL;
+  atom *OtherAtom = NULL;
+  bond *BackEdge = NULL;
+  int NumCycles = 0;
+  int MinRingSize = -1;
+
+  InitializeBFSAccounting(out, BFS, AtomCount);
+
+  //*out << Verbose(1) << "Back edge list - ";
+  //BackEdgeStack->Output(out);
 
   *out << Verbose(1) << "Analysing cycles ... " << endl;
@@ -552 +890 @@
     BackEdge = BackEdgeStack->PopFirst();
     // this is the target
-    Root = BackEdge->leftatom;
+    BFS.Root = BackEdge->leftatom;
     // this is the source point
     Walker = BackEdge->rightatom;
-    ShortestPathList[Walker->nr] = 0;
-    BFSStack->ClearStack();  // start with empty BFS stack
-    BFSStack->Push(Walker);
-    TouchedStack->Push(Walker);
+
+    ResetBFSAccounting(out, Walker, BFS);
+
     *out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
     OtherAtom = NULL;
-    do {  // look for Root
-      Walker = BFSStack->PopFirst();
-      *out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
-      for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
-        Binder = ListOfBondsPerAtom[Walker->nr][i];
-        if (Binder != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
-          OtherAtom = Binder->GetOtherAtom(Walker);
-#ifdef ADDHYDROGEN
-          if (OtherAtom->type->Z != 1) {
-#endif
-            *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
-            if (ColorList[OtherAtom->nr] == white) {
-              TouchedStack->Push(OtherAtom);
-              ColorList[OtherAtom->nr] = lightgray;
-              PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
-              ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
-              *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
-              //if (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
-                *out << Verbose(3) << "Putting OtherAtom into queue." << endl;
-                BFSStack->Push(OtherAtom);
-              //}
-            } else {
-              *out << Verbose(3) << "Not Adding, has already been visited." << endl;
-            }
-            if (OtherAtom == Root)
-              break;
-#ifdef ADDHYDROGEN
-          } else {
-            *out << Verbose(2) << "Skipping hydrogen atom " << *OtherAtom << "." << endl;
-            ColorList[OtherAtom->nr] = black;
-          }
-#endif
-        } else {
-          *out << Verbose(2) << "Bond " << *Binder << " not Visiting, is the back edge." << endl;
-        }
-      }
-      ColorList[Walker->nr] = black;
-      *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
-      if (OtherAtom == Root) { // if we have found the root, check whether this cycle wasn't already found beforehand
-        // step through predecessor list
-        while (OtherAtom != BackEdge->rightatom) {
-          if (!OtherAtom->GetTrueFather()->IsCyclic)  // if one bond in the loop is not marked as cyclic, we haven't found this cycle yet
-            break;
-          else
-            OtherAtom = PredecessorList[OtherAtom->nr];
-        }
-        if (OtherAtom == BackEdge->rightatom) { // if each atom in found cycle is cyclic, loop's been found before already
-          *out << Verbose(3) << "This cycle was already found before, skipping and removing seeker from search." << endl;\
-          do {
-            OtherAtom = TouchedStack->PopLast();
-            if (PredecessorList[OtherAtom->nr] == Walker) {
-              *out << Verbose(4) << "Removing " << *OtherAtom << " from lists and stacks." << endl;
-              PredecessorList[OtherAtom->nr] = NULL;
-              ShortestPathList[OtherAtom->nr] = -1;
-              ColorList[OtherAtom->nr] = white;
-              BFSStack->RemoveItem(OtherAtom);
-            }
-          } while ((!TouchedStack->IsEmpty()) && (PredecessorList[OtherAtom->nr] == NULL));
-          TouchedStack->Push(OtherAtom);  // last was wrongly popped
-          OtherAtom = BackEdge->rightatom; // set to not Root
-        } else
-          OtherAtom = Root;
-      }
-    } while ((!BFSStack->IsEmpty()) && (OtherAtom != Root) && (OtherAtom != NULL)); // || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr])));
-
-    if (OtherAtom == Root) {
-      // now climb back the predecessor list and thus find the cycle members
-      NumCycles++;
-      RingSize = 1;
-      Root->GetTrueFather()->IsCyclic = true;
-      *out << Verbose(1) << "Found ring contains: ";
-      Walker = Root;
-      while (Walker != BackEdge->rightatom) {
-        *out << Walker->Name << " <-> ";
-        Walker = PredecessorList[Walker->nr];
-        Walker->GetTrueFather()->IsCyclic = true;
-        RingSize++;
-      }
-      *out << Walker->Name << "  with a length of " << RingSize << "." << endl << endl;
-      // walk through all and set MinimumRingSize
-      Walker = Root;
-      MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
-      while (Walker != BackEdge->rightatom) {
-        Walker = PredecessorList[Walker->nr];
-        if (RingSize < MinimumRingSize[Walker->GetTrueFather()->nr])
-          MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
-      }
-      if ((RingSize < MinRingSize) || (MinRingSize == -1))
-        MinRingSize = RingSize;
-    } else {
-      *out << Verbose(1) << "No ring containing " << *Root << " with length equal to or smaller than " << MinimumRingSize[Walker->GetTrueFather()->nr] << " found." << endl;
-    }
-
-    // now clean the lists
-    while (!TouchedStack->IsEmpty()){
-      Walker = TouchedStack->PopFirst();
-      PredecessorList[Walker->nr] = NULL;
-      ShortestPathList[Walker->nr] = -1;
-      ColorList[Walker->nr] = white;
-    }
-  }
-  if (MinRingSize != -1) {
-    // go over all atoms
-    Root = start;
-    while(Root->next != end) {
-      Root = Root->next;
-
-      if (MinimumRingSize[Root->GetTrueFather()->nr] == AtomCount) { // check whether MinimumRingSize is set, if not BFS to next where it is
-        Walker = Root;
-        ShortestPathList[Walker->nr] = 0;
-        BFSStack->ClearStack();  // start with empty BFS stack
-        BFSStack->Push(Walker);
-        TouchedStack->Push(Walker);
-        //*out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
-        OtherAtom = Walker;
-        while (OtherAtom != NULL) {  // look for Root
-          Walker = BFSStack->PopFirst();
-          //*out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
-          for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
-            Binder = ListOfBondsPerAtom[Walker->nr][i];
-            if ((Binder != BackEdge) || (NumberOfBondsPerAtom[Walker->nr] == 1)) { // only walk along DFS spanning tree (otherwise we always find SP of 1 being backedge Binder), but terminal hydrogens may be connected via backedge, hence extra check
-              OtherAtom = Binder->GetOtherAtom(Walker);
-              //*out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
-              if (ColorList[OtherAtom->nr] == white) {
-                TouchedStack->Push(OtherAtom);
-                ColorList[OtherAtom->nr] = lightgray;
-                PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
-                ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
-                //*out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
-                if (OtherAtom->GetTrueFather()->IsCyclic) { // if the other atom is connected to a ring
-                  MinimumRingSize[Root->GetTrueFather()->nr] = ShortestPathList[OtherAtom->nr]+MinimumRingSize[OtherAtom->GetTrueFather()->nr];
-                  OtherAtom = NULL; //break;
-                  break;
-                } else
-                  BFSStack->Push(OtherAtom);
-              } else {
-                //*out << Verbose(3) << "Not Adding, has already been visited." << endl;
-              }
-            } else {
-              //*out << Verbose(3) << "Not Visiting, is a back edge." << endl;
-            }
-          }
-          ColorList[Walker->nr] = black;
-          //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
-        }
-
-        // now clean the lists
-        while (!TouchedStack->IsEmpty()){
-          Walker = TouchedStack->PopFirst();
-          PredecessorList[Walker->nr] = NULL;
-          ShortestPathList[Walker->nr] = -1;
-          ColorList[Walker->nr] = white;
-        }
-      }
-      *out << Verbose(1) << "Minimum ring size of " << *Root << " is " << MinimumRingSize[Root->GetTrueFather()->nr] << "." << endl;
-    }
-    *out << Verbose(1) << "Minimum ring size is " << MinRingSize << ", over " << NumCycles << " cycles total." << endl;
-  } else
-    *out << Verbose(1) << "No rings were detected in the molecular structure." << endl;
-
-  Free(&PredecessorList);
-  Free(&ShortestPathList);
-  Free(&ColorList);
-  delete(BFSStack);
+    CyclicStructureAnalysis_CyclicBFSFromRootToRoot(out, BackEdge, BFS);
+
+    CyclicStructureAnalysis_RetrieveCycleMembers(out, OtherAtom, BackEdge, BFS, MinimumRingSize, MinRingSize);
+
+    CleanBFSAccounting(out, BFS);
+  }
+  FinalizeBFSAccounting(out, BFS);
+
+  CyclicStructureAnalysis_AssignRingSizetoNonCycleMembers(out, MinimumRingSize, MinRingSize, NumCycles, this);
 };
 
@@ -732 +914 @@
  * \param nr number to use
  */
-void molecule::SetNextComponentNumber(atom *vertex, int nr)
-{
-  int i=0;
+void molecule::SetNextComponentNumber(atom *vertex, int nr) const
+{
+  size_t i = 0;
   if (vertex != NULL) {
-    for(;i<NumberOfBondsPerAtom[vertex->nr];i++) {
-      if (vertex->ComponentNr[i] == -1) {   // check if not yet used
+    for (; i < vertex->ListOfBonds.size(); i++) {
+      if (vertex->ComponentNr[i] == -1) { // check if not yet used
         vertex->ComponentNr[i] = nr;
         break;
-      }
-      else if (vertex->ComponentNr[i] == nr) // if number is already present, don't add another time
-        break;  // breaking here will not cause error!
-    }
-    if (i == NumberOfBondsPerAtom[vertex->nr])
+      } else if (vertex->ComponentNr[i] == nr) // if number is already present, don't add another time
+        break; // breaking here will not cause error!
+    }
+    if (i == vertex->ListOfBonds.size())
       cerr << "Error: All Component entries are already occupied!" << endl;
   } else
-      cerr << "Error: Given vertex is NULL!" << endl;
-};
-
-/** Output a list of flags, stating whether the bond was visited or not.
- * \param *out output stream for debugging
- */
-void molecule::OutputComponentNumber(ofstream *out, atom *vertex)
-{
-  for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
-    *out << vertex->ComponentNr[i] << "  ";
-};
-
-/** Allocates memory for all atom::*ComponentNr in this molecule and sets each entry to -1.
- */
-void molecule::InitComponentNumbers()
-{
-  atom *Walker = start;
-  while(Walker->next != end) {
-    Walker = Walker->next;
-    if (Walker->ComponentNr != NULL)
-      Free(&Walker->ComponentNr);
-    Walker->ComponentNr = Malloc<int>(NumberOfBondsPerAtom[Walker->nr], "molecule::InitComponentNumbers: *Walker->ComponentNr");
-    for (int i=NumberOfBondsPerAtom[Walker->nr];i--;)
-      Walker->ComponentNr[i] = -1;
-  }
-};
+    cerr << "Error: Given vertex is NULL!" << endl;
+}
+;
 
 /** Returns next unused bond for this atom \a *vertex or NULL of none exists.
@@ -778 +936 @@
  * \return bond class or NULL
  */
-bond * molecule::FindNextUnused(atom *vertex)
-{
-  for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
-    if (ListOfBondsPerAtom[vertex->nr][i]->IsUsed() == white)
-      return(ListOfBondsPerAtom[vertex->nr][i]);
+bond * molecule::FindNextUnused(atom *vertex) const
+{
+  for (BondList::const_iterator Runner = vertex->ListOfBonds.begin(); Runner != vertex->ListOfBonds.end(); (++Runner))
+    if ((*Runner)->IsUsed() == white)
+      return ((*Runner));
   return NULL;
-};
+}
+;
 
 /** Resets bond::Used flag of all bonds in this molecule.
  * \return true - success, false - -failure
  */
-void molecule::ResetAllBondsToUnused()
+void molecule::ResetAllBondsToUnused() const
 {
   bond *Binder = first;
@@ -796 +955 @@
     Binder->ResetUsed();
   }
-};
-
-/** Resets atom::nr to -1 of all atoms in this molecule.
- */
-void molecule::ResetAllAtomNumbers()
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    Walker->GraphNr  = -1;
-  }
-};
+}
+;
 
 /** Output a list of flags, stating whether the bond was visited or not.
@@ -816 +965 @@
 {
   *out << Verbose(4) << "Already Visited Bonds:\t";
-  for(int i=1;i<=list[0];i++) *out << Verbose(0) << list[i] << "  ";
+  for (int i = 1; i <= list[0]; i++)
+    *out << Verbose(0) << list[i] << "  ";
   *out << endl;
-};
-
+}
+;
 
 /** Storing the bond structure of a molecule to file.
@@ -830 +980 @@
 {
   ofstream AdjacencyFile;
-  atom *Walker = NULL;
   stringstream line;
   bool status = true;
@@ -838 +987 @@
   *out << Verbose(1) << "Saving adjacency list ... ";
   if (AdjacencyFile != NULL) {
-    Walker = start;
-    while(Walker->next != end) {
-      Walker = Walker->next;
-      AdjacencyFile << Walker->nr << "\t";
-      for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
-        AdjacencyFile << ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker)->nr << "\t";
-      AdjacencyFile << endl;
-    }
+    ActOnAllAtoms(&atom::OutputAdjacency, &AdjacencyFile);
     AdjacencyFile.close();
     *out << Verbose(1) << "done." << endl;
@@ -854 +996 @@
 
   return status;
-};
+}
+;
+
+bool CheckAdjacencyFileAgainstMolecule_Init(ofstream *out, char *path, ifstream &File, int *&CurrentBonds)
+{
+  stringstream filename;
+  filename << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
+  File.open(filename.str().c_str(), ios::out);
+  *out << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... ";
+  if (File == NULL)
+    return false;
+
+  // allocate storage structure
+  CurrentBonds = Calloc<int> (8, "molecule::CheckAdjacencyFileAgainstMolecule - CurrentBonds"); // contains parsed bonds of current atom
+  return true;
+}
+;
+
+void CheckAdjacencyFileAgainstMolecule_Finalize(ofstream *out, ifstream &File, int *&CurrentBonds)
+{
+  File.close();
+  File.clear();
+  Free(&CurrentBonds);
+}
+;
+
+void CheckAdjacencyFileAgainstMolecule_CompareBonds(ofstream *out, bool &status, int &NonMatchNumber, atom *&Walker, size_t &CurrentBondsOfAtom, int AtomNr, int *&CurrentBonds, atom **ListOfAtoms)
+{
+  size_t j = 0;
+  int id = -1;
+
+  //*out << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
+  if (CurrentBondsOfAtom == Walker->ListOfBonds.size()) {
+    for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+      id = (*Runner)->GetOtherAtom(Walker)->nr;
+      j = 0;
+      for (; (j < CurrentBondsOfAtom) && (CurrentBonds[j++] != id);)
+        ; // check against all parsed bonds
+      if (CurrentBonds[j - 1] != id) { // no match ? Then mark in ListOfAtoms
+        ListOfAtoms[AtomNr] = NULL;
+        NonMatchNumber++;
+        status = false;
+        //*out << "[" << id << "]\t";
+      } else {
+        //*out << id << "\t";
+      }
+    }
+    //*out << endl;
+  } else {
+    *out << "Number of bonds for Atom " << *Walker << " does not match, parsed " << CurrentBondsOfAtom << " against " << Walker->ListOfBonds.size() << "." << endl;
+    status = false;
+  }
+}
+;
 
 /** Checks contents of adjacency file against bond structure in structure molecule.
@@ -865 +1060 @@
 {
   ifstream File;
-  stringstream filename;
   bool status = true;
-  char *buffer = Malloc<char>(MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
-
-  filename << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
-  File.open(filename.str().c_str(), ios::out);
-  *out << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... ";
-  if (File != NULL) {
-    // allocate storage structure
-    int NonMatchNumber = 0;   // will number of atoms with differing bond structure
-    int *CurrentBonds = Malloc<int>(8, "molecule::CheckAdjacencyFileAgainstMolecule - CurrentBonds"); // contains parsed bonds of current atom
-    int CurrentBondsOfAtom;
-
-    // Parse the file line by line and count the bonds
-    while (!File.eof()) {
-      File.getline(buffer, MAXSTRINGSIZE);
-      stringstream line;
-      line.str(buffer);
-      int AtomNr = -1;
-      line >> AtomNr;
-      CurrentBondsOfAtom = -1; // we count one too far due to line end
-      // parse into structure
-      if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
-        while (!line.eof())
-          line >> CurrentBonds[ ++CurrentBondsOfAtom ];
-        // compare against present bonds
-        //cout << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
-        if (CurrentBondsOfAtom == NumberOfBondsPerAtom[AtomNr]) {
-          for(int i=0;i<NumberOfBondsPerAtom[AtomNr];i++) {
-            int id = ListOfBondsPerAtom[AtomNr][i]->GetOtherAtom(ListOfAtoms[AtomNr])->nr;
-            int j = 0;
-            for (;(j<CurrentBondsOfAtom) && (CurrentBonds[j++] != id);); // check against all parsed bonds
-            if (CurrentBonds[j-1] != id) { // no match ? Then mark in ListOfAtoms
-              ListOfAtoms[AtomNr] = NULL;
-              NonMatchNumber++;
-              status = false;
-              //out << "[" << id << "]\t";
-            } else {
-              //out << id << "\t";
-            }
-          }
-          //out << endl;
-        } else {
-          *out << "Number of bonds for Atom " << *ListOfAtoms[AtomNr] << " does not match, parsed " << CurrentBondsOfAtom << " against " << NumberOfBondsPerAtom[AtomNr] << "." << endl;
-          status = false;
-        }
-      }
-    }
-    File.close();
-    File.clear();
-    if (status) { // if equal we parse the KeySetFile
-      *out << Verbose(1) << "done: Equal." << endl;
-      status = true;
-    } else
-      *out << Verbose(1) << "done: Not equal by " << NonMatchNumber << " atoms." << endl;
-    Free(&CurrentBonds);
-  } else {
+  atom *Walker = NULL;
+  char *buffer = NULL;
+  int *CurrentBonds = NULL;
+  int NonMatchNumber = 0; // will number of atoms with differing bond structure
+  size_t CurrentBondsOfAtom = -1;
+
+  if (!CheckAdjacencyFileAgainstMolecule_Init(out, path, File, CurrentBonds)) {
     *out << Verbose(1) << "Adjacency file not found." << endl;
-    status = false;
-  }
-  *out << endl;
+    return true;
+  }
+
+  buffer = Malloc<char> (MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
+  // Parse the file line by line and count the bonds
+  while (!File.eof()) {
+    File.getline(buffer, MAXSTRINGSIZE);
+    stringstream line;
+    line.str(buffer);
+    int AtomNr = -1;
+    line >> AtomNr;
+    CurrentBondsOfAtom = -1; // we count one too far due to line end
+    // parse into structure
+    if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
+      Walker = ListOfAtoms[AtomNr];
+      while (!line.eof())
+        line >> CurrentBonds[++CurrentBondsOfAtom];
+      // compare against present bonds
+      CheckAdjacencyFileAgainstMolecule_CompareBonds(out, status, NonMatchNumber, Walker, CurrentBondsOfAtom, AtomNr, CurrentBonds, ListOfAtoms);
+    }
+  }
   Free(&buffer);
-
+  CheckAdjacencyFileAgainstMolecule_Finalize(out, File, CurrentBonds);
+
+  if (status) { // if equal we parse the KeySetFile
+    *out << Verbose(1) << "done: Equal." << endl;
+  } else
+    *out << Verbose(1) << "done: Not equal by " << NonMatchNumber << " atoms." << endl;
   return status;
-};
-
+}
+;
 
 /** Picks from a global stack with all back edges the ones in the fragment.
@@ -939 +1108 @@
  * \return true - everything ok, false - ReferenceStack was empty
  */
-bool molecule::PickLocalBackEdges(ofstream *out, atom **ListOfLocalAtoms, class StackClass<bond *> *&ReferenceStack, class StackClass<bond *> *&LocalStack)
+bool molecule::PickLocalBackEdges(ofstream *out, atom **ListOfLocalAtoms, class StackClass<bond *> *&ReferenceStack, class StackClass<bond *> *&LocalStack) const
 {
   bool status = true;
@@ -947 +1116 @@
   }
   bond *Binder = ReferenceStack->PopFirst();
-  bond *FirstBond = Binder;   // mark the first bond, so that we don't loop through the stack indefinitely
+  bond *FirstBond = Binder; // mark the first bond, so that we don't loop through the stack indefinitely
   atom *Walker = NULL, *OtherAtom = NULL;
   ReferenceStack->Push(Binder);
 
-  do {  // go through all bonds and push local ones
-    Walker = ListOfLocalAtoms[Binder->leftatom->nr];  // get one atom in the reference molecule
+  do { // go through all bonds and push local ones
+    Walker = ListOfLocalAtoms[Binder->leftatom->nr]; // get one atom in the reference molecule
     if (Walker != NULL) // if this Walker exists in the subgraph ...
-      for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {    // go through the local list of bonds
-        OtherAtom = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
-        if (OtherAtom == ListOfLocalAtoms[Binder->rightatom->nr]) { // found the bond
-          LocalStack->Push(ListOfBondsPerAtom[Walker->nr][i]);
-          *out << Verbose(3) << "Found local edge " << *(ListOfBondsPerAtom[Walker->nr][i]) << "." << endl;
+      for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+        OtherAtom = (*Runner)->GetOtherAtom(Walker);
+        if (OtherAtom == ListOfLocalAtoms[(*Runner)->rightatom->nr]) { // found the bond
+          LocalStack->Push((*Runner));
+          *out << Verbose(3) << "Found local edge " << *(*Runner) << "." << endl;
           break;
         }
       }
-    Binder = ReferenceStack->PopFirst();  // loop the stack for next item
+    Binder = ReferenceStack->PopFirst(); // loop the stack for next item
     *out << Verbose(3) << "Current candidate edge " << Binder << "." << endl;
     ReferenceStack->Push(Binder);
@@ -968 +1137 @@
 
   return status;
-};
-
+}
+;
+
+void BreadthFirstSearchAdd_Init(struct BFSAccounting &BFS, atom *&Root, int AtomCount, int BondOrder, atom **AddedAtomList = NULL)
+{
+  BFS.AtomCount = AtomCount;
+  BFS.BondOrder = BondOrder;
+  BFS.PredecessorList = Calloc<atom*> (AtomCount, "molecule::BreadthFirstSearchAdd_Init: **PredecessorList");
+  BFS.ShortestPathList = Calloc<int> (AtomCount, "molecule::BreadthFirstSearchAdd_Init: *ShortestPathList");
+  BFS.ColorList = Malloc<enum Shading> (AtomCount, "molecule::BreadthFirstSearchAdd_Init: *ColorList");
+  BFS.BFSStack = new StackClass<atom *> (AtomCount);
+
+  BFS.Root = Root;
+  BFS.BFSStack->ClearStack();
+  BFS.BFSStack->Push(Root);
+
+  // initialise each vertex as white with no predecessor, empty queue, color Root lightgray
+  for (int i = AtomCount; i--;) {
+    BFS.ShortestPathList[i] = -1;
+    if ((AddedAtomList != NULL) && (AddedAtomList[i] != NULL)) // mark already present atoms (i.e. Root and maybe others) as visited
+      BFS.ColorList[i] = lightgray;
+    else
+      BFS.ColorList[i] = white;
+  }
+  //BFS.ShortestPathList[Root->nr] = 0; //is set due to Calloc()
+}
+;
+
+void BreadthFirstSearchAdd_Free(struct BFSAccounting &BFS)
+{
+  Free(&BFS.PredecessorList);
+  Free(&BFS.ShortestPathList);
+  Free(&BFS.ColorList);
+  delete (BFS.BFSStack);
+  BFS.AtomCount = 0;
+}
+;
+
+void BreadthFirstSearchAdd_UnvisitedNode(ofstream *out, molecule *Mol, struct BFSAccounting &BFS, atom *&Walker, atom *&OtherAtom, bond *&Binder, bond *&Bond, atom **&AddedAtomList, bond **&AddedBondList, bool IsAngstroem)
+{
+  if (Binder != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem)
+    BFS.ColorList[OtherAtom->nr] = lightgray;
+  BFS.PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
+  BFS.ShortestPathList[OtherAtom->nr] = BFS.ShortestPathList[Walker->nr] + 1;
+  *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " " << ((BFS.ColorList[OtherAtom->nr] == white) ? "white" : "lightgray") << ", its predecessor is " << Walker->Name << " and its Shortest Path is " << BFS.ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
+  if ((((BFS.ShortestPathList[OtherAtom->nr] < BFS.BondOrder) && (Binder != Bond)))) { // Check for maximum distance
+    *out << Verbose(3);
+    if (AddedAtomList[OtherAtom->nr] == NULL) { // add if it's not been so far
+      AddedAtomList[OtherAtom->nr] = Mol->AddCopyAtom(OtherAtom);
+      *out << "Added OtherAtom " << OtherAtom->Name;
+      AddedBondList[Binder->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder);
+      *out << " and bond " << *(AddedBondList[Binder->nr]) << ", ";
+    } else { // this code should actually never come into play (all white atoms are not yet present in BondMolecule, that's why they are white in the first place)
+      *out << "Not adding OtherAtom " << OtherAtom->Name;
+      if (AddedBondList[Binder->nr] == NULL) {
+        AddedBondList[Binder->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder);
+        *out << ", added Bond " << *(AddedBondList[Binder->nr]);
+      } else
+        *out << ", not added Bond ";
+    }
+    *out << ", putting OtherAtom into queue." << endl;
+    BFS.BFSStack->Push(OtherAtom);
+  } else { // out of bond order, then replace
+    if ((AddedAtomList[OtherAtom->nr] == NULL) && (Binder->Cyclic))
+      BFS.ColorList[OtherAtom->nr] = white; // unmark if it has not been queued/added, to make it available via its other bonds (cyclic)
+    if (Binder == Bond)
+      *out << Verbose(3) << "Not Queueing, is the Root bond";
+    else if (BFS.ShortestPathList[OtherAtom->nr] >= BFS.BondOrder)
+      *out << Verbose(3) << "Not Queueing, is out of Bond Count of " << BFS.BondOrder;
+    if (!Binder->Cyclic)
+      *out << ", is not part of a cyclic bond, saturating bond with Hydrogen." << endl;
+    if (AddedBondList[Binder->nr] == NULL) {
+      if ((AddedAtomList[OtherAtom->nr] != NULL)) { // .. whether we add or saturate
+        AddedBondList[Binder->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder);
+      } else {
+#ifdef ADDHYDROGEN
+        if (!Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, IsAngstroem))
+        exit(1);
+#endif
+      }
+    }
+  }
+}
+;
+
+void BreadthFirstSearchAdd_VisitedNode(ofstream *out, molecule *Mol, struct BFSAccounting &BFS, atom *&Walker, atom *&OtherAtom, bond *&Binder, bond *&Bond, atom **&AddedAtomList, bond **&AddedBondList, bool IsAngstroem)
+{
+  *out << Verbose(3) << "Not Adding, has already been visited." << endl;
+  // This has to be a cyclic bond, check whether it's present ...
+  if (AddedBondList[Binder->nr] == NULL) {
+    if ((Binder != Bond) && (Binder->Cyclic) && (((BFS.ShortestPathList[Walker->nr] + 1) < BFS.BondOrder))) {
+      AddedBondList[Binder->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder);
+    } else { // if it's root bond it has to broken (otherwise we would not create the fragments)
+#ifdef ADDHYDROGEN
+      if(!Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, IsAngstroem))
+      exit(1);
+#endif
+    }
+  }
+}
+;
 
 /** Adds atoms up to \a BondCount distance from \a *Root and notes them down in \a **AddedAtomList.
@@ -985 +1253 @@
 void molecule::BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem)
 {
-  atom **PredecessorList = Malloc<atom*>(AtomCount, "molecule::BreadthFirstSearchAdd: **PredecessorList");
-  int *ShortestPathList = Malloc<int>(AtomCount, "molecule::BreadthFirstSearchAdd: *ShortestPathList");
-  enum Shading *ColorList = Malloc<enum Shading>(AtomCount, "molecule::BreadthFirstSearchAdd: *ColorList");
-  class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
+  struct BFSAccounting BFS;
   atom *Walker = NULL, *OtherAtom = NULL;
   bond *Binder = NULL;
 
   // add Root if not done yet
-  AtomStack->ClearStack();
-  if (AddedAtomList[Root->nr] == NULL)  // add Root if not yet present
+  if (AddedAtomList[Root->nr] == NULL) // add Root if not yet present
     AddedAtomList[Root->nr] = Mol->AddCopyAtom(Root);
-  AtomStack->Push(Root);
-
-  // initialise each vertex as white with no predecessor, empty queue, color Root lightgray
-  for (int i=AtomCount;i--;) {
-    PredecessorList[i] = NULL;
-    ShortestPathList[i] = -1;
-    if (AddedAtomList[i] != NULL) // mark already present atoms (i.e. Root and maybe others) as visited
-      ColorList[i] = lightgray;
-    else
-      ColorList[i] = white;
-  }
-  ShortestPathList[Root->nr] = 0;
+
+  BreadthFirstSearchAdd_Init(BFS, Root, BondOrder, AtomCount, AddedAtomList);
 
   // and go on ... Queue always contains all lightgray vertices
-  while (!AtomStack->IsEmpty()) {
+  while (!BFS.BFSStack->IsEmpty()) {
     // 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 = AtomStack->PopFirst(); // pop oldest added
-    *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << NumberOfBondsPerAtom[Walker->nr] << " bonds." << endl;
-    for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
-      Binder = ListOfBondsPerAtom[Walker->nr][i];
-      if (Binder != NULL) { // don't look at bond equal NULL
-        OtherAtom = Binder->GetOtherAtom(Walker);
-        *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
-        if (ColorList[OtherAtom->nr] == white) {
-          if (Binder != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem)
-            ColorList[OtherAtom->nr] = lightgray;
-          PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
-          ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
-          *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " " << ((ColorList[OtherAtom->nr] == white) ? "white" : "lightgray") << ", its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
-          if ((((ShortestPathList[OtherAtom->nr] < BondOrder) && (Binder != Bond))) ) { // Check for maximum distance
-            *out << Verbose(3);
-            if (AddedAtomList[OtherAtom->nr] == NULL) { // add if it's not been so far
-              AddedAtomList[OtherAtom->nr] = Mol->AddCopyAtom(OtherAtom);
-              *out << "Added OtherAtom " << OtherAtom->Name;
-              AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
-              AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
-              AddedBondList[Binder->nr]->Type = Binder->Type;
-              *out << " and bond " << *(AddedBondList[Binder->nr]) << ", ";
-            } else {  // this code should actually never come into play (all white atoms are not yet present in BondMolecule, that's why they are white in the first place)
-              *out << "Not adding OtherAtom " << OtherAtom->Name;
-              if (AddedBondList[Binder->nr] == NULL) {
-                AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
-                AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
-                AddedBondList[Binder->nr]->Type = Binder->Type;
-                *out << ", added Bond " << *(AddedBondList[Binder->nr]);
-              } else
-                *out << ", not added Bond ";
-            }
-            *out << ", putting OtherAtom into queue." << endl;
-            AtomStack->Push(OtherAtom);
-          } else { // out of bond order, then replace
-            if ((AddedAtomList[OtherAtom->nr] == NULL) && (Binder->Cyclic))
-              ColorList[OtherAtom->nr] = white; // unmark if it has not been queued/added, to make it available via its other bonds (cyclic)
-            if (Binder == Bond)
-              *out << Verbose(3) << "Not Queueing, is the Root bond";
-            else if (ShortestPathList[OtherAtom->nr] >= BondOrder)
-              *out << Verbose(3) << "Not Queueing, is out of Bond Count of " << BondOrder;
-            if (!Binder->Cyclic)
-              *out << ", is not part of a cyclic bond, saturating bond with Hydrogen." << endl;
-            if (AddedBondList[Binder->nr] == NULL) {
-              if ((AddedAtomList[OtherAtom->nr] != NULL)) { // .. whether we add or saturate
-                AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
-                AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
-                AddedBondList[Binder->nr]->Type = Binder->Type;
-              } else {
-#ifdef ADDHYDROGEN
-                if (!Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem))
-                  exit(1);
-#endif
-              }
-            }
-          }
+    Walker = BFS.BFSStack->PopFirst(); // pop oldest added
+    *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << Walker->ListOfBonds.size() << " bonds." << endl;
+    for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+      if ((*Runner) != NULL) { // don't look at bond equal NULL
+        Binder = (*Runner);
+        OtherAtom = (*Runner)->GetOtherAtom(Walker);
+        *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *(*Runner) << "." << endl;
+        if (BFS.ColorList[OtherAtom->nr] == white) {
+          BreadthFirstSearchAdd_UnvisitedNode(out, Mol, BFS, Walker, OtherAtom, Binder, Bond, AddedAtomList, AddedBondList, IsAngstroem);
         } else {
-          *out << Verbose(3) << "Not Adding, has already been visited." << endl;
-          // This has to be a cyclic bond, check whether it's present ...
-          if (AddedBondList[Binder->nr] == NULL) {
-            if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) {
-              AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
-              AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
-              AddedBondList[Binder->nr]->Type = Binder->Type;
-            } else { // if it's root bond it has to broken (otherwise we would not create the fragments)
-#ifdef ADDHYDROGEN
-              if(!Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem))
-                exit(1);
-#endif
-            }
-          }
+          BreadthFirstSearchAdd_VisitedNode(out, Mol, BFS, Walker, OtherAtom, Binder, Bond, AddedAtomList, AddedBondList, IsAngstroem);
         }
       }
     }
-    ColorList[Walker->nr] = black;
+    BFS.ColorList[Walker->nr] = black;
     *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
   }
-  Free(&PredecessorList);
-  Free(&ShortestPathList);
-  Free(&ColorList);
-  delete(AtomStack);
-};
-
-/** Adds bond structure to this molecule from \a Father molecule.
- * This basically causes this molecule to become an induced subgraph of the \a Father, i.e. for every bond in Father
- * with end points present in this molecule, bond is created in this molecule.
- * Special care was taken to ensure that this is of complexity O(N), where N is the \a Father's molecule::AtomCount.
- * \param *out output stream for debugging
- * \param *Father father molecule
- * \return true - is induced subgraph, false - there are atoms with fathers not in \a Father
- * \todo not checked, not fully working probably
- */
-bool molecule::BuildInducedSubgraph(ofstream *out, const molecule *Father)
-{
-  atom *Walker = NULL, *OtherAtom = NULL;
-  bool status = true;
-  atom **ParentList = Malloc<atom*>(Father->AtomCount, "molecule::BuildInducedSubgraph: **ParentList");
-
-  *out << Verbose(2) << "Begin of BuildInducedSubgraph." << endl;
-
+  BreadthFirstSearchAdd_Free(BFS);
+}
+;
+
+/** Adds a bond as a copy to a given one
+ * \param *left leftatom of new bond
+ * \param *right rightatom of new bond
+ * \param *CopyBond rest of fields in bond are copied from this
+ * \return pointer to new bond
+ */
+bond * molecule::CopyBond(atom *left, atom *right, bond *CopyBond)
+{
+  bond *Binder = AddBond(left, right, CopyBond->BondDegree);
+  Binder->Cyclic = CopyBond->Cyclic;
+  Binder->Type = CopyBond->Type;
+  return Binder;
+}
+;
+
+void BuildInducedSubgraph_Init(ofstream *out, atom **&ParentList, int AtomCount)
+{
   // reset parent list
+  ParentList = Calloc<atom*> (AtomCount, "molecule::BuildInducedSubgraph_Init: **ParentList");
   *out << Verbose(3) << "Resetting ParentList." << endl;
-  for (int i=Father->AtomCount;i--;)
-    ParentList[i] = NULL;
-
+}
+;
+
+void BuildInducedSubgraph_FillParentList(ofstream *out, const molecule *mol, const molecule *Father, atom **&ParentList)
+{
   // fill parent list with sons
   *out << Verbose(3) << "Filling Parent List." << endl;
-  Walker = start;
-  while (Walker->next != end) {
+  atom *Walker = mol->start;
+  while (Walker->next != mol->end) {
     Walker = Walker->next;
     ParentList[Walker->father->nr] = Walker;
     // Outputting List for debugging
-    *out << Verbose(4) << "Son["<< Walker->father->nr <<"] of " << Walker->father <<  " is " << ParentList[Walker->father->nr] << "." << endl;
-  }
-
+    *out << Verbose(4) << "Son[" << Walker->father->nr << "] of " << Walker->father << " is " << ParentList[Walker->father->nr] << "." << endl;
+  }
+
+}
+;
+
+void BuildInducedSubgraph_Finalize(ofstream *out, atom **&ParentList)
+{
+  Free(&ParentList);
+}
+;
+
+bool BuildInducedSubgraph_CreateBondsFromParent(ofstream *out, molecule *mol, const molecule *Father, atom **&ParentList)
+{
+  bool status = true;
+  atom *Walker = NULL;
+  atom *OtherAtom = NULL;
   // check each entry of parent list and if ok (one-to-and-onto matching) create bonds
   *out << Verbose(3) << "Creating bonds." << endl;
@@ -1139 +1348 @@
         status = false;
       } else {
-        for (int i=0;i<Father->NumberOfBondsPerAtom[Walker->nr];i++) {
-          OtherAtom = Father->ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
+        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+          OtherAtom = (*Runner)->GetOtherAtom(Walker);
           if (ParentList[OtherAtom->nr] != NULL) { // if otheratom is also a father of an atom on this molecule, create the bond
-            *out << Verbose(4) << "Endpoints of Bond " << Father->ListOfBondsPerAtom[Walker->nr][i] << " are both present: " << ParentList[Walker->nr]->Name << " and " << ParentList[OtherAtom->nr]->Name << "." << endl;
-            AddBond(ParentList[Walker->nr], ParentList[OtherAtom->nr], Father->ListOfBondsPerAtom[Walker->nr][i]->BondDegree);
+            *out << Verbose(4) << "Endpoints of Bond " << (*Runner) << " are both present: " << ParentList[Walker->nr]->Name << " and " << ParentList[OtherAtom->nr]->Name << "." << endl;
+            mol->AddBond(ParentList[Walker->nr], ParentList[OtherAtom->nr], (*Runner)->BondDegree);
           }
         }
@@ -1149 +1358 @@
     }
   }
-
-  Free(&ParentList);
+  return status;
+}
+;
+
+/** Adds bond structure to this molecule from \a Father molecule.
+ * This basically causes this molecule to become an induced subgraph of the \a Father, i.e. for every bond in Father
+ * with end points present in this molecule, bond is created in this molecule.
+ * Special care was taken to ensure that this is of complexity O(N), where N is the \a Father's molecule::AtomCount.
+ * \param *out output stream for debugging
+ * \param *Father father molecule
+ * \return true - is induced subgraph, false - there are atoms with fathers not in \a Father
+ * \todo not checked, not fully working probably
+ */
+bool molecule::BuildInducedSubgraph(ofstream *out, const molecule *Father)
+{
+  bool status = true;
+  atom **ParentList = NULL;
+
+  *out << Verbose(2) << "Begin of BuildInducedSubgraph." << endl;
+  BuildInducedSubgraph_Init(out, ParentList, Father->AtomCount);
+  BuildInducedSubgraph_FillParentList(out, this, Father, ParentList);
+  status = BuildInducedSubgraph_CreateBondsFromParent(out, this, Father, ParentList);
+  BuildInducedSubgraph_Finalize(out, ParentList);
   *out << Verbose(2) << "End of BuildInducedSubgraph." << endl;
   return status;
-};
-
+}
+;
 
 /** For a given keyset \a *Fragment, checks whether it is connected in the current molecule.
@@ -1173 +1403 @@
   // count number of atoms in graph
   size = 0;
-  for(KeySet::iterator runner = Fragment->begin(); runner != Fragment->end(); runner++)
+  for (KeySet::iterator runner = Fragment->begin(); runner != Fragment->end(); runner++)
     size++;
   if (size > 1)
-    for(KeySet::iterator runner = Fragment->begin(); runner != Fragment->end(); runner++) {
+    for (KeySet::iterator runner = Fragment->begin(); runner != Fragment->end(); runner++) {
       Walker = FindAtom(*runner);
       BondStatus = false;
-      for(KeySet::iterator runners = Fragment->begin(); runners != Fragment->end(); runners++) {
+      for (KeySet::iterator runners = Fragment->begin(); runners != Fragment->end(); runners++) {
         Walker2 = FindAtom(*runners);
-        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
-          if (ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker) == Walker2) {
+        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+          if ((*Runner)->GetOtherAtom(Walker) == Walker2) {
             BondStatus = true;
             break;