source: src/molecule_fragmentation.cpp@ 5034e1

/*
 * molecule_fragmentation.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"
#include "periodentafel.hpp"

/************************************* Functions for class molecule *********************************/


/** Estimates by educated guessing (using upper limit) the expected number of fragments.
 * The upper limit is
 * \f[
 *  n = N \cdot C^k
 * \f]
 * where \f$C=2^c\f$ and c is the maximum bond degree over N number of atoms.
 * \param *out output stream for debugging
 * \param order bond order k
 * \return number n of fragments
 */
int molecule::GuesstimateFragmentCount(ofstream *out, int order)
{
  int c = 0;
  int FragmentCount;
  // get maximum bond degree
  atom *Walker = start;
  while (Walker->next != end) {
    Walker = Walker->next;
    c = (NumberOfBondsPerAtom[Walker->nr] > c) ? NumberOfBondsPerAtom[Walker->nr] : c;
  }
  FragmentCount = NoNonHydrogen*(1 << (c*order));
  *out << Verbose(1) << "Upper limit for this subgraph is " << FragmentCount << " for " << NoNonHydrogen << " non-H atoms with maximum bond degree of " << c << "." << endl;
  return FragmentCount;
};

/** Scans a single line for number and puts them into \a KeySet.
 * \param *out output stream for debugging
 * \param *buffer buffer to scan
 * \param &CurrentSet filled KeySet on return
 * \return true - at least one valid atom id parsed, false - CurrentSet is empty
 */
bool ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet)
{
  stringstream line;
  int AtomNr;
  int status = 0;

  line.str(buffer);
  while (!line.eof()) {
    line >> AtomNr;
    if (AtomNr >= 0) {
      CurrentSet.insert(AtomNr);  // insert at end, hence in same order as in file!
      status++;
    } // else it's "-1" or else and thus must not be added
  }
  *out << Verbose(1) << "The scanned KeySet is ";
  for(KeySet::iterator runner = CurrentSet.begin(); runner != CurrentSet.end(); runner++) {
    *out << (*runner) << "\t";
  }
  *out << endl;
  return (status != 0);
};

/** Parses the KeySet file and fills \a *FragmentList from the known molecule structure.
 * Does two-pass scanning:
 * -# Scans the keyset file and initialises a temporary graph
 * -# Scans TEFactors file and sets the TEFactor of each key set in the temporary graph accordingly
 * Finally, the temporary graph is inserted into the given \a FragmentList for return.
 * \param *out output stream for debugging
 * \param *path path to file
 * \param *FragmentList empty, filled on return
 * \return true - parsing successfully, false - failure on parsing (FragmentList will be NULL)
 */
bool ParseKeySetFile(ofstream *out, char *path, Graph *&FragmentList)
{
  bool status = true;
  ifstream InputFile;
  stringstream line;
  GraphTestPair testGraphInsert;
  int NumberOfFragments = 0;
  char *filename = Malloc<char>(MAXSTRINGSIZE, "molecule::ParseKeySetFile - filename");

  if (FragmentList == NULL) { // check list pointer
    FragmentList = new Graph;
  }

  // 1st pass: open file and read
  *out << Verbose(1) << "Parsing the KeySet file ... " << endl;
  sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, KEYSETFILE);
  InputFile.open(filename);
  if (InputFile != NULL) {
    // each line represents a new fragment
    char *buffer = Malloc<char>(MAXSTRINGSIZE, "molecule::ParseKeySetFile - *buffer");
    // 1. parse keysets and insert into temp. graph
    while (!InputFile.eof()) {
      InputFile.getline(buffer, MAXSTRINGSIZE);
      KeySet CurrentSet;
      if ((strlen(buffer) > 0) && (ScanBufferIntoKeySet(out, buffer, CurrentSet))) {  // if at least one valid atom was added, write config
        testGraphInsert = FragmentList->insert(GraphPair (CurrentSet,pair<int,double>(NumberOfFragments++,1)));  // store fragment number and current factor
        if (!testGraphInsert.second) {
          cerr << "KeySet file must be corrupt as there are two equal key sets therein!" << endl;
        }
      }
    }
    // 2. Free and done
    InputFile.close();
    InputFile.clear();
    Free(&buffer);
    *out << Verbose(1) << "done." << endl;
  } else {
    *out << Verbose(1) << "File " << filename << " not found." << endl;
    status = false;
  }

  return status;
};

/** Parses the TE factors file and fills \a *FragmentList from the known molecule structure.
 * -# Scans TEFactors file and sets the TEFactor of each key set in the temporary graph accordingly
 * \param *out output stream for debugging
 * \param *path path to file
 * \param *FragmentList graph whose nodes's TE factors are set on return
 * \return true - parsing successfully, false - failure on parsing
 */
bool ParseTEFactorsFile(ofstream *out, char *path, Graph *FragmentList)
{
  bool status = true;
  ifstream InputFile;
  stringstream line;
  GraphTestPair testGraphInsert;
  int NumberOfFragments = 0;
  double TEFactor;
  char *filename = Malloc<char>(MAXSTRINGSIZE, "molecule::ParseTEFactorsFile - filename");

  if (FragmentList == NULL) { // check list pointer
    FragmentList = new Graph;
  }

  // 2nd pass: open TEFactors file and read
  *out << Verbose(1) << "Parsing the TEFactors file ... " << endl;
  sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, TEFACTORSFILE);
  InputFile.open(filename);
  if (InputFile != NULL) {
    // 3. add found TEFactors to each keyset
    NumberOfFragments = 0;
    for(Graph::iterator runner = FragmentList->begin();runner != FragmentList->end(); runner++) {
      if (!InputFile.eof()) {
        InputFile >> TEFactor;
        (*runner).second.second = TEFactor;
        *out << Verbose(2) << "Setting " << ++NumberOfFragments << " fragment's TEFactor to " << (*runner).second.second << "." << endl;
      } else {
        status = false;
        break;
      }
    }
    // 4. Free and done
    InputFile.close();
    *out << Verbose(1) << "done." << endl;
  } else {
    *out << Verbose(1) << "File " << filename << " not found." << endl;
    status = false;
  }

  // free memory
  Free(&filename);

  return status;
};

/** Stores key sets to file.
 * \param *out output stream for debugging
 * \param KeySetList Graph with Keysets
 * \param *path path to file
 * \return true - file written successfully, false - writing failed
 */
bool StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path)
{
  ofstream output;
  bool status =  true;
  string line;

  // open KeySet file
  line = path;
  line.append("/");
  line += FRAGMENTPREFIX;
  line += KEYSETFILE;
  output.open(line.c_str(), ios::out);
  *out << Verbose(1) << "Saving key sets of the total graph ... ";
  if(output != NULL) {
    for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++) {
      for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
        if (sprinter != (*runner).first.begin())
          output << "\t";
        output << *sprinter;
      }
      output << endl;
    }
    *out << "done." << endl;
  } else {
    cerr << "Unable to open " << line << " for writing keysets!" << endl;
    status = false;
  }
  output.close();
  output.clear();

  return status;
};


/** Stores TEFactors to file.
 * \param *out output stream for debugging
 * \param KeySetList Graph with factors
 * \param *path path to file
 * \return true - file written successfully, false - writing failed
 */
bool StoreTEFactorsFile(ofstream *out, Graph &KeySetList, char *path)
{
  ofstream output;
  bool status =  true;
  string line;

  // open TEFactors file
  line = path;
  line.append("/");
  line += FRAGMENTPREFIX;
  line += TEFACTORSFILE;
  output.open(line.c_str(), ios::out);
  *out << Verbose(1) << "Saving TEFactors of the total graph ... ";
  if(output != NULL) {
    for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++)
      output << (*runner).second.second << endl;
    *out << Verbose(1) << "done." << endl;
  } else {
    *out << Verbose(1) << "failed to open " << line << "." << endl;
    status = false;
  }
  output.close();

  return status;
};

/** For a given graph, sorts KeySets into a (index, keyset) map.
 * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
 * \return map from index to keyset
 */
map<int,KeySet> * GraphToIndexedKeySet(Graph *GlobalKeySetList)
{
  map<int,KeySet> *IndexKeySetList = new map<int,KeySet>;
  for(Graph::iterator runner = GlobalKeySetList->begin(); runner != GlobalKeySetList->end(); runner++) {
    IndexKeySetList->insert( pair<int,KeySet>(runner->second.first,runner->first) );
  }
  return IndexKeySetList;
};

