source: src/molecule_graph.cpp@ ce7cc5

/*
 * 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 "linkedcell.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.
  }
};


/** 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
 *  -# 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;
  atom *OtherWalker = NULL;
  atom **AtomMap = NULL;
  int n[NDIM];
  double distance, MinDistance, MaxDistance;
  LinkedCell *LC = NULL;
  LinkedNodes *List = NULL;
  LinkedNodes *OtherList = NULL;

  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) {
    LC = new LinkedCell(this, bonddistance);

    // create a list to map Tesselpoint::nr to atom *
    AtomMap = Malloc<atom *>(AtomCount, "molecule::CreateAdjacencyList - **AtomCount");
    Walker = start;
    while (Walker->next != end) {
      Walker = Walker->next;
      AtomMap[Walker->nr] = Walker;
    }

    // 3a. go through every cell
    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]++) {
          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::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]++) {
                    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::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;
                          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;
                          }
                        }
                      }
                    }
                  }
            }
          }
        }
    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
    CorrectBondDegree(out);

    // output bonds for debugging (if bond chain list was correctly installed)
    ActOnAllAtoms( &atom::OutputBondOfAtom, out );
  } else
    *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
  *out << Verbose(0) << "End of CreateAdjacencyList." << endl;
};

/** Prints a list of all bonds to \a *out.
 * \param output stream
 */
void molecule::OutputBondsList(ofstream *out)
{
  *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)
{
  int No = 0;

  if (BondCount != 0) {
    *out << Verbose(1) << "Correcting Bond degree of each bond ... " << endl;
    do {
      No = SumPerAtom( &atom::CorrectBondDegree, out );
    } while (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.
 * \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)
{
  NoCyclicBonds = 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)
      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) {
    case white:
      return "white";
      break;
    case lightgray:
      return "lightgray";
      break;
    case darkgray:
      return "darkgray";
      break;
    case black:
      return "black";
      break;
    default:
      return "uncolored";
      break;
  };
};

void SetWalkersGraphNr(ofstream *out, bool &BackStepping, atom *&Walker, int &CurrentGraphNr, class StackClass<atom *> *&AtomStack)
{
  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++;
  }
};

void ProbeAlongUnusedBond(ofstream *out, molecule *mol, bond *&Binder, bool &BackStepping, atom *&Walker, class StackClass<bond *> *&BackEdgeStack)
{
  atom *OtherAtom = NULL;

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

void CheckForaNewComponent(ofstream *out, molecule *mol, bool &BackStepping, atom *&Walker, atom *&Root, int &ComponentNumber, class StackClass<atom *> *&AtomStack, 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 != 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, ComponentNumber);
      *out << Verbose(3) << "(7) Walker[" << Walker->Name << "]'s Ancestor's Compont is " << ComponentNumber << "." << endl;
      mol->SetNextComponentNumber(Walker, ComponentNumber);
      *out << Verbose(3) << "(7) Walker[" << Walker->Name << "]'s Compont is " << ComponentNumber << "." << endl;
      do {
        OtherAtom = AtomStack->PopLast();
        LeafWalker->Leaf->AddCopyAtom(OtherAtom);
        mol->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;
  }
};

void CleanRootStackDownTillWalker(ofstream *out, molecule *mol, bool &BackStepping, atom *&Root, atom *&Walker, bond *&Binder, int &ComponentNumber, class StackClass<atom *> *&AtomStack, MoleculeLeafClass *&LeafWalker)
{
  atom *OtherAtom = NULL;

  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);
    mol->SetNextComponentNumber(Root, ComponentNumber);
    *out << Verbose(3) << "(9) Root[" << Root->Name << "]'s Component is " << ComponentNumber << "." << endl;
    mol->SetNextComponentNumber(Walker, ComponentNumber);
    *out << Verbose(3) << "(9) Walker[" << Walker->Name << "]'s Component is " << ComponentNumber << "." << endl;
    do {
      OtherAtom = AtomStack->PopLast();
      LeafWalker->Leaf->AddCopyAtom(OtherAtom);
      mol->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 = mol->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;
    }
  }
};


