source: src/molecule_fragmentation.cpp@ f67817f

/*
 * 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 "CodePatterns/MemDebug.hpp"

#include <cstring>

#include "atom.hpp"
#include "Bond/bond.hpp"
#include "Box.hpp"
#include "CodePatterns/Verbose.hpp"
#include "CodePatterns/Log.hpp"
#include "config.hpp"
#include "Element/element.hpp"
#include "Element/periodentafel.hpp"
#include "Fragmentation/fragmentation_helpers.hpp"
#include "Fragmentation/PowerSetGenerator.hpp"
#include "Fragmentation/UniqueFragments.hpp"
#include "Graph/BondGraph.hpp"
#include "Graph/CheckAgainstAdjacencyFile.hpp"
#include "Graph/CyclicStructureAnalysis.hpp"
#include "Graph/DepthFirstSearchAnalysis.hpp"
#include "Helpers/helpers.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "molecule.hpp"
#include "World.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) {
    const BondList& ListOfBonds = (*iter)->getListOfBonds();
    c = (ListOfBonds.size() > c) ? 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;
};

/** 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 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, 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)->getNr()] = 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)->getNr()] = true;  // include all (non-hydrogen) atoms
        if ((Order != 0) && ((*iter)->AdaptiveOrder < Order)) // && ((*iter)->AdaptiveOrder < MinimumRingSize[(*iter)->getNr()]))
          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)->getNr()]==-1,"Same SortIndex set twice");
    SortIndex[(*iter)->getNr()] = 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()->getNr()) ? (*iter)->GetTrueFather()->getNr() : 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()->getNr();
      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
 * \param &DFS contains bond structure analysis data
 * \return 1 - continue, 2 - stop (no fragmentation occured)
 */
int molecule::FragmentMolecule(int Order, std::string &prefix, DepthFirstSearchAnalysis &DFS)
{
  MoleculeListClass *BondFragments = NULL;
  int FragmentCounter;
  MoleculeLeafClass *MolecularWalker = NULL;
  MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
  fstream File;
  bool FragmentationToDo = true;
  bool CheckOrder = false;
  Graph **FragmentList = NULL;
  Graph *ParsedFragmentList = NULL;
  Graph TotalGraph;     // graph with all keysets however local numbers
  int TotalNumberOfKeySets = 0;
  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 ====

  // === compare it with adjacency file ===
  {
    std::ifstream File;
    std::string filename;
    filename = prefix + ADJACENCYFILE;
    File.open(filename.c_str(), ios::out);
    DoLog(1) && (Log() << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... " << endl);

    CheckAgainstAdjacencyFile FileChecker(World::getInstance().beginAtomSelection(), World::getInstance().endAtomSelection());
    FragmentationToDo = FragmentationToDo && FileChecker(File);
  }

  // === reset bond degree and perform CorrectBondDegree ===
  for(World::MoleculeIterator iter = World::getInstance().getMoleculeIter();
      iter != World::getInstance().moleculeEnd();
      ++iter) {
    // correct bond degree
    World::AtomComposite Set = (*iter)->getAtomSet();
    World::getInstance().getBondGraph()->CorrectBondDegree(Set);
  }

  // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
  // NOTE: We assume here that DFS has been performed and molecule structure is current.
  Subgraphs = DFS.getMoleculeStructure();

  // ===== 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 =====
  const int MolCount = World::getInstance().numMolecules();
  ListOfLocalAtoms = new atom **[MolCount];
  for (int i=0;i<MolCount;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, 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);
      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]);
      } 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);

  // ==================================== 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;
  while (Subgraphs != NULL) {
    // remove entry in fragment list
    // remove subgraph fragment
    MolecularWalker = Subgraphs->next;
    Subgraphs->Leaf = NULL;
    delete(Subgraphs);
    Subgraphs = MolecularWalker;
  }
  // free fragment list
  for (int i=0; i< FragmentCounter; ++i )
    delete(FragmentList[i]);
  delete[](FragmentList);

  DoLog(0) && (Log() << Verbose(0) << FragmentCounter << " subgraph fragments have been removed." << std::endl);

  // ===== 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);
  }
  // remove all create molecules again from the World including their atoms
  for (MoleculeList::iterator iter = BondFragments->ListOfMolecules.begin();
      !BondFragments->ListOfMolecules.empty();
      iter = BondFragments->ListOfMolecules.begin()) {
    // remove copied atoms and molecule again
    molecule *mol = *iter;
    mol->removeAtomsinMolecule();
    World::getInstance().destroyMolecule(mol);
    BondFragments->ListOfMolecules.erase(iter);
  }
  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)->getNr()];
      (*iter)->MaxOrder = MaxArray[(*iter)->getNr()];
    }
    //SetAtomValueToIndexedArray( OrderArray, &atom::getNr(), &atom::AdaptiveOrder );
    //SetAtomValueToIndexedArray( MaxArray, &atom::getNr(), &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 std::deque<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;

  // 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->getNr()] = 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->getNr()] != NULL)  {  // check if this, our father, is present in list
      // create all bonds
      const BondList& ListOfBonds = FatherOfRunner->getListOfBonds();
      for (BondList::const_iterator BondRunner = ListOfBonds.begin();
          BondRunner != ListOfBonds.end();
          ++BondRunner) {
        OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner);