/** Inserts a (\a No, \a value) pair into the list, overwriting present one.
 * Note if values are equal, No will decided on which is first
 * \param *out output stream for debugging
 * \param &AdaptiveCriteriaList list to insert into
 * \param &IndexedKeySetList list to find key set for a given index \a No
 * \param FragOrder current bond order of fragment
 * \param No index of keyset
 * \param value energy value
 */
void InsertIntoAdaptiveCriteriaList(ofstream *out, map<int, pair<double,int> > *AdaptiveCriteriaList, map<int,KeySet> &IndexKeySetList, int FragOrder, int No, double Value)
{
  map<int,KeySet>::iterator marker = IndexKeySetList.find(No);    // find keyset to Frag No.
  if (marker != IndexKeySetList.end()) {  // if found
    Value *= 1 + MYEPSILON*(*((*marker).second.begin()));     // in case of equal energies this makes them not equal without changing anything actually
    // as the smallest number in each set has always been the root (we use global id to keep the doubles away), seek smallest and insert into AtomMask
    pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList->insert( make_pair(*((*marker).second.begin()), pair<double,int>( fabs(Value), FragOrder) ));
    map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
    if (!InsertedElement.second) { // this root is already present
      if ((*PresentItem).second.second < FragOrder)  // if order there is lower, update entry with higher-order term
        //if ((*PresentItem).second.first < (*runner).first)    // as higher-order terms are not always better, we skip this part (which would always include this site into adaptive increase)
        {  // if value is smaller, update value and order
        (*PresentItem).second.first = fabs(Value);
        (*PresentItem).second.second = FragOrder;
        *out << Verbose(2) << "Updated element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
      } else {
        *out << Verbose(2) << "Did not update element " <<  (*PresentItem).first << " as " << FragOrder << " is less than or equal to " << (*PresentItem).second.second << "." << endl;
      }
    } else {
      *out << Verbose(2) << "Inserted element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
    }
  } else {
    *out << Verbose(1) << "No Fragment under No. " << No << "found." << endl;
  }
};

/** Scans the adaptive order file and insert (index, value) into map.
 * \param *out output stream for debugging
 * \param *path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
 * \param &IndexedKeySetList list to find key set for a given index \a No
 * \return adaptive criteria list from file
 */
map<int, pair<double,int> > * ScanAdaptiveFileIntoMap(ofstream *out, char *path, map<int,KeySet> &IndexKeySetList)
{
  map<int, pair<double,int> > *AdaptiveCriteriaList = new map<int, pair<double,int> >;
  int No = 0, FragOrder = 0;
  double Value = 0.;
  char *buffer = Malloc<char>(MAXSTRINGSIZE, "molecule::CheckOrderAtSite: *buffer");
  sprintf(buffer, "%s/%s%s.dat", path, FRAGMENTPREFIX, ENERGYPERFRAGMENT);
  ifstream InputFile(buffer, ios::in);

  if (CountLinesinFile(InputFile) > 0) {
    // each line represents a fragment root (Atom::nr) id and its energy contribution
    InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
    InputFile.getline(buffer, MAXSTRINGSIZE);
    while(!InputFile.eof()) {
      InputFile.getline(buffer, MAXSTRINGSIZE);
      if (strlen(buffer) > 2) {
        //*out << Verbose(2) << "Scanning: " << buffer << endl;
        stringstream line(buffer);
        line >> FragOrder;
        line >> ws >> No;
        line >> ws >> Value; // skip time entry
        line >> ws >> Value;
        No -= 1;  // indices start at 1 in file, not 0
        //*out << Verbose(2) << " - yields (" << No << "," << Value << ", " << FragOrder << ")" << endl;

        // clean the list of those entries that have been superceded by higher order terms already
        InsertIntoAdaptiveCriteriaList(out, AdaptiveCriteriaList, IndexKeySetList, FragOrder, No, Value);
      }
    }
    // close and done
    InputFile.close();
    InputFile.clear();
  }
  Free(&buffer);

  return AdaptiveCriteriaList;
};

/** Maps adaptive criteria list back onto (Value, (Root Nr., Order))
 * (i.e. sorted by value to pick the highest ones)
 * \param *out output stream for debugging
 * \param &AdaptiveCriteriaList list to insert into
 * \param *mol molecule with atoms
 * \return remapped list
 */
map<double, pair<int,int> >  * ReMapAdaptiveCriteriaListToValue(ofstream *out, map<int, pair<double,int> > *AdaptiveCriteriaList, molecule *mol)
{
  atom *Walker = mol->start;
  map<double, pair<int,int> > *FinalRootCandidates = new map<double, pair<int,int> > ;
  *out << Verbose(1) << "Root candidate list is: " << endl;
  for(map<int, pair<double,int> >::iterator runner = AdaptiveCriteriaList->begin(); runner != AdaptiveCriteriaList->end(); runner++) {
    Walker = mol->FindAtom((*runner).first);
    if (Walker != NULL) {
      //if ((*runner).second.second >= Walker->AdaptiveOrder) { // only insert if this is an "active" root site for the current order
      if (!Walker->MaxOrder) {
        *out << Verbose(2) << "(" << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "])" << endl;
        FinalRootCandidates->insert( make_pair( (*runner).second.first, pair<int,int>((*runner).first, (*runner).second.second) ) );
      } else {
        *out << Verbose(2) << "Excluding (" << *Walker << ", " << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "]), as it has reached its maximum order." << endl;
      }
    } else {
      cerr << "Atom No. " << (*runner).second.first << " was not found in this molecule." << endl;
    }
  }
  return FinalRootCandidates;
};

/** Marks all candidate sites for update if below adaptive threshold.
 * Picks a given number of highest values and set *AtomMask to true.
 * \param *out output stream for debugging
 * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site, used to activate given number of site to increment order adaptively
 * \param FinalRootCandidates list candidates to check
 * \param Order desired order
 * \param *mol molecule with atoms
 * \return true - if update is necessary, false - not
 */
bool MarkUpdateCandidates(ofstream *out, bool *AtomMask, map<double, pair<int,int> > &FinalRootCandidates, int Order, molecule *mol)
{
  atom *Walker = mol->start;
  int No = -1;
  bool status = false;
  for(map<double, pair<int,int> >::iterator runner = FinalRootCandidates.upper_bound(pow(10.,Order)); runner != FinalRootCandidates.end(); runner++) {
    No = (*runner).second.first;
    Walker = mol->FindAtom(No);
    //if (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]) {
      *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl;
      AtomMask[No] = true;
      status = true;
    //} else
      //*out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", however MinimumRingSize of " << MinimumRingSize[Walker->nr] << " does not allow further adaptive increase." << endl;
  }
  return status;
};

/** print atom mask for debugging.
 * \param *out output stream for debugging
 * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site, used to activate given number of site to increment order adaptively
 * \param AtomCount number of entries in \a *AtomMask
 */
void PrintAtomMask(ofstream *out, bool *AtomMask, int AtomCount)
{
  *out << "              ";
  for(int i=0;i<AtomCount;i++)
    *out << (i % 10);
  *out << endl << "Atom mask is: ";
  for(int i=0;i<AtomCount;i++)
    *out << (AtomMask[i] ? "t" : "f");
  *out << endl;
};

/** Checks whether the OrderAtSite is still below \a Order at some site.
 * \param *out output stream for debugging
 * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site, used to activate given number of site to increment order adaptively
 * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
 * \param Order desired Order if positive, desired exponent in threshold criteria if negative (0 is single-step)
 * \param *MinimumRingSize array of max. possible order to avoid loops
 * \param *path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
 * \return true - needs further fragmentation, false - does not need fragmentation
 */