/** 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;
  atom *Root = start->next;
  bond *Binder = NULL;
  bool BackStepping = false;

  *out << Verbose(0) << "Begin of DepthFirstSearchAnalysis" << endl;

  ResetAllBondsToUnused();
  SetAtomValueToValue( -1, &atom::GraphNr );
  ActOnAllAtoms( &atom::InitComponentNr );
  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
        SetWalkersGraphNr(out, BackStepping, Walker, CurrentGraphNr, AtomStack);

        ProbeAlongUnusedBond(out, this, Binder, BackStepping, Walker, BackEdgeStack);

        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;

      CheckForaNewComponent(out, this, BackStepping, Walker, Root,ComponentNumber, AtomStack, LeafWalker );

      CleanRootStackDownTillWalker(out, this, BackStepping, Root, Walker, Binder, ComponentNumber, AtomStack, LeafWalker);

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

  OutputGraphInfoPerAtom(out);

  OutputGraphInfoPerBond(out);

  // free all and exit
  delete(AtomStack);
  *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()
{
  NoCyclicBonds = 0;
  bond *Binder = first;
  while(Binder->next != last) {
    Binder = Binder->next;
    if (Binder->rightatom->LowpointNr == Binder->leftatom->LowpointNr) { // cyclic ??
      Binder->Cyclic = true;
      NoCyclicBonds++;
    }
  }
};

/** Output graph information per atom.
 * \param *out output stream
 */
void molecule::OutputGraphInfoPerAtom(ofstream *out)
{
  *out << Verbose(1) << "Final graph info for each atom is:" << endl;
  ActOnAllAtoms( &atom::OutputGraphInfo, out );
};

/** Output graph information per bond.
 * \param *out output stream
 */
void molecule::OutputGraphInfoPerBond(ofstream *out)
{
  *out << Verbose(1) << "Final graph info for each bond is:" << endl;
  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.";
    Binder->leftatom->OutputComponentNumber(out);
    *out << " ===  ";
    *out << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp.";
    Binder->rightatom->OutputComponentNumber(out);
    *out << ">." << endl;
    if (Binder->Cyclic) // cyclic ??
      *out << Verbose(3) << "Lowpoint at each side are equal: CYCLIC!" << endl;
  }
};

/**  initialise each vertex as white with no predecessor, empty queue, color Root lightgray.
 *
 */
void InitializeAccounting(atom **&PredecessorList, int *&ShortestPathList, enum Shading *&ColorList, int AtomCount)
{
  for (int i=AtomCount;i--;) {
    PredecessorList[i] = NULL;
    ShortestPathList[i] = -1;
    ColorList[i] = white;
  }
};

void ResetAccounting(atom *&Walker, class StackClass<atom *> *&TouchedStack, int *&ShortestPathList, class StackClass<atom *> *&BFSStack)
{
  ShortestPathList[Walker->nr] = 0;
  BFSStack->ClearStack();  // start with empty BFS stack
  BFSStack->Push(Walker);
  TouchedStack->Push(Walker);
};

void CyclicBFSFromRootToRoot(ofstream *out, atom *&Root, bond *&BackEdge, class StackClass<atom *> *&TouchedStack, int *&ShortestPathList, atom **&PredecessorList, class StackClass<atom *> *&BFSStack, enum Shading *&ColorList)
{
  atom *Walker = NULL;
  atom *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 (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 (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 " << *(*Runner) << " 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])));
};