//        Log() << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()] << " is bound to " << *OtherFather;
        if (SonList[OtherFather->getNr()] != NULL) {
//          Log() << Verbose(0) << ", whose son is " << *SonList[OtherFather->getNr()] << "." << endl;
          if (OtherFather->getNr() > FatherOfRunner->getNr()) { // 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->getNr()], (*BondRunner)->BondDegree);
//            Log() << Verbose(0) << "." << endl;
            //NumBonds[(*iter)->getNr()]++;
          } 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)->getNr()] += Binder->BondDegree;
        }
      }
    } else {
    DoeLog(1) && (eLog()<< Verbose(1) << "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->getNr()] << "!" << 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;
};

/** 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)
 * \return pointer to Graph list
 */
void molecule::FragmentBOSSANOVA(Graph *&FragmentList, KeyStack &RootStack)
{
  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;
  class PowerSetGenerator PSG;

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

  FragmentSearch.Init(FindAtom(RootKeyNr), getAtomCount());

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

  // 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()->getNr()] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->getNr()])) {
    //  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
      FragmentSearch.InitialiseSPList(Order);

      // 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.setRoot(Walker);
      NumMoleculesOfOrder[RootNr] = PSG(&FragmentSearch, Walker->AdaptiveOrder, 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
      FragmentSearch.FreeSPList(Order);
    }
  }
  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
  FragmentSearch.Cleanup();

  // 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
 */
bool molecule::ScanForPeriodicCorrection()
{
  bond *Binder = NULL;
  //bond *OtherBinder = NULL;
  atom *Walker = NULL;
  atom *OtherWalker = NULL;
  RealSpaceMatrix matrix = World::getInstance().getDomain().getM();
  enum GraphEdge::Shading *ColorList = NULL;
  double tmp;
  //bool LastBond = true; // only needed to due list construct
  Vector Translationvector;
  //std::deque<atom *> *CompStack = NULL;
  std::deque<atom *> *AtomStack = new std::deque<atom *>; // (getAtomCount());
  bool flag = true;
  BondGraph *BG = World::getInstance().getBondGraph();

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

  ColorList = new enum GraphEdge::Shading[getAtomCount()];
  for (int i=0;i<getAtomCount();i++)
    ColorList[i] = (enum GraphEdge::Shading)0;
  if (flag) {
    // remove bonds that are beyond bonddistance
    Translationvector.Zero();
    // scan all bonds
    flag = false;
    for(molecule::iterator AtomRunner = begin(); (!flag) && (AtomRunner != end()); ++AtomRunner) {
      const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
      for(BondList::const_iterator BondRunner = ListOfBonds.begin();
          (!flag) && (BondRunner != 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;
          const range<double> MinMaxDistance(
              BG->getMinMaxDistance(Binder->leftatom, Binder->rightatom));
          if (!MinMaxDistance.isInRange(tmp)) {
            DoLog(2) && (Log() << Verbose(2) << "Correcting at bond " << *Binder << "." << endl);
            flag = true;
            break;
          }
        }
      }
    }
    //if (flag) {
    if (0) {
      // create translation vector from their periodically modified distance
      for (int i=NDIM;i--;) {
        tmp = Binder->leftatom->at(i) - Binder->rightatom->at(i);
        const range<double> MinMaxDistance(
            BG->getMinMaxDistance(Binder->leftatom, Binder->rightatom));
        if (fabs(tmp) > MinMaxDistance.last)  // check against Min is not useful for components
          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] = GraphEdge::white;
      AtomStack->push_front(Binder->leftatom);
      while (!AtomStack->empty()) {
        Walker = AtomStack->front();
        AtomStack->pop_front();
        //Log() << Verbose (3) << "Current Walker is: " << *Walker << "." << endl;
        ColorList[Walker->getNr()] = GraphEdge::black;    // mark as explored
        *Walker += Translationvector; // translate
        const BondList& ListOfBonds = Walker->getListOfBonds();
        for (BondList::const_iterator Runner = ListOfBonds.begin();
            Runner != ListOfBonds.end();
            ++Runner) {
          if ((*Runner) != Binder) {
            OtherWalker = (*Runner)->GetOtherAtom(Walker);
            if (ColorList[OtherWalker->getNr()] == GraphEdge::white) {
              AtomStack->push_front(OtherWalker); // push if yet unexplored
            }
          }
        }
      }
//      // re-add bond
//      if (OtherBinder == NULL) { // is the only bond?
//        //Do nothing
//      } else {
//        if (!LastBond) {
//          link(Binder, OtherBinder); // no more implemented bond::previous ...
//        } else {
//          link(OtherBinder, Binder); // no more implemented bond::previous ...
//        }
//      }
    } 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);

  return flag;
};