bool molecule::CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path)
{
  atom *Walker = start;
  bool status = false;

  // initialize mask list
  for(int i=AtomCount;i--;)
    AtomMask[i] = false;

  if (Order < 0) { // adaptive increase of BondOrder per site
    if (AtomMask[AtomCount] == true)  // break after one step
      return false;

    // transmorph graph keyset list into indexed KeySetList
    if (GlobalKeySetList == NULL) {
      cout << Verbose(1) << "ERROR: Given global key set list (graph) is NULL!" << endl;
      return false;
    }
    map<int,KeySet> *IndexKeySetList = GraphToIndexedKeySet(GlobalKeySetList);

    // parse the EnergyPerFragment file
    map<int, pair<double,int> > *AdaptiveCriteriaList = ScanAdaptiveFileIntoMap(out, path, *IndexKeySetList); // (Root No., (Value, Order)) !
    if (AdaptiveCriteriaList->empty()) {
      cerr << "Unable to parse file, incrementing all." << endl;
      while (Walker->next != end) {
        Walker = Walker->next;
    #ifdef ADDHYDROGEN
        if (Walker->type->Z != 1) // skip hydrogen
    #endif
        {
          AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
          status = true;
        }
      }
    }
    // then map back onto (Value, (Root Nr., Order)) (i.e. sorted by value to pick the highest ones)
    map<double, pair<int,int> > *FinalRootCandidates = ReMapAdaptiveCriteriaListToValue(out, AdaptiveCriteriaList, this);

    // pick the ones still below threshold and mark as to be adaptively updated
    MarkUpdateCandidates(out, AtomMask, *FinalRootCandidates, Order, this);

    Free(&IndexKeySetList);
    Free(&AdaptiveCriteriaList);
    Free(&FinalRootCandidates);
  } else { // global increase of Bond Order
    while (Walker->next != end) {
      Walker = Walker->next;
  #ifdef ADDHYDROGEN
      if (Walker->type->Z != 1) // skip hydrogen
  #endif
      {
        AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
        if ((Order != 0) && (Walker->AdaptiveOrder < Order)) // && (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]))
          status = true;
      }
    }
    if ((Order == 0) && (AtomMask[AtomCount] == false))  // single stepping, just check
      status = true;

    if (!status) {
      if (Order == 0)
        *out << Verbose(1) << "Single stepping done." << endl;
      else
        *out << Verbose(1) << "Order at every site is already equal or above desired order " << Order << "." << endl;
    }
  }

  PrintAtomMask(out, AtomMask, AtomCount); // for debugging

  return status;
};

/** Create a SortIndex to map from atomic labels to the sequence in which the atoms are given in the config file.
 * \param *out output stream for debugging
 * \param *&SortIndex Mapping array of size molecule::AtomCount
 * \return true - success, false - failure of SortIndex alloc
 */
bool molecule::CreateMappingLabelsToConfigSequence(ofstream *out, int *&SortIndex)
{
  if (SortIndex != NULL) {
    *out << Verbose(1) << "SortIndex is " << SortIndex << " and not NULL as expected." << endl;
    return false;
  }
  SortIndex = Malloc<int>(AtomCount, "molecule::CreateMappingLabelsToConfigSequence: *SortIndex");
  for(int i=AtomCount;i--;)
    SortIndex[i] = -1;

  int AtomNo = 0;
  SetIndexedArrayForEachAtomTo( SortIndex, &atom::nr, &IncrementalAbsoluteValue, AtomNo );

  return true;
};

/** Performs a many-body bond order analysis for a given bond order.
 * -# parses adjacency, keysets and orderatsite files
 * -# performs DFS to find connected subgraphs (to leave this in was a design decision: might be useful later)
 * -# RootStack is created for every subgraph (here, later we implement the "update 10 sites with highest energ
y contribution", and that's why this consciously not done in the following loop)
 * -# in a loop over all subgraphs
 *  -# calls FragmentBOSSANOVA with this RootStack and within the subgraph molecule structure
 *  -# creates molecule (fragment)s from the returned keysets (StoreFragmentFromKeySet)
 * -# combines the generated molecule lists from all subgraphs
 * -# saves to disk: fragment configs, adjacency, orderatsite, keyset files
 * Note that as we split "this" molecule up into a list of subgraphs, i.e. a MoleculeListClass, we have two sets
 * of vertex indices: Global always means the index in "this" molecule, whereas local refers to the molecule or
 * subgraph in the MoleculeListClass.
 * \param *out output stream for debugging
 * \param Order up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
 * \param *configuration configuration for writing config files for each fragment
 * \return 1 - continue, 2 - stop (no fragmentation occured)
 */
int molecule::FragmentMolecule(ofstream *out, int Order, config *configuration)
{
  MoleculeListClass *BondFragments = NULL;
  int *SortIndex = NULL;
  int *MinimumRingSize = new int[AtomCount];
  int FragmentCounter;
  MoleculeLeafClass *MolecularWalker = NULL;
  MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
  fstream File;
  bool FragmentationToDo = true;
  class StackClass<bond *> *BackEdgeStack = NULL, *LocalBackEdgeStack = NULL;
  bool CheckOrder = false;
  Graph **FragmentList = NULL;
  Graph *ParsedFragmentList = NULL;
  Graph TotalGraph;     // graph with all keysets however local numbers
  int TotalNumberOfKeySets = 0;
  atom **ListOfAtoms = NULL;
  atom ***ListOfLocalAtoms = NULL;
  bool *AtomMask = NULL;

  *out << endl;
#ifdef ADDHYDROGEN
  *out << Verbose(0) << "I will treat hydrogen special and saturate dangling bonds with it." << endl;
#else
  *out << Verbose(0) << "Hydrogen is treated just like the rest of the lot." << endl;
#endif

  // ++++++++++++++++++++++++++++ INITIAL STUFF: Bond structure analysis, file parsing, ... ++++++++++++++++++++++++++++++++++++++++++

  // ===== 1. Check whether bond structure is same as stored in files ====

  // fill the adjacency list
  CreateListOfBondsPerAtom(out);

  // create lookup table for Atom::nr
  FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(out, start, end, ListOfAtoms, AtomCount);

  // === compare it with adjacency file ===
  FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(out, configuration->configpath, ListOfAtoms);
  Free(&ListOfAtoms);

  // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
  Subgraphs = DepthFirstSearchAnalysis(out, BackEdgeStack);
  // fill the bond structure of the individually stored subgraphs
  Subgraphs->next->FillBondStructureFromReference(out, this, (FragmentCounter = 0), ListOfLocalAtoms, false);  // we want to keep the created ListOfLocalAtoms
  // analysis of the cycles (print rings, get minimum cycle length) for each subgraph
  for(int i=AtomCount;i--;)
    MinimumRingSize[i] = AtomCount;
  MolecularWalker = Subgraphs;
  FragmentCounter = 0;
  while (MolecularWalker->next != NULL) {
    MolecularWalker = MolecularWalker->next;
    *out << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
    LocalBackEdgeStack = new StackClass<bond *> (MolecularWalker->Leaf->BondCount);
//    // check the list of local atoms for debugging
//    *out << Verbose(0) << "ListOfLocalAtoms for this subgraph is:" << endl;
//    for (int i=0;i<AtomCount;i++)
//      if (ListOfLocalAtoms[FragmentCounter][i] == NULL)
//        *out << "\tNULL";
//      else
//        *out << "\t" << ListOfLocalAtoms[FragmentCounter][i]->Name;
    *out << Verbose(0) << "Gathering local back edges for subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
    MolecularWalker->Leaf->PickLocalBackEdges(out, ListOfLocalAtoms[FragmentCounter++], BackEdgeStack, LocalBackEdgeStack);
    *out << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
    MolecularWalker->Leaf->CyclicStructureAnalysis(out, LocalBackEdgeStack, MinimumRingSize);
    *out << Verbose(0) << "Done with Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
    delete(LocalBackEdgeStack);
  }

  // ===== 3. if structure still valid, parse key set file and others =====
  FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList);

  // ===== 4. check globally whether there's something to do actually (first adaptivity check)
  FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(out, configuration->configpath);

  // =================================== Begin of FRAGMENTATION ===============================
  // ===== 6a. assign each keyset to its respective subgraph =====
  Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), true);

  // ===== 6b. prepare and go into the adaptive (Order<0), single-step (Order==0) or incremental (Order>0) cycle
  KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
  AtomMask = new bool[AtomCount+1];
  AtomMask[AtomCount] = false;
  FragmentationToDo = false;  // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards
  while ((CheckOrder = CheckOrderAtSite(out, AtomMask, ParsedFragmentList, Order, MinimumRingSize, configuration->configpath))) {
    FragmentationToDo = FragmentationToDo || CheckOrder;
    AtomMask[AtomCount] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
    // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
    Subgraphs->next->FillRootStackForSubgraphs(out, RootStack, AtomMask, (FragmentCounter = 0));

    // ===== 7. fill the bond fragment list =====
    FragmentCounter = 0;
    MolecularWalker = Subgraphs;
    while (MolecularWalker->next != NULL) {
      MolecularWalker = MolecularWalker->next;
      *out << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
      //MolecularWalker->Leaf->OutputListOfBonds(out);  // output ListOfBondsPerAtom for debugging
      if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
        // call BOSSANOVA method
        *out << Verbose(0) << endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= " << endl;
        MolecularWalker->Leaf->FragmentBOSSANOVA(out, FragmentList[FragmentCounter], RootStack[FragmentCounter], MinimumRingSize);
      } else {
        cerr << "Subgraph " << MolecularWalker << " has no atoms!" << endl;
      }
      FragmentCounter++;  // next fragment list
    }
  }
  delete[](RootStack);
  delete[](AtomMask);
  delete(ParsedFragmentList);
  delete[](MinimumRingSize);


  // ==================================== End of FRAGMENTATION ============================================

  // ===== 8a. translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
  Subgraphs->next->TranslateIndicesToGlobalIDs(out, FragmentList, (FragmentCounter = 0), TotalNumberOfKeySets, TotalGraph);

  // free subgraph memory again
  FragmentCounter = 0;
  if (Subgraphs != NULL) {
    while (Subgraphs->next != NULL) {
      Subgraphs = Subgraphs->next;
      delete(FragmentList[FragmentCounter++]);
      delete(Subgraphs->previous);
    }
    delete(Subgraphs);
  }
  Free(&FragmentList);

  // ===== 8b. gather keyset lists (graphs) from all subgraphs and transform into MoleculeListClass =====
  //if (FragmentationToDo) {    // we should always store the fragments again as coordination might have changed slightly without changing bond structure
    // allocate memory for the pointer array and transmorph graphs into full molecular fragments
    BondFragments = new MoleculeListClass();
    int k=0;
    for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
      KeySet test = (*runner).first;
      *out << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl;
      BondFragments->insert(StoreFragmentFromKeySet(out, test, configuration));
      k++;
    }
    *out << k << "/" << BondFragments->ListOfMolecules.size() << " fragments generated from the keysets." << endl;

    // ===== 9. Save fragments' configuration and keyset files et al to disk ===
    if (BondFragments->ListOfMolecules.size() != 0) {
      // create the SortIndex from BFS labels to order in the config file
      CreateMappingLabelsToConfigSequence(out, SortIndex);

      *out << Verbose(1) << "Writing " << BondFragments->ListOfMolecules.size() << " possible bond fragmentation configs" << endl;
      if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
        *out << Verbose(1) << "All configs written." << endl;
      else
        *out << Verbose(1) << "Some config writing failed." << endl;

      // store force index reference file
      BondFragments->StoreForcesFile(out, configuration->configpath, SortIndex);

      // store keysets file
      StoreKeySetFile(out, TotalGraph, configuration->configpath);

      // store Adjacency file
      StoreAdjacencyToFile(out, configuration->configpath);

      // store Hydrogen saturation correction file
      BondFragments->AddHydrogenCorrection(out, configuration->configpath);

      // store adaptive orders into file
      StoreOrderAtSiteFile(out, configuration->configpath);

      // restore orbital and Stop values
      CalculateOrbitals(*configuration);

      // free memory for bond part
      *out << Verbose(1) << "Freeing bond memory" << endl;
      delete(FragmentList); // remove bond molecule from memory
      Free(&SortIndex);
    } else
      *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
  //} else
  //  *out << Verbose(1) << "No fragments to store." << endl;
  *out << Verbose(0) << "End of bond fragmentation." << endl;

  return ((int)(!FragmentationToDo)+1);    // 1 - continue, 2 - stop (no fragmentation occured)
};


