source: src/molecule_graph.cpp@ 07051c

/*
 * molecule_graph.cpp
 *
 *  Created on: Oct 5, 2009
 *      Author: heber
 */

#include "atom.hpp"
#include "bond.hpp"
#include "config.hpp"
#include "element.hpp"
#include "helpers.hpp"
#include "lists.hpp"
#include "memoryallocator.hpp"
#include "molecule.hpp"

/************************************* 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::CreateAdjacencyListFromDbondFile(ofstream *out, ifstream *input)
{

  // 1 We will parse bonds out of the dbond file created by tremolo.
  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.
  }

  CreateListOfBondsPerAtom(out);
};


/** Creates an adjacency list of the molecule.
 * Generally, we use the CSD approach to bond recognition, that is the the distance
 * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
 * a threshold t = 0.4 Angstroem.
 * To make it O(N log N) the function uses the linked-cell technique as follows:
 * The procedure is step-wise:
 *  -# Remove every bond in list
 *  -# Count the atoms in the molecule with CountAtoms()
 *  -# partition cell into smaller linked cells of size \a bonddistance
 *  -# 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
 * \param *out out stream for printing the matrix, NULL if no output
 * \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;

  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);
  }

  // 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
    Walker = start;
    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;


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



    // 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;

    // 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;
  } else
    *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
  *out << Verbose(0) << "End of CreateAdjacencyList." << endl;
  Free(&matrix);

};

/** Counts all cyclic bonds and returns their number.
 * \note Hydrogen bonds can never by cyclic, thus no check for that
 * \param *out output stream for debugging
 * \return number opf cyclic bonds
 */
int molecule::CountCyclicBonds(ofstream *out)
{
  int No = 0;
  int *MinimumRingSize = NULL;
  MoleculeLeafClass *Subgraphs = NULL;
  class StackClass<bond *> *BackEdgeStack = NULL;
  bond *Binder = first;
  if ((Binder->next != last) && (Binder->next->Type == Undetermined)) {
    *out << Verbose(0) << "No Depth-First-Search analysis performed so far, calling ..." << endl;
    Subgraphs = DepthFirstSearchAnalysis(out, BackEdgeStack);
    while (Subgraphs->next != NULL) {
      Subgraphs = Subgraphs->next;
      delete(Subgraphs->previous);
    }
    delete(Subgraphs);
    delete[](MinimumRingSize);
  }
  while(Binder->next != last) {
    Binder = Binder->next;
    if (Binder->Cyclic)
      No++;
  }
  delete(BackEdgeStack);
  return No;
};
/** Returns Shading as a char string.
 * \param color the Shading
 * \return string of the flag
 */
string molecule::GetColor(enum Shading color)
{
  switch(color) {
    case white:
      return "white";
      break;
    case lightgray:
      return "lightgray";
      break;
    case darkgray:
      return "darkgray";
      break;
    case black:
      return "black";
      break;
    default:
      return "uncolored";
      break;
  };
};

