Changes in src/Fragmentation/Fragmentation.cpp [94d5ac6:3aa8a5]

File:

• src/Fragmentation/Fragmentation.cpp (modified) (19 diffs)

src/Fragmentation/Fragmentation.cpp

@@ -46 +46 @@
 #include "Element/periodentafel.hpp"
 #include "Fragmentation/AdaptivityMap.hpp"
+#include "Fragmentation/AtomMask.hpp"
 #include "Fragmentation/fragmentation_helpers.hpp"
 #include "Fragmentation/Graph.hpp"
+#include "Fragmentation/helpers.hpp"
 #include "Fragmentation/KeySet.hpp"
 #include "Fragmentation/PowerSetGenerator.hpp"
 #include "Fragmentation/UniqueFragments.hpp"
 #include "Graph/BondGraph.hpp"
-#include "Graph/CheckAgainstAdjacencyFile.hpp"
+#include "Graph/AdjacencyList.hpp"
+#include "Graph/ListOfLocalAtoms.hpp"
 #include "molecule.hpp"
-#include "MoleculeLeafClass.hpp"
-#include "MoleculeListClass.hpp"
-#include "Parser/FormatParserStorage.hpp"
 #include "World.hpp"
 
@@ -63 +63 @@
  *
  * \param _mol molecule for internal use (looking up atoms)
+ * \param _FileChecker instance contains adjacency parsed from elsewhere
  * \param _saturation whether to treat hydrogen special and saturate dangling bonds or not
  */
-Fragmentation::Fragmentation(molecule *_mol, const enum HydrogenSaturation _saturation) :
+Fragmentation::Fragmentation(molecule *_mol, AdjacencyList &_FileChecker, const enum HydrogenSaturation _saturation) :
   mol(_mol),
-  saturation(_saturation)
+  saturation(_saturation),
+  FileChecker(_FileChecker)
 {}
 
@@ -89 +91 @@
  * of vertex indices: Global always means the index in "this" molecule, whereas local refers to the molecule or
  * subgraph in the MoleculeListClass.
+ * \param atomids atomic ids (local to Fragmentation::mol) to fragment, used in AtomMask
  * \param Order up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
  * \param prefix prefix string for every fragment file name (may include path)
@@ -94 +97 @@
  * \return 1 - continue, 2 - stop (no fragmentation occured)
  */