/** Stores pairs (Atom::nr, Atom::AdaptiveOrder) into file.
 * Atoms not present in the file get "-1".
 * \param *out output stream for debugging
 * \param *path path to file ORDERATSITEFILE
 * \return true - file writable, false - not writable
 */
bool molecule::StoreOrderAtSiteFile(ofstream *out, char *path)
{
  stringstream line;
  ofstream file;

  line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
  file.open(line.str().c_str());
  *out << Verbose(1) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl;
  if (file != NULL) {
    ActOnAllAtoms( &atom::OutputOrder, &file );
    file.close();
    *out << Verbose(1) << "done." << endl;
    return true;
  } else {
    *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
    return false;
  }
};

/** Parses pairs(Atom::nr, Atom::AdaptiveOrder) from file and stores in molecule's Atom's.
 * Atoms not present in the file get "0".
 * \param *out output stream for debugging
 * \param *path path to file ORDERATSITEFILEe
 * \return true - file found and scanned, false - file not found
 * \sa ParseKeySetFile() and CheckAdjacencyFileAgainstMolecule() as this is meant to be used in conjunction with the two
 */
bool molecule::ParseOrderAtSiteFromFile(ofstream *out, char *path)
{
  unsigned char *OrderArray = Malloc<unsigned char>(AtomCount, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
  bool *MaxArray = Malloc<bool>(AtomCount, "molecule::ParseOrderAtSiteFromFile - *MaxArray");
  bool status;
  int AtomNr, value;
  stringstream line;
  ifstream file;

  *out << Verbose(1) << "Begin of ParseOrderAtSiteFromFile" << endl;
  for(int i=AtomCount;i--;)
    OrderArray[i] = 0;
  line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
  file.open(line.str().c_str());
  if (file != NULL) {
    for (int i=AtomCount;i--;) { // initialise with 0
      OrderArray[i] = 0;
      MaxArray[i] = 0;
    }
    while (!file.eof()) { // parse from file
      AtomNr = -1;
      file >> AtomNr;
      if (AtomNr != -1) {   // test whether we really parsed something (this is necessary, otherwise last atom is set twice and to 0 on second time)
        file >> value;
        OrderArray[AtomNr] = value;
        file >> value;
        MaxArray[AtomNr] = value;
        //*out << Verbose(2) << "AtomNr " << AtomNr << " with order " << (int)OrderArray[AtomNr] << " and max order set to " << (int)MaxArray[AtomNr] << "." << endl;
      }
    }
    file.close();

    // set atom values
    SetAtomValueToIndexedArray( OrderArray, &atom::nr, &atom::AdaptiveOrder );
    SetAtomValueToIndexedArray( MaxArray, &atom::nr, &atom::MaxOrder );

    *out << Verbose(1) << "done." << endl;
    status = true;
  } else {
    *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
    status = false;
  }
  Free(&OrderArray);
  Free(&MaxArray);

  *out << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl;
  return status;
};



/** Looks through a StackClass<atom *> and returns the likeliest removal candiate.
 * \param *out output stream for debugging messages
 * \param *&Leaf KeySet to look through
 * \param *&ShortestPathList list of the shortest path to decide which atom to suggest as removal candidate in the end
 * \param index of the atom suggested for removal
 */
int molecule::LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList)
{
  atom *Runner = NULL;
  int SP, Removal;

  *out << Verbose(2) << "Looking for removal candidate." << endl;
  SP = -1; //0;  // not -1, so that Root is never removed
  Removal = -1;
  for (KeySet::iterator runner = Leaf->begin(); runner != Leaf->end(); runner++) {
    Runner = FindAtom((*runner));
    if (Runner->type->Z != 1) { // skip all those added hydrogens when re-filling snake stack
      if (ShortestPathList[(*runner)] > SP) {  // remove the oldest one with longest shortest path
        SP = ShortestPathList[(*runner)];
        Removal = (*runner);
      }
    }
  }
  return Removal;
};

/** Initializes some value for putting fragment of \a *mol into \a *Leaf.
 * \param *mol total molecule
 * \param *Leaf fragment molecule
 * \param &Leaflet pointer to KeySet structure
 * \param **SonList list which atom of \a *Leaf is a son of which atom in \a *mol
 * \return number of atoms in fragment
 */
int StoreFragmentFromKeySet_Init(molecule *mol, molecule *Leaf, KeySet &Leaflet, atom **SonList)
{
  atom *FatherOfRunner = NULL;

  Leaf->BondDistance = mol->BondDistance;
  for(int i=NDIM*2;i--;)
    Leaf->cell_size[i] = mol->cell_size[i];

  // initialise SonList (indicates when we need to replace a bond with hydrogen instead)
  for(int i=mol->AtomCount;i--;)
    SonList[i] = NULL;

  // first create the minimal set of atoms from the KeySet
  int size = 0;
  for(KeySet::iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
    FatherOfRunner = mol->FindAtom((*runner));  // find the id
    SonList[FatherOfRunner->nr] = Leaf->AddCopyAtom(FatherOfRunner);
    size++;
  }
  return size;
};

/** Creates an induced subgraph out of a fragmental key set, adding bonds and hydrogens (if treated specially).
 * \param *out output stream for debugging messages
 * \param *mol total molecule
 * \param *Leaf fragment molecule
 * \param IsAngstroem whether we have Ansgtroem or bohrradius
 * \param **SonList list which atom of \a *Leaf is a son of which atom in \a *mol
 */