/** Performs a Depth-First search on this molecule.
 * Marks bonds in molecule as cyclic, bridge, ... and atoms as
 * articulations points, ...
 * We use the algorithm from [Even, Graph Algorithms, p.62].
 * \param *out output stream for debugging
 * \param *&BackEdgeStack NULL pointer to StackClass with all the found back edges, allocated and filled on return
 * \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);
  BackEdgeStack = new StackClass<bond *> (BondCount);
  MoleculeLeafClass *SubGraphs = new MoleculeLeafClass(NULL);
  MoleculeLeafClass *LeafWalker = SubGraphs;
  int CurrentGraphNr = 0, OldGraphNr;
  int ComponentNumber = 0;
  atom *Walker = NULL, *OtherAtom = NULL, *Root = start->next;
  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();

    // 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;
    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)
        if (Binder == NULL) {
          *out << Verbose(2) << "No more Unused Bonds." << endl;
          break;
        } else
          Binder = NULL;
      } 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))
        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

    // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph
    *out << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << CurrentGraphNr << "." << endl;
    LeafWalker->Leaf->Output(out);
    *out << endl;

    // step on to next root
    while ((Root != end) && (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;
    }
  }
  // set cyclic bond criterium on "same LP" basis
  Binder = first;
  while(Binder->next != last) {
    Binder = Binder->next;
    if (Binder->rightatom->LowpointNr == Binder->leftatom->LowpointNr) { // cyclic ??
      Binder->Cyclic = true;
      NoCyclicBonds++;
    }
  }


  *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;
  }

  *out << Verbose(1) << "Final graph info for each bond is:" << endl;
  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);
    *out << " ===  ";
    *out << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp.";
    OutputComponentNumber(out, Binder->rightatom);
    *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;
};

/** Analyses the cycles found and returns minimum of all cycle lengths.
 * We begin with a list of Back edges found during DepthFirstSearchAnalysis(). We go through this list - one end is the Root,
 * the other our initial Walker - and do a Breadth First Search for the Root. We mark down each Predecessor and as soon as
 * we have found the Root via BFS, we may climb back the closed cycle via the Predecessors. Thereby we mark atoms and bonds
 * as cyclic and print out the cycles.
 * \param *out output stream for debugging
 * \param *BackEdgeStack stack with all back edges found during DFS scan. Beware: This stack contains the bonds from the total molecule, not from the subgraph!
 * \param *&MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found, if set is maximum search distance
 * \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);

  *out << Verbose(1) << "Analysing cycles ... " << endl;
  NumCycles = 0;
  while (!BackEdgeStack->IsEmpty()) {
    BackEdge = BackEdgeStack->PopFirst();
    // this is the target
    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);
    *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);
};

/** Sets the next component number.
 * This is O(N) as the number of bonds per atom is bound.
 * \param *vertex atom whose next atom::*ComponentNr is to be set
 * \param nr number to use
 */
void molecule::SetNextComponentNumber(atom *vertex, int nr)
{
  int i=0;
  if (vertex != NULL) {
    for(;i<NumberOfBondsPerAtom[vertex->nr];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])
      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;
  }
};

/** Returns next unused bond for this atom \a *vertex or NULL of none exists.
 * \param *vertex atom to regard
 * \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]);
  return NULL;
};

/** Resets bond::Used flag of all bonds in this molecule.
 * \return true - success, false - -failure
 */
void molecule::ResetAllBondsToUnused()
{
  bond *Binder = first;
  while (Binder->next != last) {
    Binder = Binder->next;
    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.
 * \param *out output stream for debugging
 * \param *list
 */
void OutputAlreadyVisited(ofstream *out, int *list)
{
  *out << Verbose(4) << "Already Visited Bonds:\t";
  for(int i=1;i<=list[0];i++) *out << Verbose(0) << list[i] << "  ";
  *out << endl;
};


/** Storing the bond structure of a molecule to file.
 * Simply stores Atom::nr and then the Atom::nr of all bond partners per line.
 * \param *out output stream for debugging
 * \param *path path to file
 * \return true - file written successfully, false - writing failed
 */
bool molecule::StoreAdjacencyToFile(ofstream *out, char *path)
{
  ofstream AdjacencyFile;
  atom *Walker = NULL;
  stringstream line;
  bool status = true;

  line << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
  AdjacencyFile.open(line.str().c_str(), ios::out);
  *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;
    }
    AdjacencyFile.close();
    *out << Verbose(1) << "done." << endl;
  } else {
    *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
    status = false;
  }

  return status;
};

/** Checks contents of adjacency file against bond structure in structure molecule.
 * \param *out output stream for debugging
 * \param *path path to file
 * \param **ListOfAtoms allocated (molecule::AtomCount) and filled lookup table for ids (Atom::nr) to *Atom
 * \return true - structure is equal, false - not equivalence
 */
bool molecule::CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms)
{
  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 {
    *out << Verbose(1) << "Adjacency file not found." << endl;
    status = false;
  }
  *out << endl;
  Free(&buffer);

  return status;
};