void RetrieveCycleMembers(ofstream *out, atom *&Root, atom *&OtherAtom, bond *&BackEdge, atom **&PredecessorList, int *&MinimumRingSize, int &MinRingSize)
{
  atom *Walker = NULL;
  int NumCycles = 0;
  int RingSize = -1;

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

void CleanAccounting(class StackClass<atom *> *&TouchedStack, atom **&PredecessorList, int *&ShortestPathList, enum Shading *&ColorList)
{
  atom *Walker = NULL;
  while (!TouchedStack->IsEmpty()){
    Walker = TouchedStack->PopFirst();
    PredecessorList[Walker->nr] = NULL;
    ShortestPathList[Walker->nr] = -1;
    ColorList[Walker->nr] = white;
  }
};


void BFSToNextCycle(ofstream *out, atom *&Root, atom *&Walker, bond *&BackEdge, int *&MinimumRingSize, int AtomCount)
{
  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 *OtherAtom = Walker;

  InitializeAccounting(PredecessorList, ShortestPathList, ColorList, AtomCount);

  ResetAccounting(Walker, TouchedStack, ShortestPathList, BFSStack);
  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 (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
      if (((*Runner) != BackEdge) || (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 (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;
  }
  //CleanAccountingLists(TouchedStack, PredecessorList, ShortestPathList, ColorList);

  Free(&PredecessorList);
  Free(&ShortestPathList);
  Free(&ColorList);
  delete(BFSStack);
};

void AssignRingSizetoNonCycleMembers(ofstream *out, int *&MinimumRingSize, int &MinRingSize, int &NumCycles, molecule *mol, bond *&BackEdge)
{
  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;
        BFSToNextCycle(out, Root, Walker, BackEdge, 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.
 * 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;
  atom *OtherAtom = NULL;
  atom *Root = NULL;
  bond *BackEdge = NULL;
  int NumCycles = 0;
  int MinRingSize = -1;

  InitializeAccounting(PredecessorList, ShortestPathList, ColorList, AtomCount);

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

    ResetAccounting(Walker, TouchedStack, ShortestPathList, BFSStack);

    *out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
    OtherAtom = NULL;
    CyclicBFSFromRootToRoot(out, Root, BackEdge, TouchedStack, ShortestPathList, PredecessorList, BFSStack, ColorList);

    RetrieveCycleMembers(out, Root, OtherAtom, BackEdge, PredecessorList, MinimumRingSize, MinRingSize);

    CleanAccounting(TouchedStack, PredecessorList, ShortestPathList, ColorList);
  }
  Free(&PredecessorList);
  Free(&ShortestPathList);
  Free(&ColorList);
  delete(BFSStack);

  AssignRingSizetoNonCycleMembers(out, MinimumRingSize, MinRingSize, NumCycles, this, BackEdge);

};

/** 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)
{
  size_t i=0;
  if (vertex != NULL) {
    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 == vertex->ListOfBonds.size())
      cerr << "Error: All Component entries are already occupied!" << endl;
  } else
      cerr << "Error: Given vertex is NULL!" << endl;
};

/** 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 (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()
{
  bond *Binder = first;
  while (Binder->next != last) {
    Binder = Binder->next;
    Binder->ResetUsed();
  }
};

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

  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 = Malloc<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.
 * \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;
  bool status = true;
  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;
    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.
 * \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 (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
    *out << Verbose(3) << "Current candidate edge " << Binder << "." << endl;
    ReferenceStack->Push(Binder);
  } while (FirstBond != Binder);

  return status;
};

struct BFSAccounting {
  atom **PredecessorList;
  int *ShortestPathList;
  enum Shading *ColorList;
  class StackClass<atom *> *AtomStack;
  int AtomCount;
  int BondOrder;
  atom *Root;
};

void BreadthFirstSearchAdd_Init(struct BFSAccounting &BFS, atom *&Root, int AtomCount, int BondOrder, atom **AddedAtomList = NULL)
{
  BFS.AtomCount = AtomCount;
  BFS.BondOrder = BondOrder;
  BFS.PredecessorList = Malloc<atom*>(AtomCount, "molecule::BreadthFirstSearchAdd_Init: **PredecessorList");
  BFS.ShortestPathList = Malloc<int>(AtomCount, "molecule::BreadthFirstSearchAdd_Init: *ShortestPathList");
  BFS.ColorList = Malloc<enum Shading>(AtomCount, "molecule::BreadthFirstSearchAdd_Init: *ColorList");
  BFS.AtomStack = new StackClass<atom *>(AtomCount);

  BFS.Root = Root;
  BFS.AtomStack->ClearStack();
  BFS.AtomStack->Push(Root);

  // initialise each vertex as white with no predecessor, empty queue, color Root lightgray
  for (int i=AtomCount;i--;) {
    BFS.PredecessorList[i] = NULL;
    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;
};

void BreadthFirstSearchAdd_Free(struct BFSAccounting &BFS)
{
  Free(&BFS.PredecessorList);
  Free(&BFS.ShortestPathList);
  Free(&BFS.ColorList);
  delete(BFS.AtomStack);
  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.AtomStack->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.
 * 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)
{
  struct BFSAccounting BFS;
  atom *Walker = NULL, *OtherAtom = NULL;
  bond *Binder = NULL;

  // add Root if not done yet
  if (AddedAtomList[Root->nr] == NULL)  // add Root if not yet present
    AddedAtomList[Root->nr] = Mol->AddCopyAtom(Root);

  BreadthFirstSearchAdd_Init(BFS, Root, BondOrder, AtomCount, AddedAtomList);

  // and go on ... Queue always contains all lightgray vertices
  while (!BFS.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 = BFS.AtomStack->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 {
          BreadthFirstSearchAdd_VisitedNode(out, Mol, BFS, Walker, OtherAtom, Binder, Bond, AddedAtomList, AddedBondList, IsAngstroem);
        }
      }
    }
    BFS.ColorList[Walker->nr] = black;
    *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << 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;
};


/** 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 (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 " << (*Runner) << " are both present: " << ParentList[Walker->nr]->Name << " and " << ParentList[OtherAtom->nr]->Name << "." << endl;
            AddBond(ParentList[Walker->nr], ParentList[OtherAtom->nr], (*Runner)->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 (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
          if ((*Runner)->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;
}