void CreateInducedSubgraphOfFragment(ofstream *out, molecule *mol, molecule *Leaf, atom **SonList, bool IsAngstroem)
{
  bool LonelyFlag = false;
  atom *OtherFather = NULL;
  atom *FatherOfRunner = NULL;
  Leaf->CountAtoms(out);

  atom *Runner = Leaf->start;
  while (Runner->next != Leaf->end) {
    Runner = Runner->next;
    LonelyFlag = true;
    FatherOfRunner = Runner->father;
    if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
      // create all bonds
      for (int i=0;i<mol->NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father
        OtherFather = mol->ListOfBondsPerAtom[FatherOfRunner->nr][i]->GetOtherAtom(FatherOfRunner);
//        *out << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->nr] << " is bound to " << *OtherFather;
        if (SonList[OtherFather->nr] != NULL) {
//          *out << ", whose son is " << *SonList[OtherFather->nr] << "." << endl;
          if (OtherFather->nr > FatherOfRunner->nr) { // add bond (nr check is for adding only one of both variants: ab, ba)
//            *out << Verbose(3) << "Adding Bond: ";
//            *out <<
            Leaf->AddBond(Runner, SonList[OtherFather->nr], mol->ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree);
//            *out << "." << endl;
            //NumBonds[Runner->nr]++;
          } else {
//            *out << Verbose(3) << "Not adding bond, labels in wrong order." << endl;
          }
          LonelyFlag = false;
        } else {
//          *out << ", who has no son in this fragment molecule." << endl;
#ifdef ADDHYDROGEN
          //*out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
          if(!Leaf->AddHydrogenReplacementAtom(out, mol->ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, mol->ListOfBondsPerAtom[FatherOfRunner->nr],mol->NumberOfBondsPerAtom[FatherOfRunner->nr], IsAngstroem))
            exit(1);
#endif
          //NumBonds[Runner->nr] += ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree;
        }
      }
    } else {
      *out << Verbose(0) << "ERROR: Son " << Runner->Name << " has father " << FatherOfRunner->Name << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl;
    }
    if ((LonelyFlag) && (Leaf->AtomCount > 1)) {
      *out << Verbose(0) << *Runner << "has got bonds only to hydrogens!" << endl;
    }
#ifdef ADDHYDROGEN
    while ((Runner->next != Leaf->end) && (Runner->next->type->Z == 1)) // skip added hydrogen
      Runner = Runner->next;
#endif
  }
};

/** Stores a fragment from \a KeySet into \a molecule.
 * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
 * molecule and adds missing hydrogen where bonds were cut.
 * \param *out output stream for debugging messages
 * \param &Leaflet pointer to KeySet structure
 * \param IsAngstroem whether we have Ansgtroem or bohrradius
 * \return pointer to constructed molecule
 */
molecule * molecule::StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem)
{
  atom **SonList = Malloc<atom*>(AtomCount, "molecule::StoreFragmentFromStack: **SonList");
  molecule *Leaf = new molecule(elemente);

//  *out << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
  StoreFragmentFromKeySet_Init(this, Leaf, Leaflet, SonList);
  // create the bonds between all: Make it an induced subgraph and add hydrogen
//  *out << Verbose(2) << "Creating bonds from father graph (i.e. induced subgraph creation)." << endl;
  CreateInducedSubgraphOfFragment(out, this, Leaf, SonList, IsAngstroem);

  Leaf->CreateListOfBondsPerAtom(out);
  //Leaflet->Leaf->ScanForPeriodicCorrection(out);
  Free(&SonList);
//  *out << Verbose(1) << "End of StoreFragmentFromKeyset." << endl;
  return Leaf;
};

/** From a given set of Bond sorted by Shortest Path distance, create all possible fragments of size \a SetDimension.
 * -# loops over every possible combination (2^dimension of edge set)
 *  -# inserts current set, if there's still space left
 *    -# yes: calls SPFragmentGenerator with structure, created new edge list and size respective to root dist
ance+1
 *    -# no: stores fragment into keyset list by calling InsertFragmentIntoGraph
 *  -# removes all items added into the snake stack (in UniqueFragments structure) added during level (root
distance) and current set
 * \param *out output stream for debugging
 * \param FragmentSearch UniqueFragments structure with all values needed
 * \param RootDistance current shortest path level, whose set of edges is represented by **BondsSet
 * \param SetDimension Number of possible bonds on this level (i.e. size of the array BondsSet[])
 * \param SubOrder remaining number of allowed vertices to add
 */
void molecule::SPFragmentGenerator(ofstream *out, struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder)
{
  atom *OtherWalker = NULL;
  int verbosity = 0; //FragmentSearch->ANOVAOrder-SubOrder;
  int NumCombinations;
  bool bit;
  int bits, TouchedIndex, SubSetDimension, SP, Added;
  int Removal;
  int SpaceLeft;
  int *TouchedList = Malloc<int>(SubOrder + 1, "molecule::SPFragmentGenerator: *TouchedList");
  bond *Binder = NULL;
  bond **BondsList = NULL;
  KeySetTestPair TestKeySetInsert;

  NumCombinations = 1 << SetDimension;

  // Hier muessen von 1 bis NumberOfBondsPerAtom[Walker->nr] alle Kombinationen
  // von Endstuecken (aus den Bonds) hinzugefuegt werden und fuer verbleibende ANOVAOrder
  // rekursiv GraphCrawler in der naechsten Ebene aufgerufen werden

  *out << Verbose(1+verbosity) << "Begin of SPFragmentGenerator." << endl;
  *out << Verbose(1+verbosity) << "We are " << RootDistance << " away from Root, which is " << *FragmentSearch->Root << ", SubOrder is " << SubOrder << ", SetDimension is " << SetDimension << " and this means " <<  NumCombinations-1 << " combination(s)." << endl;

  // initialised touched list (stores added atoms on this level)
  *out << Verbose(1+verbosity) << "Clearing touched list." << endl;
  for (TouchedIndex=SubOrder+1;TouchedIndex--;)  // empty touched list
    TouchedList[TouchedIndex] = -1;
  TouchedIndex = 0;

  // create every possible combination of the endpieces
  *out << Verbose(1+verbosity) << "Going through all combinations of the power set." << endl;
  for (int i=1;i<NumCombinations;i++) {  // sweep through all power set combinations (skip empty set!)
    // count the set bit of i
    bits = 0;
    for (int j=SetDimension;j--;)
      bits += (i & (1 << j)) >> j;

    *out << Verbose(1+verbosity) << "Current set is " << Binary(i | (1 << SetDimension)) << ", number of bits is " << bits << "." << endl;
    if (bits <= SubOrder) { // if not greater than additional atoms allowed on stack, continue
      // --1-- add this set of the power set of bond partners to the snake stack
      Added = 0;
      for (int j=0;j<SetDimension;j++) {  // pull out every bit by shifting
        bit = ((i & (1 << j)) != 0);  // mask the bit for the j-th bond
        if (bit) {  // if bit is set, we add this bond partner
          OtherWalker = BondsSet[j]->rightatom;  // rightatom is always the one more distant, i.e. the one to add
          //*out << Verbose(1+verbosity) << "Current Bond is " << ListOfBondsPerAtom[Walker->nr][i] << ", checking on " << *OtherWalker << "." << endl;
          *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
          TestKeySetInsert = FragmentSearch->FragmentSet->insert(OtherWalker->nr);
          if (TestKeySetInsert.second) {
            TouchedList[TouchedIndex++] = OtherWalker->nr;  // note as added
            Added++;
          } else {
            *out << Verbose(2+verbosity) << "This was item was already present in the keyset." << endl;
          }
            //FragmentSearch->UsedList[OtherWalker->nr][i] = true;
          //}
        } else {
          *out << Verbose(2+verbosity) << "Not adding." << endl;
        }
      }

      SpaceLeft = SubOrder - Added ;// SubOrder - bits; // due to item's maybe being already present, this does not work anymore
      if (SpaceLeft > 0) {
        *out << Verbose(1+verbosity) << "There's still some space left on stack: " << SpaceLeft << "." << endl;
        if (SubOrder > 1) {    // Due to Added above we have to check extra whether we're not already reaching beyond the desired Order
          // --2-- look at all added end pieces of this combination, construct bond subsets and sweep through a power set of these by recursion
          SP = RootDistance+1;  // this is the next level
          // first count the members in the subset
          SubSetDimension = 0;
          Binder = FragmentSearch->BondsPerSPList[2*SP];    // start node for this level
          while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {    // compare to end node of this level
            Binder = Binder->next;
            for (int k=TouchedIndex;k--;) {
              if (Binder->Contains(TouchedList[k]))   // if we added this very endpiece
                SubSetDimension++;
            }
          }
          // then allocate and fill the list
          BondsList = Malloc<bond*>(SubSetDimension, "molecule::SPFragmentGenerator: **BondsList");
          SubSetDimension = 0;
          Binder = FragmentSearch->BondsPerSPList[2*SP];
          while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {
            Binder = Binder->next;
            for (int k=0;k<TouchedIndex;k++) {
              if (Binder->leftatom->nr == TouchedList[k])   // leftatom is always the close one
                BondsList[SubSetDimension++] = Binder;
            }
          }
          *out << Verbose(2+verbosity) << "Calling subset generator " << SP << " away from root " << *FragmentSearch->Root << " with sub set dimension " << SubSetDimension << "." << endl;
          SPFragmentGenerator(out, FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits);
          Free(&BondsList);
        }
      } else {
        // --2-- otherwise store the complete fragment
        *out << Verbose(1+verbosity) << "Enough items on stack for a fragment!" << endl;
        // store fragment as a KeySet
        *out << Verbose(2) << "Found a new fragment[" << FragmentSearch->FragmentCounter << "], local nr.s are: ";
        for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
          *out << (*runner) << " ";
        *out << endl;
        //if (!CheckForConnectedSubgraph(out, FragmentSearch->FragmentSet))
          //*out << Verbose(0) << "ERROR: The found fragment is not a connected subgraph!" << endl;
        InsertFragmentIntoGraph(out, FragmentSearch);
        //Removal = LookForRemovalCandidate(out, FragmentSearch->FragmentSet, FragmentSearch->ShortestPathList);
        //Removal = StoreFragmentFromStack(out, FragmentSearch->Root, FragmentSearch->Leaflet, FragmentSearch->FragmentStack, FragmentSearch->ShortestPathList, &FragmentSearch->FragmentCounter, FragmentSearch->configuration);
      }

      // --3-- remove all added items in this level from snake stack
      *out << Verbose(1+verbosity) << "Removing all items that were added on this SP level " << RootDistance << "." << endl;
      for(int j=0;j<TouchedIndex;j++) {
        Removal = TouchedList[j];
        *out << Verbose(2+verbosity) << "Removing item nr. " << Removal << " from snake stack." << endl;
        FragmentSearch->FragmentSet->erase(Removal);
        TouchedList[j] = -1;
      }
      *out << Verbose(2) << "Remaining local nr.s on snake stack are: ";
      for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
        *out << (*runner) << " ";
      *out << endl;
      TouchedIndex = 0; // set Index to 0 for list of atoms added on this level
    } else {
      *out << Verbose(2+verbosity) << "More atoms to add for this set (" << bits << ") than space left on stack " << SubOrder << ", skipping this set." << endl;
    }
  }
  Free(&TouchedList);
  *out << Verbose(1+verbosity) << "End of SPFragmentGenerator, " << RootDistance << " away from Root " << *FragmentSearch->Root << " and SubOrder is " << SubOrder << "." << endl;
};



