source: src/molecule_fragmentation.cpp@ c38826

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

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

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "Helpers/MemDebug.hpp"

#include <cstring>

#include "World.hpp"
#include "atom.hpp"
#include "bond.hpp"
#include "config.hpp"
#include "element.hpp"
#include "Helpers/helpers.hpp"
#include "lists.hpp"
#include "Helpers/Verbose.hpp"
#include "Helpers/Log.hpp"
#include "molecule.hpp"
#include "periodentafel.hpp"
#include "World.hpp"
#include "LinearAlgebra/Matrix.hpp"
#include "Box.hpp"
#include "stackclass.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(int order)
{
  size_t c = 0;
  int FragmentCount;
  // get maximum bond degree
  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    c = ((*iter)->ListOfBonds.size() > c) ? (*iter)->ListOfBonds.size() : c;
  }
  FragmentCount = NoNonHydrogen*(1 << (c*order));
  DoLog(1) && (Log() << 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(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
  }
  DoLog(1) && (Log() << Verbose(1) << "The scanned KeySet is ");
  for(KeySet::iterator runner = CurrentSet.begin(); runner != CurrentSet.end(); runner++) {
    DoLog(0) && (Log() << Verbose(0) << (*runner) << "\t");
  }
  DoLog(0) && (Log() << Verbose(0) << 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 &path path to file
 * \param *FragmentList empty, filled on return
 * \return true - parsing successfully, false - failure on parsing (FragmentList will be NULL)
 */
bool ParseKeySetFile(std::string &path, Graph *&FragmentList)
{
  bool status = true;
  ifstream InputFile;
  stringstream line;
  GraphTestPair testGraphInsert;
  int NumberOfFragments = 0;
  string filename;

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

  // 1st pass: open file and read
  DoLog(1) && (Log() << Verbose(1) << "Parsing the KeySet file ... " << endl);
  filename = path + KEYSETFILE;
  InputFile.open(filename.c_str());
  if (InputFile.good()) {
    // each line represents a new fragment
    char buffer[MAXSTRINGSIZE];
    // 1. parse keysets and insert into temp. graph
    while (!InputFile.eof()) {
      InputFile.getline(buffer, MAXSTRINGSIZE);
      KeySet CurrentSet;
      if ((strlen(buffer) > 0) && (ScanBufferIntoKeySet(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) {
          DoeLog(0) && (eLog()<< Verbose(0) << "KeySet file must be corrupt as there are two equal key sets therein!" << endl);
          performCriticalExit();
        }
      }
    }
    // 2. Free and done
    InputFile.close();
    InputFile.clear();
    DoLog(1) && (Log() << Verbose(1) << "\t ... done." << endl);
  } else {
    DoLog(1) && (Log() << Verbose(1) << "\t ... 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(char *path, Graph *FragmentList)
{
  bool status = true;
  ifstream InputFile;
  stringstream line;
  GraphTestPair testGraphInsert;
  int NumberOfFragments = 0;
  double TEFactor;
  char filename[MAXSTRINGSIZE];

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

  // 2nd pass: open TEFactors file and read
  DoLog(1) && (Log() << 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;
        DoLog(2) && (Log() << Verbose(2) << "Setting " << ++NumberOfFragments << " fragment's TEFactor to " << (*runner).second.second << "." << endl);
      } else {
        status = false;
        break;
      }
    }
    // 4. Free and done
    InputFile.close();
    DoLog(1) && (Log() << Verbose(1) << "done." << endl);
  } else {
    DoLog(1) && (Log() << Verbose(1) << "File " << filename << " not found." << endl);
    status = false;
  }

  return status;
};

/** Stores key sets to file.
 * \param KeySetList Graph with Keysets
 * \param &path path to file
 * \return true - file written successfully, false - writing failed
 */
bool StoreKeySetFile(Graph &KeySetList, std::string &path)
{
  bool status =  true;
  string line = path + KEYSETFILE;
  ofstream output(line.c_str());

  // open KeySet file
  DoLog(1) && (Log() << Verbose(1) << "Saving key sets of the total graph ... ");
  if(output.good()) {
    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;
    }
    DoLog(0) && (Log() << Verbose(0) << "done." << endl);
  } else {
    DoeLog(0) && (eLog()<< Verbose(0) << "Unable to open " << line << " for writing keysets!" << endl);
    performCriticalExit();
    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(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);
  DoLog(1) && (Log() << 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;
    DoLog(1) && (Log() << Verbose(1) << "done." << endl);
  } else {
    DoLog(1) && (Log() << 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(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;
        DoLog(2) && (Log() << Verbose(2) << "Updated element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl);
      } else {
        DoLog(2) && (Log() << Verbose(2) << "Did not update element " <<  (*PresentItem).first << " as " << FragOrder << " is less than or equal to " << (*PresentItem).second.second << "." << endl);
      }
    } else {
      DoLog(2) && (Log() << Verbose(2) << "Inserted element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl);
    }
  } else {
    DoLog(1) && (Log() << Verbose(1) << "No Fragment under No. " << No << "found." << endl);
  }
};

/** Scans the adaptive order file and insert (index, value) into map.
 * \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(std::string &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[MAXSTRINGSIZE];
  string filename = path + ENERGYPERFRAGMENT;
  ifstream InputFile(filename.c_str());

  if (InputFile.fail()) {
    DoeLog(1) && (eLog() << Verbose(1) << "Cannot find file " << filename << "." << endl);
    return AdaptiveCriteriaList;
  }

  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) {
        //Log() << 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
        //Log() << Verbose(2) << " - yields (" << No << "," << Value << ", " << FragOrder << ")" << endl;

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

  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(map<int, pair<double,int> > *AdaptiveCriteriaList, molecule *mol)
{
  atom *Walker = NULL;
  map<double, pair<int,int> > *FinalRootCandidates = new map<double, pair<int,int> > ;
  DoLog(1) && (Log() << 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) {
        DoLog(2) && (Log() << 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 {
        DoLog(2) && (Log() << Verbose(2) << "Excluding (" << *Walker << ", " << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "]), as it has reached its maximum order." << endl);
      }
    } else {
      DoeLog(0) && (eLog()<< Verbose(0) << "Atom No. " << (*runner).second.first << " was not found in this molecule." << endl);
      performCriticalExit();
    }
  }
  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(bool *AtomMask, map<double, pair<int,int> > &FinalRootCandidates, int Order, molecule *mol)
{
  atom *Walker = NULL;
  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]) {
      DoLog(2) && (Log() << 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
      //Log() << 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(bool *AtomMask, int AtomCount)
{
  DoLog(2) && (Log() << Verbose(2) << "              ");
  for(int i=0;i<AtomCount;i++)
    DoLog(0) && (Log() << Verbose(0) << (i % 10));
  DoLog(0) && (Log() << Verbose(0) << endl);
  DoLog(2) && (Log() << Verbose(2) << "Atom mask is: ");
  for(int i=0;i<AtomCount;i++)
    DoLog(0) && (Log() << Verbose(0) << (AtomMask[i] ? "t" : "f"));
  DoLog(0) && (Log() << Verbose(0) << endl);
};

/** Checks whether the OrderAtSite is still below \a Order at some site.
 * \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(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, std::string path)
{
  bool status = false;

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

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

    // transmorph graph keyset list into indexed KeySetList
    if (GlobalKeySetList == NULL) {
      DoeLog(1) && (eLog()<< Verbose(1) << "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(path, *IndexKeySetList); // (Root No., (Value, Order)) !
    if (AdaptiveCriteriaList->empty()) {
      DoeLog(2) && (eLog()<< Verbose(2) << "Unable to parse file, incrementing all." << endl);
      for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    #ifdef ADDHYDROGEN
        if ((*iter)->getType()->getAtomicNumber() != 1) // skip hydrogen
    #endif
        {
          AtomMask[(*iter)->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(AdaptiveCriteriaList, this);

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

    delete[](IndexKeySetList);
    delete[](AdaptiveCriteriaList);
    delete[](FinalRootCandidates);
  } else { // global increase of Bond Order
    for(molecule::const_iterator iter = begin(); iter != end(); ++iter) {
  #ifdef ADDHYDROGEN
      if ((*iter)->getType()->getAtomicNumber() != 1) // skip hydrogen
  #endif
      {
        AtomMask[(*iter)->nr] = true;  // include all (non-hydrogen) atoms
        if ((Order != 0) && ((*iter)->AdaptiveOrder < Order)) // && ((*iter)->AdaptiveOrder < MinimumRingSize[(*iter)->nr]))
          status = true;
      }
    }
    if ((!Order) && (!AtomMask[getAtomCount()]))  // single stepping, just check
      status = true;

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

  PrintAtomMask(AtomMask, getAtomCount()); // 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(int *&SortIndex)
{
  if (SortIndex != NULL) {
    DoLog(1) && (Log() << Verbose(1) << "SortIndex is " << SortIndex << " and not NULL as expected." << endl);
    return false;
  }
  SortIndex = new int[getAtomCount()];
  for(int i=getAtomCount();i--;)
    SortIndex[i] = -1;

  int AtomNo = 0;
  for(internal_iterator iter=atoms.begin();iter!=atoms.end();++iter){
    ASSERT(SortIndex[(*iter)->nr]==-1,"Same SortIndex set twice");
    SortIndex[(*iter)->nr] = AtomNo++;
  }
  //SetIndexedArrayForEachAtomTo( SortIndex, &atom::nr, &IncrementalAbsoluteValue, AtomNo );

  return true;
};



/** Creates a lookup table for true father's Atom::Nr -> atom ptr.
 * \param *start begin of list (STL iterator, i.e. first item)
 * \paran *end end of list (STL iterator, i.e. one past last item)
 * \param **Lookuptable pointer to return allocated lookup table (should be NULL on start)
 * \param count optional predetermined size for table (otherwise we set the count to highest true father id)
 * \return true - success, false - failure
 */
bool molecule::CreateFatherLookupTable(atom **&LookupTable, int count)
{
  bool status = true;
  int AtomNo;

  if (LookupTable != NULL) {
    Log() << Verbose(0) << "Pointer for Lookup table is not NULL! Aborting ..." <<endl;
    return false;
  }

  // count them
  if (count == 0) {
    for (molecule::iterator iter = begin(); iter != end(); ++iter) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron
      count = (count < (*iter)->GetTrueFather()->nr) ? (*iter)->GetTrueFather()->nr : count;
    }
  }
  if (count <= 0) {
    Log() << Verbose(0) << "Count of lookup list is 0 or less." << endl;
    return false;
  }

  // allocate and fill
  LookupTable = new atom *[count];
  if (LookupTable == NULL) {
    eLog() << Verbose(0) << "LookupTable memory allocation failed!" << endl;
    performCriticalExit();
    status = false;
  } else {
    for (int i=0;i<count;i++)
      LookupTable[i] = NULL;
    for (molecule::iterator iter = begin(); iter != end(); ++iter) {
      AtomNo = (*iter)->GetTrueFather()->nr;
      if ((AtomNo >= 0) && (AtomNo < count)) {
        //*out << "Setting LookupTable[" << AtomNo << "] to " << *(*iter) << endl;
        LookupTable[AtomNo] = (*iter);
      } else {
        Log() << Verbose(0) << "Walker " << *(*iter) << " exceeded range of nuclear ids [0, " << count << ")." << endl;
        status = false;
        break;
      }
    }
  }

  return status;
};

/** 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 Order up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
 * \param &prefix path and prefix of the bond order configs to be written
 * \return 1 - continue, 2 - stop (no fragmentation occured)
 */
int molecule::FragmentMolecule(int Order, std::string &prefix)
{
  MoleculeListClass *BondFragments = NULL;
  int *MinimumRingSize = new int[getAtomCount()];
  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;

  DoLog(0) && (Log() << Verbose(0) << endl);
#ifdef ADDHYDROGEN
  DoLog(0) && (Log() << Verbose(0) << "I will treat hydrogen special and saturate dangling bonds with it." << endl);
#else
  DoLog(0) && (Log() << 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 ====

  // create lookup table for Atom::nr
  FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(ListOfAtoms, getAtomCount());

  // === compare it with adjacency file ===
  FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(prefix, ListOfAtoms);
  delete[](ListOfAtoms);

  // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
  Subgraphs = DepthFirstSearchAnalysis(BackEdgeStack);

  // analysis of the cycles (print rings, get minimum cycle length) for each subgraph
  for(int i=getAtomCount();i--;)
    MinimumRingSize[i] = getAtomCount();
  MolecularWalker = Subgraphs;
  const int LeafCount = Subgraphs->next->Count();
  FragmentCounter = 0;
  while (MolecularWalker->next != NULL) {
    MolecularWalker = MolecularWalker->next;
    // fill the bond structure of the individually stored subgraphs
    ListOfAtoms = NULL;
    MolecularWalker->FillBondStructureFromReference(this, ListOfAtoms, false);  // we want to keep the created ListOfLocalAtoms
    DoLog(0) && (Log() << 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
//    Log() << Verbose(0) << "ListOfLocalAtoms for this subgraph is:" << endl;
//    for (int i=0;i<getAtomCount();i++)
//      if (ListOfLocalAtoms[FragmentCounter][i] == NULL)
//        Log() << Verbose(0) << "\tNULL";
//      else
//        Log() << Verbose(0) << "\t" << ListOfLocalAtoms[FragmentCounter][i]->Name;
    DoLog(0) && (Log() << Verbose(0) << "Gathering local back edges for subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl);
    MolecularWalker->Leaf->PickLocalBackEdges(ListOfAtoms, BackEdgeStack, LocalBackEdgeStack);
    DoLog(0) && (Log() << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl);
    MolecularWalker->Leaf->CyclicStructureAnalysis(LocalBackEdgeStack, MinimumRingSize);
    DoLog(0) && (Log() << Verbose(0) << "Done with Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl);
    delete(LocalBackEdgeStack);
    delete(ListOfAtoms);
    FragmentCounter++;
  }
  delete(BackEdgeStack);

  // ===== 3. if structure still valid, parse key set file and others =====
  FragmentationToDo = FragmentationToDo && ParseKeySetFile(prefix, ParsedFragmentList);

  // ===== 4. check globally whether there's something to do actually (first adaptivity check)
  FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(prefix);

  // =================================== Begin of FRAGMENTATION ===============================
  // ===== 6a. assign each keyset to its respective subgraph =====
  ListOfLocalAtoms = new atom **[LeafCount];
  for (int i=0;i<LeafCount;i++)
    ListOfLocalAtoms[i] = NULL;
  FragmentCounter = 0;
  Subgraphs->next->AssignKeySetsToFragment(this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, FragmentCounter, true);
  delete[](ListOfLocalAtoms);

  // ===== 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[getAtomCount()+1];
  AtomMask[getAtomCount()] = false;
  FragmentationToDo = false;  // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards
  while ((CheckOrder = CheckOrderAtSite(AtomMask, ParsedFragmentList, Order, MinimumRingSize, prefix))) {
    FragmentationToDo = FragmentationToDo || CheckOrder;
    AtomMask[getAtomCount()] = 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(RootStack, AtomMask, (FragmentCounter = 0));

    // ===== 7. fill the bond fragment list =====
    FragmentCounter = 0;
    MolecularWalker = Subgraphs;
    while (MolecularWalker->next != NULL) {
      MolecularWalker = MolecularWalker->next;
      DoLog(1) && (Log() << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl);
      //MolecularWalker->Leaf->OutputListOfBonds(out);  // output atom::ListOfBonds for debugging
      if (MolecularWalker->Leaf->hasBondStructure()) {
        // call BOSSANOVA method
        DoLog(0) && (Log() << Verbose(0) << endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= " << endl);
        MolecularWalker->Leaf->FragmentBOSSANOVA(FragmentList[FragmentCounter], RootStack[FragmentCounter], MinimumRingSize);
      } else {
        DoeLog(1) && (eLog()<< Verbose(1) << "Subgraph " << MolecularWalker << " has no atoms!" << endl);
      }
      FragmentCounter++;  // next fragment list
    }
  }
  DoLog(2) && (Log() << Verbose(2) << "CheckOrder is " << CheckOrder << "." << endl);
  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(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);
  }
  delete[](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(World::getPointer());
  int k=0;
  for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
    KeySet test = (*runner).first;
    DoLog(0) && (Log() << Verbose(0) << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl);
    BondFragments->insert(StoreFragmentFromKeySet(test, World::getInstance().getConfig()));
    k++;
  }
  DoLog(0) && (Log() << Verbose(0) << 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
    int *SortIndex = NULL;
    CreateMappingLabelsToConfigSequence(SortIndex);

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

    // store force index reference file
    BondFragments->StoreForcesFile(prefix, SortIndex);

    // store keysets file
    StoreKeySetFile(TotalGraph, prefix);

    {
      // store Adjacency file
      std::string filename = prefix + ADJACENCYFILE;
      StoreAdjacencyToFile(filename);
    }

    // store Hydrogen saturation correction file
    BondFragments->AddHydrogenCorrection(prefix);

    // store adaptive orders into file
    StoreOrderAtSiteFile(prefix);

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

    // free memory for bond part
    DoLog(1) && (Log() << Verbose(1) << "Freeing bond memory" << endl);
    delete[](SortIndex);
  } else {
    DoLog(1) && (Log() << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl);
  }
  delete(BondFragments);
  DoLog(0) && (Log() << 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 &path path to file ORDERATSITEFILE
 * \return true - file writable, false - not writable
 */
bool molecule::StoreOrderAtSiteFile(std::string &path)
{
  string line;
  ofstream file;

  line = path + ORDERATSITEFILE;
  file.open(line.c_str());
  DoLog(1) && (Log() << Verbose(1) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl);
  if (file.good()) {
    for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputOrder), &file));
    file.close();
    DoLog(1) && (Log() << Verbose(1) << "done." << endl);
    return true;
  } else {
    DoLog(1) && (Log() << Verbose(1) << "failed to open file " << line << "." << 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 &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(std::string &path)
{
  unsigned char *OrderArray = new unsigned char[getAtomCount()];
  bool *MaxArray = new bool[getAtomCount()];
  bool status;
  int AtomNr, value;
  string line;
  ifstream file;

  for(int i=0;i<getAtomCount();i++) {
    OrderArray[i] = 0;
    MaxArray[i] = false;
  }

  DoLog(1) && (Log() << Verbose(1) << "Begin of ParseOrderAtSiteFromFile" << endl);
  line = path + ORDERATSITEFILE;
  file.open(line.c_str());
  if (file.good()) {
    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;
        //Log() << Verbose(2) << "AtomNr " << AtomNr << " with order " << (int)OrderArray[AtomNr] << " and max order set to " << (int)MaxArray[AtomNr] << "." << endl;
      }
    }
    file.close();

    // set atom values
    for(internal_iterator iter=atoms.begin();iter!=atoms.end();++iter){
      (*iter)->AdaptiveOrder = OrderArray[(*iter)->nr];
      (*iter)->MaxOrder = MaxArray[(*iter)->nr];
    }
    //SetAtomValueToIndexedArray( OrderArray, &atom::nr, &atom::AdaptiveOrder );
    //SetAtomValueToIndexedArray( MaxArray, &atom::nr, &atom::MaxOrder );

    DoLog(1) && (Log() << Verbose(1) << "\t ... done." << endl);
    status = true;
  } else {
    DoLog(1) && (Log() << Verbose(1) << "\t ... failed to open file " << line << "." << endl);
    status = false;
  }
  delete[](OrderArray);
  delete[](MaxArray);

  DoLog(1) && (Log() << 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(KeySet *&Leaf, int *&ShortestPathList)
{
  atom *Runner = NULL;
  int SP, Removal;

  DoLog(2) && (Log() << 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->getType()->getAtomicNumber() != 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 calloc'd 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;

  // 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(molecule *mol, molecule *Leaf, atom **SonList, bool IsAngstroem)
{
  bool LonelyFlag = false;
  atom *OtherFather = NULL;
  atom *FatherOfRunner = NULL;

#ifdef ADDHYDROGEN
  molecule::const_iterator runner;
#endif
  // we increment the iter just before skipping the hydrogen
  for (molecule::const_iterator iter = Leaf->begin(); iter != Leaf->end();) {
    LonelyFlag = true;
    FatherOfRunner = (*iter)->father;
    ASSERT(FatherOfRunner,"Atom without father found");
    if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
      // create all bonds
      for (BondList::const_iterator BondRunner = FatherOfRunner->ListOfBonds.begin(); BondRunner != FatherOfRunner->ListOfBonds.end(); (++BondRunner)) {
        OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner);
//        Log() << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->nr] << " is bound to " << *OtherFather;
        if (SonList[OtherFather->nr] != NULL) {
//          Log() << Verbose(0) << ", 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)
//            Log() << Verbose(3) << "Adding Bond: ";
//            Log() << Verbose(0) <<
            Leaf->AddBond((*iter), SonList[OtherFather->nr], (*BondRunner)->BondDegree);
//            Log() << Verbose(0) << "." << endl;
            //NumBonds[(*iter)->nr]++;
          } else {
//            Log() << Verbose(3) << "Not adding bond, labels in wrong order." << endl;
          }
          LonelyFlag = false;
        } else {
//          Log() << Verbose(0) << ", who has no son in this fragment molecule." << endl;
#ifdef ADDHYDROGEN
          //Log() << Verbose(3) << "Adding Hydrogen to " << (*iter)->Name << " and a bond in between." << endl;
          if(!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem))
            exit(1);
#endif
          //NumBonds[(*iter)->nr] += Binder->BondDegree;
        }
      }
    } else {
    DoeLog(1) && (eLog()<< Verbose(1) << "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl);
    }
    if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) {
      DoLog(0) && (Log() << Verbose(0) << **iter << "has got bonds only to hydrogens!" << endl);
    }
    ++iter;
#ifdef ADDHYDROGEN
    while ((iter != Leaf->end()) && ((*iter)->getType()->getAtomicNumber() == 1)){ // skip added hydrogen
      iter++;
    }
#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(KeySet &Leaflet, bool IsAngstroem)
{
  atom **SonList = new atom*[getAtomCount()];
  molecule *Leaf = World::getInstance().createMolecule();

  for(int i=0;i<getAtomCount();i++)
    SonList[i] = NULL;

//  Log() << 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
//  Log() << Verbose(2) << "Creating bonds from father graph (i.e. induced subgraph creation)." << endl;
  CreateInducedSubgraphOfFragment(this, Leaf, SonList, IsAngstroem);

  //Leaflet->Leaf->ScanForPeriodicCorrection(out);
  delete[](SonList);
//  Log() << Verbose(1) << "End of StoreFragmentFromKeyset." << endl;
  return Leaf;
};


/** Clears the touched list
 * \param *out output stream for debugging
 * \param verbosity verbosity level
 * \param *&TouchedList touched list
 * \param SubOrder current suborder
 * \param TouchedIndex currently touched
 */
void SPFragmentGenerator_ClearingTouched(int verbosity, int *&TouchedList, int SubOrder, int &TouchedIndex)
{
  Log() << Verbose(1+verbosity) << "Clearing touched list." << endl;
  for (TouchedIndex=SubOrder+1;TouchedIndex--;)  // empty touched list
    TouchedList[TouchedIndex] = -1;
  TouchedIndex = 0;

}

/** Adds the current combination of the power set to the snake stack.
 * \param *out output stream for debugging
 * \param verbosity verbosity level
 * \param CurrentCombination
 * \param SetDimension maximum number of bits in power set
 * \param *FragmentSet snake stack to remove from
 * \param *&TouchedList touched list
 * \param TouchedIndex currently touched
 * \return number of set bits
 */
int AddPowersetToSnakeStack(int verbosity, int CurrentCombination, int SetDimension, KeySet *FragmentSet, bond **BondsSet, int *&TouchedList, int &TouchedIndex)
{
  atom *OtherWalker = NULL;
  bool bit = false;
  KeySetTestPair TestKeySetInsert;

  int Added = 0;
  for (int j=0;j<SetDimension;j++) {  // pull out every bit by shifting
    bit = ((CurrentCombination & (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
      //Log() << Verbose(1+verbosity) << "Current Bond is " << BondsSet[j] << ", checking on " << *OtherWalker << "." << endl;
      Log() << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
      TestKeySetInsert = FragmentSet->insert(OtherWalker->nr);
      if (TestKeySetInsert.second) {
        TouchedList[TouchedIndex++] = OtherWalker->nr;  // note as added
        Added++;
      } else {
        Log() << Verbose(2+verbosity) << "This was item was already present in the keyset." << endl;
      }
    } else {
      Log() << Verbose(2+verbosity) << "Not adding." << endl;
    }
  }
  return Added;
};

/** Counts the number of elements in a power set.
 * \param *SetFirst
 * \param *SetLast
 * \param *&TouchedList touched list
 * \param TouchedIndex currently touched
 * \return number of elements
 */
int CountSetMembers(bond *SetFirst, bond *SetLast, int *&TouchedList, int TouchedIndex)
{
  int SetDimension = 0;
  bond *Binder = SetFirst;    // start node for this level
  while (Binder->next != SetLast) {    // 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
        SetDimension++;
    }
  }
  return SetDimension;
};

/** Counts the number of elements in a power set.
 * \param *BondsList bonds list to fill
 * \param *SetFirst
 * \param *SetLast
 * \param *&TouchedList touched list
 * \param TouchedIndex currently touched
 * \return number of elements
 */
int FillBondsList(bond **BondsList, bond *SetFirst, bond *SetLast, int *&TouchedList, int TouchedIndex)
{
  int SetDimension = 0;
  bond *Binder = SetFirst;    // start node for this level
  while (Binder->next != SetLast) {    // compare to end node of this level
    Binder = Binder->next;
    for (int k=0;k<TouchedIndex;k++) {
      if (Binder->leftatom->nr == TouchedList[k])   // leftatom is always the close one
        BondsList[SetDimension++] = Binder;
    }
  }
  return SetDimension;
};

/** Remove all items that were added on this SP level.
 * \param *out output stream for debugging
 * \param verbosity verbosity level
 * \param *FragmentSet snake stack to remove from
 * \param *&TouchedList touched list
 * \param TouchedIndex currently touched
 */
void RemoveAllTouchedFromSnakeStack(int verbosity, KeySet *FragmentSet, int *&TouchedList, int &TouchedIndex)
{
  int Removal = 0;
  for(int j=0;j<TouchedIndex;j++) {
    Removal = TouchedList[j];
    Log() << Verbose(2+verbosity) << "Removing item nr. " << Removal << " from snake stack." << endl;
    FragmentSet->erase(Removal);
    TouchedList[j] = -1;
  }
  DoLog(2) && (Log() << Verbose(2) << "Remaining local nr.s on snake stack are: ");
  for(KeySet::iterator runner = FragmentSet->begin(); runner != FragmentSet->end(); runner++)
    DoLog(0) && (Log() << Verbose(0) << (*runner) << " ");
  DoLog(0) && (Log() << Verbose(0) << endl);
  TouchedIndex = 0; // set Index to 0 for list of atoms added on this level
};

/** 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(struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder)
{
  int verbosity = 0; //FragmentSearch->ANOVAOrder-SubOrder;
  int NumCombinations;
  int bits, TouchedIndex, SubSetDimension, SP, Added;
  int SpaceLeft;
  int *TouchedList = new int[SubOrder + 1];
  KeySetTestPair TestKeySetInsert;

  NumCombinations = 1 << SetDimension;

  // here for all bonds of Walker all combinations of end pieces (from the bonds)
  // have to be added and for the remaining ANOVA order GraphCrawler be called
  // recursively for the next level

  Log() << Verbose(1+verbosity) << "Begin of SPFragmentGenerator." << endl;
  Log() << 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)
  SPFragmentGenerator_ClearingTouched(verbosity, TouchedList, SubOrder, TouchedIndex);

  // create every possible combination of the endpieces
  Log() << 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;

    Log() << 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 = AddPowersetToSnakeStack(verbosity, i, SetDimension, FragmentSearch->FragmentSet, BondsSet, TouchedList, TouchedIndex);

      SpaceLeft = SubOrder - Added ;// SubOrder - bits; // due to item's maybe being already present, this does not work anymore
      if (SpaceLeft > 0) {
        Log() << 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 = CountSetMembers(FragmentSearch->BondsPerSPList[2*SP], FragmentSearch->BondsPerSPList[2*SP+1], TouchedList, TouchedIndex);

          // then allocate and fill the list
          bond *BondsList[SubSetDimension];
          SubSetDimension = FillBondsList(BondsList, FragmentSearch->BondsPerSPList[2*SP], FragmentSearch->BondsPerSPList[2*SP+1], TouchedList, TouchedIndex);

          // then iterate
          Log() << Verbose(2+verbosity) << "Calling subset generator " << SP << " away from root " << *FragmentSearch->Root << " with sub set dimension " << SubSetDimension << "." << endl;
          SPFragmentGenerator(FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits);
        }
      } else {
        // --2-- otherwise store the complete fragment
        Log() << Verbose(1+verbosity) << "Enough items on stack for a fragment!" << endl;
        // store fragment as a KeySet
        DoLog(2) && (Log() << Verbose(2) << "Found a new fragment[" << FragmentSearch->FragmentCounter << "], local nr.s are: ");
        for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
          DoLog(0) && (Log() << Verbose(0) << (*runner) << " ");
        DoLog(0) && (Log() << Verbose(0) << endl);
        //if (!CheckForConnectedSubgraph(FragmentSearch->FragmentSet))
          //DoeLog(1) && (eLog()<< Verbose(1) << "The found fragment is not a connected subgraph!" << endl);
        InsertFragmentIntoGraph(FragmentSearch);
      }

      // --3-- remove all added items in this level from snake stack
      Log() << Verbose(1+verbosity) << "Removing all items that were added on this SP level " << RootDistance << "." << endl;
      RemoveAllTouchedFromSnakeStack(verbosity, FragmentSearch->FragmentSet, TouchedList, TouchedIndex);
    } else {
      Log() << Verbose(2+verbosity) << "More atoms to add for this set (" << bits << ") than space left on stack " << SubOrder << ", skipping this set." << endl;
    }
  }
  delete[](TouchedList);
  Log() << Verbose(1+verbosity) << "End of SPFragmentGenerator, " << RootDistance << " away from Root " << *FragmentSearch->Root << " and SubOrder is " << SubOrder << "." << endl;
};

/** Allocates memory for UniqueFragments::BondsPerSPList.
 * \param *out output stream
 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
 * \param FragmentSearch UniqueFragments
 * \sa FreeSPList()
 */
void InitialiseSPList(int Order, struct UniqueFragments &FragmentSearch)
{
  FragmentSearch.BondsPerSPList = new bond* [Order * 2];
  FragmentSearch.BondsPerSPCount = new int[Order];
  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;
  }
};

/** Free's memory for for UniqueFragments::BondsPerSPList.
 * \param *out output stream
 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
 * \param FragmentSearch UniqueFragments\
 * \sa InitialiseSPList()
 */
void FreeSPList(int Order, struct UniqueFragments &FragmentSearch)
{
  delete[](FragmentSearch.BondsPerSPCount);
  for (int i=Order;i--;) {
    delete(FragmentSearch.BondsPerSPList[2*i]);
    delete(FragmentSearch.BondsPerSPList[2*i+1]);
  }
  delete[](FragmentSearch.BondsPerSPList);
};

/** Sets FragmenSearch to initial value.
 * Sets UniqueFragments::ShortestPathList entries to zero, UniqueFragments::BondsPerSPCount to zero (except zero level to 1) and
 * adds initial bond UniqueFragments::Root to UniqueFragments::Root to UniqueFragments::BondsPerSPList
 * \param *out output stream
 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
 * \param FragmentSearch UniqueFragments
 * \sa FreeSPList()
 */
void SetSPList(int Order, struct UniqueFragments &FragmentSearch)
{
  // 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=Order;i--;)
    FragmentSearch.BondsPerSPCount[i] = 0;
  FragmentSearch.BondsPerSPCount[0] = 1;
  bond *Binder = new bond(FragmentSearch.Root, FragmentSearch.Root);
  add(Binder, FragmentSearch.BondsPerSPList[1]);
};

/** Resets UniqueFragments::ShortestPathList and cleans bonds from UniqueFragments::BondsPerSPList.
 * \param *out output stream
 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
 * \param FragmentSearch UniqueFragments
 * \sa InitialiseSPList()
 */
void ResetSPList(int Order, struct UniqueFragments &FragmentSearch)
{
  bond *Binder = NULL;
  DoLog(0) && (Log() << Verbose(0) << "Free'ing all found lists. and resetting index lists" << endl);
  for(int i=Order;i--;) {
    DoLog(1) && (Log() << Verbose(1) << "Current SP level is " << i << ": ");
    Binder = FragmentSearch.BondsPerSPList[2*i];
    while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
      Binder = Binder->next;
      // Log() << Verbose(0) << "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
    DoLog(0) && (Log() << Verbose(0) << "cleaned." << endl);
  }
};


/** Fills the Bonds per Shortest Path List and set the vertex labels.
 * \param *out output stream
 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
 * \param FragmentSearch UniqueFragments
 * \param *mol molecule with atoms and bonds
 * \param RestrictedKeySet Restricted vertex set to use in context of molecule
 */
void FillSPListandLabelVertices(int Order, struct UniqueFragments &FragmentSearch, molecule *mol, KeySet RestrictedKeySet)
{
  // 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.
  int AtomKeyNr = -1;
  atom *Walker = NULL;
  atom *OtherWalker = NULL;
  atom *Predecessor = NULL;
  bond *CurrentEdge = NULL;
  bond *Binder = NULL;
  int RootKeyNr = FragmentSearch.Root->GetTrueFather()->nr;
  int RemainingWalkers = -1;
  int SP = -1;

  DoLog(0) && (Log() << Verbose(0) << "Starting BFS analysis ..." << endl);
  for (SP = 0; SP < (Order-1); SP++) {
    DoLog(1) && (Log() << Verbose(1) << "New SP level reached: " << SP << ", creating new SP list with " << FragmentSearch.BondsPerSPCount[SP] << " item(s)");
    if (SP > 0) {
      DoLog(0) && (Log() << Verbose(0) << ", old level closed with " << FragmentSearch.BondsPerSPCount[SP-1] << " item(s)." << endl);
      FragmentSearch.BondsPerSPCount[SP] = 0;
    } else
      DoLog(0) && (Log() << Verbose(0) << "." << 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;
      DoLog(0) && (Log() << 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
      DoLog(1) && (Log() << Verbose(1) << "Going through all bonds of Walker." << endl);
      for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
        OtherWalker = (*Runner)->GetOtherAtom(Walker);
        if ((RestrictedKeySet.find(OtherWalker->nr) != RestrictedKeySet.end())
  #ifdef ADDHYDROGEN
         && (OtherWalker->getType()->getAtomicNumber() != 1)
  #endif
                                                              ) {  // skip hydrogens and restrict to fragment
          DoLog(2) && (Log() << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *(*Runner) << "." << 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;
            DoLog(3) && (Log() << 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]++;
            DoLog(3) && (Log() << Verbose(3) << "Added its bond to SP list, having now " << FragmentSearch.BondsPerSPCount[SP+1] << " item(s)." << endl);
          } else {
            if (OtherWalker != Predecessor)
              DoLog(3) && (Log() << Verbose(3) << "Not passing on, as index of " << *OtherWalker << " " << OtherWalker->GetTrueFather()->nr << " is smaller than that of Root " << RootKeyNr << "." << endl);
            else
              DoLog(3) && (Log() << Verbose(3) << "This is my predecessor " << *Predecessor << "." << endl);
          }
        } else Log() << Verbose(2) << "Is not in the restricted keyset or skipping hydrogen " << *OtherWalker << "." << endl;
      }
    }
  }
};

/** prints the Bonds per Shortest Path list in UniqueFragments.
 * \param *out output stream
 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
 * \param FragmentSearch UniqueFragments
 */
void OutputSPList(int Order, struct UniqueFragments &FragmentSearch)
{
  bond *Binder = NULL;
  DoLog(0) && (Log() << 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];
    DoLog(1) && (Log() << Verbose(1) << "Current SP level is " << i << "." << endl);
    while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
      Binder = Binder->next;
      DoLog(2) && (Log() << Verbose(2) << *Binder << endl);
    }
  }
};

/** Simply counts all bonds in all UniqueFragments::BondsPerSPList lists.
 * \param *out output stream
 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
 * \param FragmentSearch UniqueFragments
 */
int CountNumbersInBondsList(int Order, struct UniqueFragments &FragmentSearch)
{
  bond *Binder = NULL;
  int 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++;
    }
  }
  return SP;
};

/** 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(int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet)
{
  int Counter = FragmentSearch.FragmentCounter; // mark current value of counter

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

  SetSPList(Order, FragmentSearch);

  // do a BFS search to fill the SP lists and label the found vertices
  FillSPListandLabelVertices(Order, FragmentSearch, this, RestrictedKeySet);

  // outputting all list for debugging
  OutputSPList(Order, FragmentSearch);

  // creating fragments with the found edge sets  (may be done in reverse order, faster)
  int SP = CountNumbersInBondsList(Order, FragmentSearch);
  DoLog(0) && (Log() << Verbose(0) << "Total number of edges is " << SP << "." << endl);
  if (SP >= (Order-1)) {
    // start with root (push on fragment stack)
    DoLog(0) && (Log() << Verbose(0) << "Starting fragment generation with " << *FragmentSearch.Root << ", local nr is " << FragmentSearch.Root->nr << "." << endl);
    FragmentSearch.FragmentSet->clear();
    DoLog(0) && (Log() << Verbose(0) << "Preparing subset for this root and calling generator." << endl);

    // prepare the subset and call the generator
    bond* BondsList[FragmentSearch.BondsPerSPCount[0]];
    for(int i=0;i<FragmentSearch.BondsPerSPCount[0];i++)
      BondsList[i] = NULL;
    BondsList[0] = FragmentSearch.BondsPerSPList[0]->next;  // on SP level 0 there's only the root bond

    SPFragmentGenerator(&FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order);
  } else {
    DoLog(0) && (Log() << 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
  DoLog(0) && (Log() << Verbose(0) << "Removing root again from stack." << endl);
  FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);

  // free'ing the bonds lists
  ResetSPList(Order, FragmentSearch);

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

bool KeyCompare::operator() (const KeySet SubgraphA, const KeySet SubgraphB) const
{
  //Log() << Verbose(0) << "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
};


/** Combines all KeySets from all orders into single ones (with just unique entries).
 * \param *out output stream for debugging
 * \param *&FragmentList list to fill
 * \param ***FragmentLowerOrdersList
 * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
 * \param *mol molecule with atoms and bonds
 */
int CombineAllOrderListIntoOne(Graph *&FragmentList, Graph ***FragmentLowerOrdersList, KeyStack &RootStack, molecule *mol)
{
  int RootNr = 0;
  int RootKeyNr = 0;
  int StartNr = 0;
  int counter = 0;
  int NumLevels = 0;
  atom *Walker = NULL;

  DoLog(0) && (Log() << 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();
  }

  StartNr = RootStack.back();
  do {
    RootKeyNr = RootStack.front();
    RootStack.pop_front();
    Walker = mol->FindAtom(RootKeyNr);
    NumLevels = 1 << (Walker->AdaptiveOrder - 1);
    for(int i=0;i<NumLevels;i++) {
      if (FragmentLowerOrdersList[RootNr][i] != NULL) {
        InsertGraphIntoGraph(*FragmentList, (*FragmentLowerOrdersList[RootNr][i]), &counter);
      }
    }
    RootStack.push_back(Walker->nr);
    RootNr++;
  } while (RootKeyNr != StartNr);
  return counter;
};

/** Free's memory allocated for all KeySets from all orders.
 * \param *out output stream for debugging
 * \param ***FragmentLowerOrdersList
 * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
 * \param *mol molecule with atoms and bonds
 */
void FreeAllOrdersList(Graph ***FragmentLowerOrdersList, KeyStack &RootStack, molecule *mol)
{
  DoLog(1) && (Log() << Verbose(1) << "Free'ing the lists of all orders per order." << endl);
  int RootNr = 0;
  int RootKeyNr = 0;
  int NumLevels = 0;
  atom *Walker = NULL;
  while (!RootStack.empty()) {
    RootKeyNr = RootStack.front();
    RootStack.pop_front();
    Walker = mol->FindAtom(RootKeyNr);
    NumLevels = 1 << (Walker->AdaptiveOrder - 1);
    for(int i=0;i<NumLevels;i++) {
      if (FragmentLowerOrdersList[RootNr][i] != NULL) {
        delete(FragmentLowerOrdersList[RootNr][i]);
      }
    }
    delete[](FragmentLowerOrdersList[RootNr]);
    RootNr++;
  }
  delete[](FragmentLowerOrdersList);
};


/** 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(Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize)
{
  Graph ***FragmentLowerOrdersList = NULL;
  int NumLevels = 0;
  int NumMolecules = 0;
  int TotalNumMolecules = 0;
  int *NumMoleculesOfOrder = NULL;
  int Order = 0;
  int UpgradeCount = RootStack.size();
  KeyStack FragmentRootStack;
  int RootKeyNr = 0;
  int RootNr = 0;
  struct UniqueFragments FragmentSearch;

  DoLog(0) && (Log() << 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 = new int[UpgradeCount];
  FragmentLowerOrdersList = new Graph**[UpgradeCount];

  for(int i=0;i<UpgradeCount;i++) {
    NumMoleculesOfOrder[i] = 0;
    FragmentLowerOrdersList[i] = NULL;
  }

  // initialise the fragments structure
  FragmentSearch.FragmentCounter = 0;
  FragmentSearch.FragmentSet = new KeySet;
  FragmentSearch.Root = FindAtom(RootKeyNr);
  FragmentSearch.ShortestPathList = new int[getAtomCount()];
  for (int i=getAtomCount();i--;) {
    FragmentSearch.ShortestPathList[i] = -1;
  }

  // Construct the complete KeySet which we need for topmost level only (but for all Roots)
  KeySet CompleteMolecule;
  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    CompleteMolecule.insert((*iter)->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();
    atom *Walker = FindAtom(RootKeyNr);
    // check cyclic lengths
    //if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) {
    //  Log() << 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
      InitialiseSPList(Order, FragmentSearch);

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

      // create top order where nothing is reduced
      DoLog(0) && (Log() << Verbose(0) << "==============================================================================================================" << endl);
      DoLog(0) && (Log() << 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(Walker->AdaptiveOrder, FragmentSearch, CompleteMolecule);

      // output resulting number
      DoLog(1) && (Log() << Verbose(1) << "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << "." << endl);
      if (NumMoleculesOfOrder[RootNr] != 0) {
        NumMolecules = 0;
      } else {
        Walker->GetTrueFather()->MaxOrder = true;
      }
      // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
      //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
      TotalNumMolecules += NumMoleculesOfOrder[RootNr];
//      Log() << 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
      FreeSPList(Order, FragmentSearch);
    }
  }
  DoLog(0) && (Log() << Verbose(0) << "==============================================================================================================" << endl);
  DoLog(1) && (Log() << Verbose(1) << "Total number of resulting molecules is: " << TotalNumMolecules << "." << endl);
  DoLog(0) && (Log() << Verbose(0) << "==============================================================================================================" << endl);

  // cleanup FragmentSearch structure
  delete[](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)
  CombineAllOrderListIntoOne(FragmentList, FragmentLowerOrdersList, RootStack, this);
  FreeAllOrdersList(FragmentLowerOrdersList, RootStack, this);
  delete[](NumMoleculesOfOrder);

  DoLog(0) && (Log() << 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()
{
  bond *Binder = NULL;
  bond *OtherBinder = NULL;
  atom *Walker = NULL;
  atom *OtherWalker = NULL;
  Matrix matrix = World::getInstance().getDomain().getM();
  enum Shading *ColorList = NULL;
  double tmp;
  Vector Translationvector;
  //class StackClass<atom *> *CompStack = NULL;
  class StackClass<atom *> *AtomStack = new StackClass<atom *>(getAtomCount());
  bool flag = true;

  DoLog(2) && (Log() << Verbose(2) << "Begin of ScanForPeriodicCorrection." << endl);

  ColorList = new enum Shading[getAtomCount()];
  for (int i=0;i<getAtomCount();i++)
    ColorList[i] = (enum Shading)0;
  while (flag) {
    // remove bonds that are beyond bonddistance
    Translationvector.Zero();
    // scan all bonds
    flag = false;
    for(molecule::iterator AtomRunner = begin(); (!flag) && (AtomRunner != end()); ++AtomRunner)
      for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); (!flag) && (BondRunner != (*AtomRunner)->ListOfBonds.end()); ++BondRunner) {
        Binder = (*BondRunner);
        for (int i=NDIM;i--;) {
          tmp = fabs(Binder->leftatom->at(i) - Binder->rightatom->at(i));
          //Log() << 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
            DoLog(2) && (Log() << 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->at(i) - Binder->rightatom->at(i);
        if (fabs(tmp) > BondDistance)
          Translationvector[i] = (tmp < 0) ? +1. : -1.;
      }
      Translationvector *= matrix;
      //Log() << Verbose(3) << "Translation vector is ";
      Log() << Verbose(0) << Translationvector <<  endl;
      // apply to all atoms of first component via BFS
      for (int i=getAtomCount();i--;)
        ColorList[i] = white;
      AtomStack->Push(Binder->leftatom);
      while (!AtomStack->IsEmpty()) {
        Walker = AtomStack->PopFirst();
        //Log() << Verbose (3) << "Current Walker is: " << *Walker << "." << endl;
        ColorList[Walker->nr] = black;    // mark as explored
        *Walker += Translationvector; // translate
        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
          if ((*Runner) != Binder) {
            OtherWalker = (*Runner)->GetOtherAtom(Walker);
            if (ColorList[OtherWalker->nr] == white) {
              AtomStack->Push(OtherWalker); // push if yet unexplored
            }
          }
        }
      }
      // re-add bond
      link(Binder, OtherBinder);
    } else {
      DoLog(3) && (Log() << Verbose(3) << "No corrections for this fragment." << endl);
    }
    //delete(CompStack);
  }
  // free allocated space from ReturnFullMatrixforSymmetric()
  delete(AtomStack);
  delete[](ColorList);
  DoLog(2) && (Log() << Verbose(2) << "End of ScanForPeriodicCorrection." << endl);
};