/** Picks from a global stack with all back edges the ones in the fragment.
 * \param *out output stream for debugging
 * \param **ListOfLocalAtoms array of father atom::nr to local atom::nr (reverse of atom::father)
 * \param *ReferenceStack stack with all the back egdes
 * \param *LocalStack stack to be filled
 * \return true - everything ok, false - ReferenceStack was empty
 */
bool molecule::PickLocalBackEdges(ofstream *out, atom **ListOfLocalAtoms, class StackClass<bond *> *&ReferenceStack, class StackClass<bond *> *&LocalStack)
{
  bool status = true;
  if (ReferenceStack->IsEmpty()) {
    cerr << "ReferenceStack is empty!" << endl;
    return false;
  }
  bond *Binder = ReferenceStack->PopFirst();
  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
    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;
          break;
        }
      }
    Binder = ReferenceStack->PopFirst();  // loop the stack for next item
    *out << Verbose(3) << "Current candidate edge " << Binder << "." << endl;
    ReferenceStack->Push(Binder);
  } while (FirstBond != Binder);

  return status;
};


/** Adds atoms up to \a BondCount distance from \a *Root and notes them down in \a **AddedAtomList.
 * Gray vertices are always enqueued in an StackClass<atom *> FIFO queue, the rest is usual BFS with adding vertices found was
 * white and putting into queue.
 * \param *out output stream for debugging
 * \param *Mol Molecule class to add atoms to
 * \param **AddedAtomList list with added atom pointers, index is atom father's number
 * \param **AddedBondList list with added bond pointers, index is bond father's number
 * \param *Root root vertex for BFS
 * \param *Bond bond not to look beyond
 * \param BondOrder maximum distance for vertices to add
 * \param IsAngstroem lengths are in angstroem or bohrradii
 */
void 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);
  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
    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;

  // and go on ... Queue always contains all lightgray vertices
  while (!AtomStack->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
              }
            }
          }
        } 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
            }
          }
        }
      }
    }
    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;

  // reset parent list
  *out << Verbose(3) << "Resetting ParentList." << endl;
  for (int i=Father->AtomCount;i--;)
    ParentList[i] = NULL;

  // fill parent list with sons
  *out << Verbose(3) << "Filling Parent List." << endl;
  Walker = start;
  while (Walker->next != 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;
  }

  // check each entry of parent list and if ok (one-to-and-onto matching) create bonds
  *out << Verbose(3) << "Creating bonds." << endl;
  Walker = Father->start;
  while (Walker->next != Father->end) {
    Walker = Walker->next;
    if (ParentList[Walker->nr] != NULL) {
      if (ParentList[Walker->nr]->father != Walker) {
        status = false;
      } else {
        for (int i=0;i<Father->NumberOfBondsPerAtom[Walker->nr];i++) {
          OtherAtom = Father->ListOfBondsPerAtom[Walker->nr][i]->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);
          }
        }
      }
    }
  }

  Free(&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.
 * \param *out output stream for debugging
 * \param *Fragment Keyset of fragment's vertices
 * \return true - connected, false - disconnected
 * \note this is O(n^2) for it's just a bug checker not meant for permanent use!
 */
bool molecule::CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment)
{
  atom *Walker = NULL, *Walker2 = NULL;
  bool BondStatus = false;
  int size;

  *out << Verbose(1) << "Begin of CheckForConnectedSubgraph" << endl;
  *out << Verbose(2) << "Disconnected atom: ";

  // count number of atoms in graph
  size = 0;
  for(KeySet::iterator runner = Fragment->begin(); runner != Fragment->end(); runner++)
    size++;
  if (size > 1)
    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++) {
        Walker2 = FindAtom(*runners);
        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
          if (ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker) == Walker2) {
            BondStatus = true;
            break;
          }
          if (BondStatus)
            break;
        }
      }
      if (!BondStatus) {
        *out << (*Walker) << endl;
        return false;
      }
    }
  else {
    *out << "none." << endl;
    return true;
  }
  *out << "none." << endl;

  *out << Verbose(1) << "End of CheckForConnectedSubgraph" << endl;

  return true;
}