/** Creates a list of all unique fragments of certain vertex size from a given graph \a Fragment for a given root vertex in the context of \a this molecule.
 * -# initialises UniqueFragments structure
 * -# fills edge list via BFS
 * -# creates the fragment by calling recursive function SPFragmentGenerator with UniqueFragments structure, 0 as
 root distance, the edge set, its dimension and the current suborder
 * -# Free'ing structure
 * Note that we may use the fact that the atoms are SP-ordered on the atomstack. I.e. when popping always the last, we first get all
 * with SP of 2, then those with SP of 3, then those with SP of 4 and so on.
 * \param *out output stream for debugging
 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
 * \param FragmentSearch UniqueFragments structure containing TEFactor, root atom and so on
 * \param RestrictedKeySet Restricted vertex set to use in context of molecule
 * \return number of inserted fragments
 * \note ShortestPathList in FragmentSearch structure is probably due to NumberOfAtomsSPLevel and SP not needed anymore
 */
int molecule::PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet)
{
  int SP, AtomKeyNr;
  atom *Walker = NULL, *OtherWalker = NULL, *Predecessor = NULL;
  bond *Binder = NULL;
  bond *CurrentEdge = NULL;
  bond **BondsList = NULL;
  int RootKeyNr = FragmentSearch.Root->GetTrueFather()->nr;
  int Counter = FragmentSearch.FragmentCounter;
  int RemainingWalkers;

  *out << endl;
  *out << Verbose(0) << "Begin of PowerSetGenerator with order " << Order << " at Root " << *FragmentSearch.Root << "." << endl;

  // prepare Label and SP arrays of the BFS search
  FragmentSearch.ShortestPathList[FragmentSearch.Root->nr] = 0;

  // prepare root level (SP = 0) and a loop bond denoting Root
  for (int i=1;i<Order;i++)
    FragmentSearch.BondsPerSPCount[i] = 0;
  FragmentSearch.BondsPerSPCount[0] = 1;
  Binder = new bond(FragmentSearch.Root, FragmentSearch.Root);
  add(Binder, FragmentSearch.BondsPerSPList[1]);

  // do a BFS search to fill the SP lists and label the found vertices
  // Actually, we should construct a spanning tree vom the root atom and select all edges therefrom and put them into
  // according shortest path lists. However, we don't. Rather we fill these lists right away, as they do form a spanning
  // tree already sorted into various SP levels. That's why we just do loops over the depth (CurrentSP) and breadth
  // (EdgeinSPLevel) of this tree ...
  // In another picture, the bonds always contain a direction by rightatom being the one more distant from root and hence
  // naturally leftatom forming its predecessor, preventing the BFS"seeker" from continuing in the wrong direction.
  *out << endl;
  *out << Verbose(0) << "Starting BFS analysis ..." << endl;
  for (SP = 0; SP < (Order-1); SP++) {
    *out << Verbose(1) << "New SP level reached: " << SP << ", creating new SP list with " << FragmentSearch.BondsPerSPCount[SP] << " item(s)";
    if (SP > 0) {
      *out << ", old level closed with " << FragmentSearch.BondsPerSPCount[SP-1] << " item(s)." << endl;
      FragmentSearch.BondsPerSPCount[SP] = 0;
    } else
      *out << "." << endl;

    RemainingWalkers = FragmentSearch.BondsPerSPCount[SP];
    CurrentEdge = FragmentSearch.BondsPerSPList[2*SP];    /// start of this SP level's list
    while (CurrentEdge->next != FragmentSearch.BondsPerSPList[2*SP+1]) {    /// end of this SP level's list
      CurrentEdge = CurrentEdge->next;
      RemainingWalkers--;
      Walker = CurrentEdge->rightatom;    // rightatom is always the one more distant
      Predecessor = CurrentEdge->leftatom;    // ... and leftatom is predecessor
      AtomKeyNr = Walker->nr;
      *out << Verbose(0) << "Current Walker is: " << *Walker << " with nr " << Walker->nr << " and SP of " << SP << ", with " << RemainingWalkers << " remaining walkers on this level." << endl;
      // check for new sp level
      // go through all its bonds
      *out << Verbose(1) << "Going through all bonds of Walker." << endl;
      for (int i=0;i<NumberOfBondsPerAtom[AtomKeyNr];i++) {
        Binder = ListOfBondsPerAtom[AtomKeyNr][i];
        OtherWalker = Binder->GetOtherAtom(Walker);
        if ((RestrictedKeySet.find(OtherWalker->nr) != RestrictedKeySet.end())
  #ifdef ADDHYDROGEN
         && (OtherWalker->type->Z != 1)
  #endif
                                                              ) {  // skip hydrogens and restrict to fragment
          *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *Binder << "." << endl;
          // set the label if not set (and push on root stack as well)
          if ((OtherWalker != Predecessor) && (OtherWalker->GetTrueFather()->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
            FragmentSearch.ShortestPathList[OtherWalker->nr] = SP+1;
            *out << Verbose(3) << "Set Shortest Path to " << FragmentSearch.ShortestPathList[OtherWalker->nr] << "." << endl;
            // add the bond in between to the SP list
            Binder = new bond(Walker, OtherWalker); // create a new bond in such a manner, that bond::rightatom is always the one more distant
            add(Binder, FragmentSearch.BondsPerSPList[2*(SP+1)+1]);
            FragmentSearch.BondsPerSPCount[SP+1]++;
            *out << Verbose(3) << "Added its bond to SP list, having now " << FragmentSearch.BondsPerSPCount[SP+1] << " item(s)." << endl;
          } else {
            if (OtherWalker != Predecessor)
              *out << Verbose(3) << "Not passing on, as index of " << *OtherWalker << " " << OtherWalker->GetTrueFather()->nr << " is smaller than that of Root " << RootKeyNr << "." << endl;
            else
              *out << Verbose(3) << "This is my predecessor " << *Predecessor << "." << endl;
          }
        } else *out << Verbose(2) << "Is not in the restricted keyset or skipping hydrogen " << *OtherWalker << "." << endl;
      }
    }
  }

  // outputting all list for debugging
  *out << Verbose(0) << "Printing all found lists." << endl;
  for(int i=1;i<Order;i++) {    // skip the root edge in the printing
    Binder = FragmentSearch.BondsPerSPList[2*i];
    *out << Verbose(1) << "Current SP level is " << i << "." << endl;
    while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
      Binder = Binder->next;
      *out << Verbose(2) << *Binder << endl;
    }
  }

  // creating fragments with the found edge sets  (may be done in reverse order, faster)
  SP = -1;  // the Root <-> Root edge must be subtracted!
  for(int i=Order;i--;) { // sum up all found edges
    Binder = FragmentSearch.BondsPerSPList[2*i];
    while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
      Binder = Binder->next;
      SP ++;
    }
  }
  *out << Verbose(0) << "Total number of edges is " << SP << "." << endl;
  if (SP >= (Order-1)) {
    // start with root (push on fragment stack)
    *out << Verbose(0) << "Starting fragment generation with " << *FragmentSearch.Root << ", local nr is " << FragmentSearch.Root->nr << "." << endl;
    FragmentSearch.FragmentSet->clear();
    *out << Verbose(0) << "Preparing subset for this root and calling generator." << endl;
    // prepare the subset and call the generator
    BondsList = Malloc<bond*>(FragmentSearch.BondsPerSPCount[0], "molecule::PowerSetGenerator: **BondsList");
    BondsList[0] = FragmentSearch.BondsPerSPList[0]->next;  // on SP level 0 there's only the root bond

    SPFragmentGenerator(out, &FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order);

    Free(&BondsList);
  } else {
    *out << Verbose(0) << "Not enough total number of edges to build " << Order << "-body fragments." << endl;
  }

  // as FragmentSearch structure is used only once, we don't have to clean it anymore
  // remove root from stack
  *out << Verbose(0) << "Removing root again from stack." << endl;
  FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);

  // free'ing the bonds lists
  *out << Verbose(0) << "Free'ing all found lists. and resetting index lists" << endl;
  for(int i=Order;i--;) {
    *out << Verbose(1) << "Current SP level is " << i << ": ";
    Binder = FragmentSearch.BondsPerSPList[2*i];
    while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
      Binder = Binder->next;
      // *out << "Removing atom " << Binder->leftatom->nr << " and " << Binder->rightatom->nr << "." << endl; // make sure numbers are local
      FragmentSearch.ShortestPathList[Binder->leftatom->nr] = -1;
      FragmentSearch.ShortestPathList[Binder->rightatom->nr] = -1;
    }
    // delete added bonds
    cleanup(FragmentSearch.BondsPerSPList[2*i], FragmentSearch.BondsPerSPList[2*i+1]);
    // also start and end node
    *out << "cleaned." << endl;
  }

  // return list
  *out << Verbose(0) << "End of PowerSetGenerator." << endl;
  return (FragmentSearch.FragmentCounter - Counter);
};