-int Fragmentation::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;
+int Fragmentation::FragmentMolecule(const std::vector<atomId_t> &atomids, int Order, std::string prefix, DepthFirstSearchAnalysis &DFS)
+{
+  std::fstream File;
   bool CheckOrder = false;
-  Graph **FragmentList = NULL;
-  Graph TotalGraph;     // graph with all keysets however local numbers
   int TotalNumberOfKeySets = 0;
-  atom ***ListOfLocalAtoms = NULL;
-  bool *AtomMask = NULL;
-
-  LOG(0, endl);
+
+  LOG(0, std::endl);
   switch (saturation) {
     case DoSaturate:
@@ -123 +117 @@
 
   // ++++++++++++++++++++++++++++ INITIAL STUFF: Bond structure analysis, file parsing, ... ++++++++++++++++++++++++++++++++++++++++++
+  bool FragmentationToDo = true;
 
   // ===== 1. Check whether bond structure is same as stored in files ====
@@ -128 +123 @@
   // === compare it with adjacency file ===
   {
-    std::ifstream File;
-    std::string filename;
-    filename = prefix + ADJACENCYFILE;
-    File.open(filename.c_str(), ios::out);
-    LOG(1, "Looking at bond structure stored in adjacency file and comparing to present one ... ");
-
-    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();
+    const std::vector<atomId_t> globalids = getGlobalIdsFromLocalIds(*mol, atomids);
+    AdjacencyList WorldAdjacency(globalids);
+    FragmentationToDo = FragmentationToDo && (FileChecker > WorldAdjacency);
+  }
+
+  // ===== 2. create AtomMask that takes Saturation condition into account
+  AtomMask_t AtomMask(atomids);
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    // remove in hydrogen and we do saturate
+    if ((saturation == DoSaturate) && ((*iter)->getType()->getAtomicNumber() == 1)) // skip hydrogen
+      AtomMask.setFalse((*iter)->getNr());
+  }
 
   // ===== 3. if structure still valid, parse key set file and others =====
@@ -156 +141 @@
 
   // ===== 4. check globally whether there's something to do actually (first adaptivity check)
-  FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(prefix);
+  FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(atomids, 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(mol, &ParsedFragmentList, ListOfLocalAtoms, FragmentList, FragmentCounter, true);
-  delete[](ListOfLocalAtoms);
+  ListOfLocalAtoms_t ListOfLocalAtoms;
+  Graph FragmentList;
+  AssignKeySetsToFragment(ParsedFragmentList, ListOfLocalAtoms, FragmentList, true);
 
   // ===== 6b. prepare and go into the adaptive (Order<0), single-step (Order==0) or incremental (Order>0) cycle
-  KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
-  AtomMask = new bool[mol->getAtomCount()+1];
-  AtomMask[mol->getAtomCount()] = false;
+  KeyStack RootStack;
   FragmentationToDo = false;  // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards
-  while ((CheckOrder = CheckOrderAtSite(AtomMask, &ParsedFragmentList, Order, prefix))) {
+  bool LoopDoneAlready = false;
+  while ((CheckOrder = CheckOrderAtSite(AtomMask, ParsedFragmentList, Order, prefix, LoopDoneAlready))) {
     FragmentationToDo = FragmentationToDo || CheckOrder;
-    AtomMask[mol->getAtomCount()] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
+    LoopDoneAlready = 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), saturation);
-
-    // ===== 7. fill the bond fragment list =====
-    FragmentCounter = 0;
-    MolecularWalker = Subgraphs;
-    while (MolecularWalker->next != NULL) {
-      MolecularWalker = MolecularWalker->next;
-      LOG(1, "Fragmenting subgraph " << MolecularWalker << ".");
-      if (MolecularWalker->Leaf->hasBondStructure()) {
-        // call BOSSANOVA method
-        LOG(0, endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= ");
-        FragmentBOSSANOVA(MolecularWalker->Leaf, FragmentList[FragmentCounter], RootStack[FragmentCounter]);
-      } else {
-        ELOG(1, "Subgraph " << MolecularWalker << " has no atoms!");
-      }
-      FragmentCounter++;  // next fragment list
-    }
+    FillRootStackForSubgraphs(RootStack, AtomMask);
+
+    // call BOSSANOVA method
+    FragmentBOSSANOVA(mol, FragmentList, RootStack);
   }
   LOG(2, "CheckOrder is " << CheckOrder << ".");
-  delete[](RootStack);
-  delete[](AtomMask);
 
   // ==================================== 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);
-
-  LOG(0, FragmentCounter << " subgraph fragments have been removed.");
-
-  // ===== 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;
-    LOG(0, "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << ".");
-    BondFragments->insert(StoreFragmentFromKeySet(test, World::getInstance().getConfig()));
-    k++;
-  }
-  LOG(0, k << "/" << BondFragments->ListOfMolecules.size() << " fragments generated from the keysets.");
-
-  // ===== 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);
-
-    LOG(1, "Writing " << BondFragments->ListOfMolecules.size() << " possible bond fragmentation configs");
-    bool write_status = true;
-    for (std::vector<std::string>::const_iterator iter = typelist.begin();
-        iter != typelist.end();
-        ++iter) {
-      LOG(2, "INFO: Writing bond fragments for type " << (*iter) << ".");
-      write_status = write_status
-      && BondFragments->OutputConfigForListOfFragments(
-          prefix,
-          SortIndex,
-          FormatParserStorage::getInstance().getTypeFromName(*iter));
-    }
-    if (write_status)
-      LOG(1, "All configs written.");
-    else
-      LOG(1, "Some config writing failed.");
-
-    // store force index reference file
-    BondFragments->StoreForcesFile(prefix, SortIndex);
-
-    // store keysets file
-    TotalGraph.StoreKeySetFile(prefix);
-
-    {
-      // store Adjacency file
-      std::string filename = prefix + ADJACENCYFILE;
-      mol->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
-    LOG(1, "Freeing bond memory");
-    delete[](SortIndex);
-  } else {
-    LOG(1, "FragmentList is zero on return, splitting failed.");
-  }
-  // 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);
+  TranslateIndicesToGlobalIDs(FragmentList, TotalNumberOfKeySets, TotalGraph);
+
+  LOG(1, "STATUS: We have created " << TotalGraph.size() << " fragments.");
+
+
+  // store adaptive orders into file
+  StoreOrderAtSiteFile(prefix);
+
   LOG(0, "End of bond fragmentation.");
-
   return ((int)(!FragmentationToDo)+1);    // 1 - continue, 2 - stop (no fragmentation occured)
 };
@@ -317 +196 @@
  * \return pointer to Graph list
  */
-void Fragmentation::FragmentBOSSANOVA(molecule *mol, Graph *&FragmentList, KeyStack &RootStack)
-{
+void Fragmentation::FragmentBOSSANOVA(molecule *mol, Graph &FragmentList, KeyStack &RootStack)
+{
+  Info FunctionInfo(__func__);
   Graph ***FragmentLowerOrdersList = NULL;
   int NumLevels = 0;
@@ -330 +210 @@
   int RootNr = 0;
   UniqueFragments FragmentSearch;
-
-  LOG(0, "Begin of FragmentBOSSANOVA.");
 
   // FragmentLowerOrdersList is a 2D-array of pointer to MoleculeListClass objects, one dimension represents the ANOVA expansion of a single order (i.e. 5)
@@ -420 +298 @@
   FreeAllOrdersList(FragmentLowerOrdersList, RootStack, mol);
   delete[](NumMoleculesOfOrder);
-
-  LOG(0, "End of FragmentBOSSANOVA.");
 };
-
-/** 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 * Fragmentation::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
-{
-  Info info(__func__);
-  atom **SonList = new atom*[mol->getAtomCount()+1];
-  molecule *Leaf = World::getInstance().createMolecule();
-
-  for(int i=0;i<=mol->getAtomCount();i++)
-    SonList[i] = NULL;
-
-  StoreFragmentFromKeySet_Init(mol, Leaf, Leaflet, SonList);
-  // create the bonds between all: Make it an induced subgraph and add hydrogen
-//  LOG(2, "Creating bonds from father graph (i.e. induced subgraph creation).");
-  CreateInducedSubgraphOfFragment(mol, Leaf, SonList, IsAngstroem);
-
-  //Leaflet->Leaf->ScanForPeriodicCorrection(out);
-  delete[](SonList);
-  return Leaf;
-};
-
 
 /** Estimates by educated guessing (using upper limit) the expected number of fragments.
@@ -471 +319 @@
     c = (ListOfBonds.size() > c) ? ListOfBonds.size() : c;
   }
-  FragmentCount = mol->getNoNonHydrogen()*(1 << (c*order));
+  FragmentCount = (saturation == DoSaturate ? mol->getNoNonHydrogen() : mol->getAtomCount()) *(1 << (c*order));
   LOG(1, "Upper limit for this subgraph is " << FragmentCount << " for "
       << mol->getNoNonHydrogen() << " non-H atoms with maximum bond degree of " << c << ".");
@@ -479 +327 @@
 
 /** 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 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)
+ * \param LoopDoneAlready indicate whether we have done a fragmentation loop already
  * \return true - needs further fragmentation, false - does not need fragmentation
  */
-bool Fragmentation::CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, std::string path)
+bool Fragmentation::CheckOrderAtSite(AtomMask_t &AtomMask, const Graph &GlobalKeySetList, int Order, std::string path, bool LoopDoneAlready)
 {
   bool status = false;
 
-  // initialize mask list
-  for(int i=mol->getAtomCount();i--;)
-    AtomMask[i] = false;
-
   if (Order < 0) { // adaptive increase of BondOrder per site
-    if (AtomMask[mol->getAtomCount()] == true)  // break after one step
+    if (LoopDoneAlready)  // break after one step
       return false;
 
     // transmorph graph keyset list into indexed KeySetList
-    if (GlobalKeySetList == NULL) {
-      ELOG(1, "Given global key set list (graph) is NULL!");
-      return false;
-    }
-    AdaptivityMap * IndexKeySetList = GlobalKeySetList->GraphToAdaptivityMap();
+    AdaptivityMap * IndexKeySetList = GlobalKeySetList.GraphToAdaptivityMap();
 
     // parse the EnergyPerFragment file
@@ -512 +353 @@
     if (IndexKeySetList->IsAdaptiveCriteriaListEmpty()) {
       ELOG(2, "Unable to parse file, incrementing all.");
-      for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
-        if ((saturation == DontSaturate) || ((*iter)->getType()->getAtomicNumber() != 1)) // skip hydrogen
-        {
-          AtomMask[(*iter)->getNr()] = true;  // include all (non-hydrogen) atoms
-          status = true;
-        }
-      }
+      status = true;
     } else {
-      IndexKeySetList->MarkUpdateCandidates(AtomMask, Order, mol);
+      // mark as false all sites that are below threshold already
+      status = IndexKeySetList->MarkUpdateCandidates(AtomMask, Order, mol);
     }
 
@@ -526 +362 @@
   } else { // global increase of Bond Order
     for(molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
-      if ((saturation == DontSaturate) || ((*iter)->getType()->getAtomicNumber() != 1)) // skip hydrogen
-      {
-        AtomMask[(*iter)->getNr()] = true;  // include all (non-hydrogen) atoms
-        if ((Order != 0) && ((*iter)->AdaptiveOrder < Order)) // && ((*iter)->AdaptiveOrder < MinimumRingSize[(*iter)->getNr()]))
+      if (AtomMask.isTrue((*iter)->getNr())) { // skip masked out
+        // remove all that have reached desired order
+        if ((Order != 0) && ((*iter)->AdaptiveOrder >= Order)) // && ((*iter)->AdaptiveOrder < MinimumRingSize[(*iter)->getNr()]))
+          AtomMask.setFalse((*iter)->getNr());
+        else
           status = true;
       }
     }
-    if ((!Order) && (!AtomMask[mol->getAtomCount()]))  // single stepping, just check
+    if ((!Order) && (!LoopDoneAlready))  // single stepping, just check
       status = true;
 
@@ -576 +413 @@
 /** Parses pairs(Atom::Nr, Atom::AdaptiveOrder) from file and stores in molecule's Atom's.
  * Atoms not present in the file get "0".
+ * \param atomids atoms to fragment, used in AtomMask
  * \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 Fragmentation::ParseOrderAtSiteFromFile(std::string &path)
-{
-  unsigned char *OrderArray = new unsigned char[mol->getAtomCount()];
-  bool *MaxArray = new bool[mol->getAtomCount()];
+bool Fragmentation::ParseOrderAtSiteFromFile(const std::vector<atomId_t> &atomids, std::string &path)
+{
+  Info FunctionInfo(__func__);
+  typedef unsigned char order_t;
+  typedef std::map<atomId_t, order_t> OrderArray_t;
+  OrderArray_t OrderArray;
+  AtomMask_t MaxArray(atomids);
   bool status;
   int AtomNr, value;
@@ -589 +430 @@
   ifstream file;
 
-  for(int i=0;i<mol->getAtomCount();i++) {
-    OrderArray[i] = 0;
-    MaxArray[i] = false;
-  }
-
-  LOG(1, "Begin of ParseOrderAtSiteFromFile");
   line = path + ORDERATSITEFILE;
   file.open(line.c_str());
@@ -605 +440 @@
         OrderArray[AtomNr] = value;
         file >> value;
-        MaxArray[AtomNr] = value;
+        MaxArray.setValue(AtomNr, (bool)value);
         //LOG(2, "AtomNr " << AtomNr << " with order " << (int)OrderArray[AtomNr] << " and max order set to " << (int)MaxArray[AtomNr] << ".");
       }
@@ -614 +449 @@
     for(molecule::iterator iter=mol->begin();iter!=mol->end();++iter){
       (*iter)->AdaptiveOrder = OrderArray[(*iter)->getNr()];
-      (*iter)->MaxOrder = MaxArray[(*iter)->getNr()];
+      (*iter)->MaxOrder = MaxArray.isTrue((*iter)->getNr());
     }
     //SetAtomValueToIndexedArray( OrderArray, &atom::getNr(), &atom::AdaptiveOrder );
@@ -625 +460 @@
     status = false;
   }
-  delete[](OrderArray);
-  delete[](MaxArray);
-
-  LOG(1, "End of ParseOrderAtSiteFromFile");
+
   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 Fragmentation::CreateMappingLabelsToConfigSequence(int *&SortIndex)
-{
-  if (SortIndex != NULL) {
-    LOG(1, "SortIndex is " << SortIndex << " and not NULL as expected.");
+/** Fills the root stack for sites to be used as root in fragmentation depending on order or adaptivity criteria
+ * Again, as in \sa FillBondStructureFromReference steps recursively through each Leaf in this chain list of molecule's.
+ * \param &RootStack stack to be filled
+ * \param AtomMask defines true/false per global Atom::Nr to mask in/out each nuclear site
+ * \return true - stack is non-empty, fragmentation necessary, false - stack is empty, no more sites to update
+ */
+void Fragmentation::FillRootStackForSubgraphs(KeyStack &RootStack, const AtomMask_t &AtomMask)
+{
+  for(molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    const atom * const Father = (*iter)->GetTrueFather();
+    if (AtomMask.isTrue(Father->getNr())) // apply mask
+      if ((saturation == DontSaturate) || ((*iter)->getType()->getAtomicNumber() != 1)) // skip hydrogen
+        RootStack.push_front((*iter)->getNr());
+  }
+}
+
+/** The indices per keyset are compared to the respective father's Atom::Nr in each subgraph and thus put into \a **&FragmentList.
+ * \param *KeySetList list with all keysets
+ * \param ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
+ * \param **&FragmentList list to be allocated and returned
+ * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
+ * \retuen true - success, false - failure
+ */
+bool Fragmentation::AssignKeySetsToFragment(Graph &KeySetList, ListOfLocalAtoms_t &ListOfLocalAtoms, Graph &FragmentList, bool FreeList)
+{
+  Info FunctionInfo(__func__);
+  bool status = true;
+  size_t KeySetCounter = 0;
+
+  // fill ListOfLocalAtoms if NULL was given
+  if (!mol->FillListOfLocalAtoms(ListOfLocalAtoms, mol->getAtomCount())) {
+    LOG(1, "Filling of ListOfLocalAtoms failed.");
     return false;
   }
-  SortIndex = new int[mol->getAtomCount()+1];
-  for(int i=mol->getAtomCount()+1;i--;)
-    SortIndex[i] = -1;
-
-  int AtomNo = 0;
-  for(molecule::const_iterator iter=mol->begin();iter!=mol->end();++iter){
-    ASSERT(SortIndex[(*iter)->getNr()]==-1,"Same SortIndex set twice");
-    SortIndex[(*iter)->getNr()] = AtomNo++;
-  }
-
-  return true;
-};
-
-
-/** 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 Fragmentation::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 Fragmentation::CreateInducedSubgraphOfFragment(molecule *mol, molecule *Leaf, atom **SonList, bool IsAngstroem)
-{
-  bool LonelyFlag = false;
-  atom *OtherFather = NULL;
-  atom *FatherOfRunner = NULL;
-
-  // we increment the iter just before skipping the hydrogen
-  // as we use AddBond, we cannot have a const_iterator here
-  for (molecule::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);
-        if (SonList[OtherFather->getNr()] != NULL) {
-//          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
-//              << " is bound to " << *OtherFather << ", whose son is "
-//              << *SonList[OtherFather->getNr()] << ".");
-          if (OtherFather->getNr() > FatherOfRunner->getNr()) { // add bond (Nr check is for adding only one of both variants: ab, ba)
-            std::stringstream output;
-//            output << "ACCEPT: Adding Bond: "
-            output << Leaf->AddBond((*iter), SonList[OtherFather->getNr()], (*BondRunner)->BondDegree);
-//            LOG(3, output.str());
-            //NumBonds[(*iter)->getNr()]++;
-          } else {
-//            LOG(3, "REJECY: Not adding bond, labels in wrong order.");
-          }
-          LonelyFlag = false;
-        } else {
-//          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
-//              << " is bound to " << *OtherFather << ", who has no son in this fragment molecule.");
-          if (saturation == DoSaturate) {
-//          LOG(3, "ACCEPT: Adding Hydrogen to " << (*iter)->Name << " and a bond in between.");
-            if (!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem))
-              exit(1);
-          }
-          //NumBonds[(*iter)->getNr()] += Binder->BondDegree;
-        }
+
+  if (KeySetList.size() != 0) { // if there are some scanned keysets at all
+    // assign scanned keysets
+    KeySet TempSet;
+    for (Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++) { // key sets contain global numbers!
+      if (ListOfLocalAtoms[mol->FindAtom(*((*runner).first.begin()))->getNr()] != NULL) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
+        // translate keyset to local numbers
+        for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
+          TempSet.insert(ListOfLocalAtoms[mol->FindAtom(*sprinter)->getNr()]->getNr());
+        // insert into FragmentList
+        FragmentList.insert(GraphPair(TempSet, pair<int, double> (KeySetCounter++, (*runner).second.second)));
       }
-    } else {
-    ELOG(1, "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->getNr()] << "!");
-    }
-    if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) {
-      LOG(0, **iter << "has got bonds only to hydrogens!");
-    }
-    ++iter;
-    if (saturation == DoSaturate) {
-      while ((iter != Leaf->end()) && ((*iter)->getType()->getAtomicNumber() == 1)){ // skip added hydrogen
-        iter++;
-      }
-    }
-  }
-};
-
-/** Sets the desired output types of the fragment configurations.
- *
- * @param types vector of desired types.
- */
-void Fragmentation::setOutputTypes(const std::vector<std::string> &types)
-{
-  typelist = types;
+      TempSet.clear();
+    }
+  } else
+    LOG(1, "KeySetList is NULL or empty.");
+
+  if (FreeList) {
+    // free the index lookup list
+    ListOfLocalAtoms.clear();
+  }
+  return status;
 }
+
+/** Translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
+ * \param &FragmentList Graph with local numbers per fragment
+ * \param &TotalNumberOfKeySets global key set counter
+ * \param &TotalGraph Graph to be filled with global numbers
+ */
+void Fragmentation::TranslateIndicesToGlobalIDs(Graph &FragmentList, int &TotalNumberOfKeySets, Graph &TotalGraph)
+{
+  Info FunctionInfo(__func__);
+  for (Graph::iterator runner = FragmentList.begin(); runner != FragmentList.end(); runner++) {
+    KeySet TempSet;
+    for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
+      TempSet.insert((mol->FindAtom(*sprinter))->GetTrueFather()->getId());
+    TotalGraph.insert(GraphPair(TempSet, pair<int, double> (TotalNumberOfKeySets++, (*runner).second.second)));
+  }
+}
+