bool KeyCompare::operator() (const KeySet SubgraphA, const KeySet SubgraphB) const
{
  //cout << "my check is used." << endl;
  if (SubgraphA.size() < SubgraphB.size()) {
    return true;
  } else {
    if (SubgraphA.size() > SubgraphB.size()) {
      return false;
    } else {
      KeySet::iterator IteratorA = SubgraphA.begin();
      KeySet::iterator IteratorB = SubgraphB.begin();
      while ((IteratorA != SubgraphA.end()) && (IteratorB != SubgraphB.end())) {
        if ((*IteratorA) <  (*IteratorB))
          return true;
        else if ((*IteratorA) > (*IteratorB)) {
            return false;
          } // else, go on to next index
        IteratorA++;
        IteratorB++;
      } // end of while loop
    }// end of check in case of equal sizes
  }
  return false; // if we reach this point, they are equal
};


/** Performs BOSSANOVA decomposition at selected sites, increasing the cutoff by one at these sites.
 * -# constructs a complete keyset of the molecule
 * -# In a loop over all possible roots from the given rootstack
 *  -# increases order of root site
 *  -# calls PowerSetGenerator with this order, the complete keyset and the rootkeynr
 *  -# for all consecutive lower levels PowerSetGenerator is called with the suborder, the higher order keyset
as the restricted one and each site in the set as the root)
 *  -# these are merged into a fragment list of keysets
 * -# All fragment lists (for all orders, i.e. from all destination fields) are merged into one list for return
 * Important only is that we create all fragments, it is not important if we create them more than once
 * as these copies are filtered out via use of the hash table (KeySet).
 * \param *out output stream for debugging
 * \param Fragment&*List list of already present keystacks (adaptive scheme) or empty list
 * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
 * \param *MinimumRingSize minimum ring size for each atom (molecule::Atomcount)
 * \return pointer to Graph list
 */
void molecule::FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize)
{
  Graph ***FragmentLowerOrdersList = NULL;
  int NumLevels, NumMolecules, TotalNumMolecules = 0, *NumMoleculesOfOrder = NULL;
  int counter = 0, Order;
  int UpgradeCount = RootStack.size();
  KeyStack FragmentRootStack;
  int RootKeyNr, RootNr;
  struct UniqueFragments FragmentSearch;

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

  // FragmentLowerOrdersList is a 2D-array of pointer to MoleculeListClass objects, one dimension represents the ANOVA expansion of a single order (i.e. 5)
  // with all needed lower orders that are subtracted, the other dimension is the BondOrder (i.e. from 1 to 5)
  NumMoleculesOfOrder = Malloc<int>(UpgradeCount, "molecule::FragmentBOSSANOVA: *NumMoleculesOfOrder");
  FragmentLowerOrdersList = Malloc<Graph**>(UpgradeCount, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");

  // initialise the fragments structure
  FragmentSearch.ShortestPathList = Malloc<int>(AtomCount, "molecule::PowerSetGenerator: *ShortestPathList");
  FragmentSearch.FragmentCounter = 0;
  FragmentSearch.FragmentSet = new KeySet;
  FragmentSearch.Root = FindAtom(RootKeyNr);
  for (int i=AtomCount;i--;) {
    FragmentSearch.ShortestPathList[i] = -1;
  }

  // Construct the complete KeySet which we need for topmost level only (but for all Roots)
  atom *Walker = start;
  KeySet CompleteMolecule;
  while (Walker->next != end) {
    Walker = Walker->next;
    CompleteMolecule.insert(Walker->GetTrueFather()->nr);
  }

  // this can easily be seen: if Order is 5, then the number of levels for each lower order is the total sum of the number of levels above, as
  // each has to be split up. E.g. for the second level we have one from 5th, one from 4th, two from 3th (which in turn is one from 5th, one from 4th),
  // hence we have overall four 2th order levels for splitting. This also allows for putting all into a single array (FragmentLowerOrdersList[])
  // with the order along the cells as this: 5433222211111111 for BondOrder 5 needing 16=pow(2,5-1) cells (only we use bit-shifting which is faster)
  RootNr = 0;   // counts through the roots in RootStack
  while ((RootNr < UpgradeCount) && (!RootStack.empty())) {
    RootKeyNr = RootStack.front();
    RootStack.pop_front();
    Walker = FindAtom(RootKeyNr);
    // check cyclic lengths
    //if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) {
    //  *out << Verbose(0) << "Bond order " << Walker->GetTrueFather()->AdaptiveOrder << " of Root " << *Walker << " greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
    //} else
    {
      // increase adaptive order by one
      Walker->GetTrueFather()->AdaptiveOrder++;
      Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;

      // initialise Order-dependent entries of UniqueFragments structure
      FragmentSearch.BondsPerSPList = Malloc<bond*>(Order * 2, "molecule::PowerSetGenerator: ***BondsPerSPList");
      FragmentSearch.BondsPerSPCount = Malloc<int>(Order, "molecule::PowerSetGenerator: *BondsPerSPCount");
      for (int i=Order;i--;) {
        FragmentSearch.BondsPerSPList[2*i] = new bond();    // start node
        FragmentSearch.BondsPerSPList[2*i+1] = new bond();  // end node
        FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1];     // intertwine these two
        FragmentSearch.BondsPerSPList[2*i+1]->previous = FragmentSearch.BondsPerSPList[2*i];
        FragmentSearch.BondsPerSPCount[i] = 0;
      }

      // allocate memory for all lower level orders in this 1D-array of ptrs
      NumLevels = 1 << (Order-1); // (int)pow(2,Order);
      FragmentLowerOrdersList[RootNr] = Malloc<Graph*>(NumLevels, "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
      for (int i=0;i<NumLevels;i++)
        FragmentLowerOrdersList[RootNr][i] = NULL;

      // create top order where nothing is reduced
      *out << Verbose(0) << "==============================================================================================================" << endl;
      *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", " << (RootStack.size()-RootNr) << " Roots remaining." << endl; // , NumLevels is " << NumLevels << "

      // Create list of Graphs of current Bond Order (i.e. F_{ij})
      FragmentLowerOrdersList[RootNr][0] =  new Graph;
      FragmentSearch.TEFactor = 1.;
      FragmentSearch.Leaflet = FragmentLowerOrdersList[RootNr][0];      // set to insertion graph
      FragmentSearch.Root = Walker;
      NumMoleculesOfOrder[RootNr] = PowerSetGenerator(out, Walker->AdaptiveOrder, FragmentSearch, CompleteMolecule);
      *out << Verbose(1) << "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
      if (NumMoleculesOfOrder[RootNr] != 0) {
        NumMolecules = 0;

        // we don't have to dive into suborders! These keysets are all already created on lower orders!
        // this was all ancient stuff, when we still depended on the TEFactors (and for those the suborders were needed)

//        if ((NumLevels >> 1) > 0) {
//          // create lower order fragments
//          *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
//          Order = Walker->AdaptiveOrder;
//          for (int source=0;source<(NumLevels >> 1);source++) { // 1-terms don't need any more splitting, that's why only half is gone through (shift again)
//            // step down to next order at (virtual) boundary of powers of 2 in array
//            while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))
//              Order--;
//            *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
//            for (int SubOrder=Order-1;SubOrder>0;SubOrder--) {
//              int dest = source + (1 << (Walker->AdaptiveOrder-(SubOrder+1)));
//              *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
//              *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
//
//              // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
//              //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
//              //NumMolecules = 0;
//              FragmentLowerOrdersList[RootNr][dest] = new Graph;
//              for(Graph::iterator runner = (*FragmentLowerOrdersList[RootNr][source]).begin();runner != (*FragmentLowerOrdersList[RootNr][source]).end(); runner++) {
//                for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
//                  Graph TempFragmentList;
//                  FragmentSearch.TEFactor = -(*runner).second.second;
//                  FragmentSearch.Leaflet = &TempFragmentList;      // set to insertion graph
//                  FragmentSearch.Root = FindAtom(*sprinter);
//                  NumMoleculesOfOrder[RootNr] += PowerSetGenerator(out, SubOrder, FragmentSearch, (*runner).first);
//                  // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
//                  *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
//                  InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], TempFragmentList, &NumMolecules);
//                }
//              }
//              *out << Verbose(1) << "Number of resulting molecules for SubOrder " << SubOrder << " is: " << NumMolecules << "." << endl;
//            }
//          }
//        }
      } else {
        Walker->GetTrueFather()->MaxOrder = true;
//        *out << Verbose(1) << "Hence, we don't dive into SubOrders ... " << endl;
      }
      // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
      //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
      TotalNumMolecules += NumMoleculesOfOrder[RootNr];
//      *out << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
      RootStack.push_back(RootKeyNr); // put back on stack
      RootNr++;

      // free Order-dependent entries of UniqueFragments structure for next loop cycle
      Free(&FragmentSearch.BondsPerSPCount);
      for (int i=Order;i--;) {
        delete(FragmentSearch.BondsPerSPList[2*i]);
        delete(FragmentSearch.BondsPerSPList[2*i+1]);
      }
      Free(&FragmentSearch.BondsPerSPList);
    }
  }
  *out << Verbose(0) << "==============================================================================================================" << endl;
  *out << Verbose(1) << "Total number of resulting molecules is: " << TotalNumMolecules << "." << endl;
  *out << Verbose(0) << "==============================================================================================================" << endl;

  // cleanup FragmentSearch structure
  Free(&FragmentSearch.ShortestPathList);
  delete(FragmentSearch.FragmentSet);

  // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein)
  // 5433222211111111
  // 43221111
  // 3211
  // 21
  // 1

  // Subsequently, we combine all into a single list (FragmentList)

  *out << Verbose(0) << "Combining the lists of all orders per order and finally into a single one." << endl;
  if (FragmentList == NULL) {
    FragmentList = new Graph;
    counter = 0;
  } else {
    counter = FragmentList->size();
  }
  RootNr = 0;
  while (!RootStack.empty()) {
    RootKeyNr = RootStack.front();
    RootStack.pop_front();
    Walker = FindAtom(RootKeyNr);
    NumLevels = 1 << (Walker->AdaptiveOrder - 1);
    for(int i=0;i<NumLevels;i++) {
      if (FragmentLowerOrdersList[RootNr][i] != NULL) {
        InsertGraphIntoGraph(out, *FragmentList, (*FragmentLowerOrdersList[RootNr][i]), &counter);
        delete(FragmentLowerOrdersList[RootNr][i]);
      }
    }
    Free(&FragmentLowerOrdersList[RootNr]);
    RootNr++;
  }
  Free(&FragmentLowerOrdersList);
  Free(&NumMoleculesOfOrder);

  *out << Verbose(0) << "End of FragmentBOSSANOVA." << endl;
};

/** Corrects the nuclei position if the fragment was created over the cell borders.
 * Scans all bonds, checks the distance, if greater than typical, we have a candidate for the correction.
 * We remove the bond whereafter the graph probably separates. Then, we translate the one component periodically
 * and re-add the bond. Looping on the distance check.
 * \param *out ofstream for debugging messages
 */
void molecule::ScanForPeriodicCorrection(ofstream *out)
{
  bond *Binder = NULL;
  bond *OtherBinder = NULL;
  atom *Walker = NULL;
  atom *OtherWalker = NULL;
  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
  enum Shading *ColorList = NULL;
  double tmp;
  Vector Translationvector;
  //class StackClass<atom *> *CompStack = NULL;
  class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
  bool flag = true;

  *out << Verbose(2) << "Begin of ScanForPeriodicCorrection." << endl;

  ColorList = Malloc<enum Shading>(AtomCount, "molecule::ScanForPeriodicCorrection: *ColorList");
  while (flag) {
    // remove bonds that are beyond bonddistance
    for(int i=NDIM;i--;)
      Translationvector.x[i] = 0.;
    // scan all bonds
    Binder = first;
    flag = false;
    while ((!flag) && (Binder->next != last)) {
      Binder = Binder->next;
      for (int i=NDIM;i--;) {
        tmp = fabs(Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i]);
        //*out << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
        if (tmp > BondDistance) {
          OtherBinder = Binder->next; // note down binding partner for later re-insertion
          unlink(Binder);   // unlink bond
          *out << Verbose(2) << "Correcting at bond " << *Binder << "." << endl;
          flag = true;
          break;
        }
      }
    }
    if (flag) {
      // create translation vector from their periodically modified distance
      for (int i=NDIM;i--;) {
        tmp = Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i];
        if (fabs(tmp) > BondDistance)
          Translationvector.x[i] = (tmp < 0) ? +1. : -1.;
      }
      Translationvector.MatrixMultiplication(matrix);
      //*out << Verbose(3) << "Translation vector is ";
      Translationvector.Output(out);
      *out << endl;
      // apply to all atoms of first component via BFS
      for (int i=AtomCount;i--;)
        ColorList[i] = white;
      AtomStack->Push(Binder->leftatom);
      while (!AtomStack->IsEmpty()) {
        Walker = AtomStack->PopFirst();
        //*out << Verbose (3) << "Current Walker is: " << *Walker << "." << endl;
        ColorList[Walker->nr] = black;    // mark as explored
        Walker->x.AddVector(&Translationvector); // translate
        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {  // go through all binding partners
          if (ListOfBondsPerAtom[Walker->nr][i] != Binder) {
            OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
            if (ColorList[OtherWalker->nr] == white) {
              AtomStack->Push(OtherWalker); // push if yet unexplored
            }
          }
        }
      }
      // re-add bond
      link(Binder, OtherBinder);
    } else {
      *out << Verbose(3) << "No corrections for this fragment." << endl;
    }
    //delete(CompStack);
  }

  // free allocated space from ReturnFullMatrixforSymmetric()
  delete(AtomStack);
  Free(&ColorList);
  Free(&matrix);
  *out << Verbose(2) << "End of ScanForPeriodicCorrection." << endl;
};