Index: src/analyzer.cpp
===================================================================
--- src/analyzer.cpp	(revision 260934983a76c25770b6a1f5706e2fc4eb1eecec)
+++ src/analyzer.cpp	(revision db942effae4d25e57224c90078ac288434a4a783)
@@ -58,8 +58,8 @@
 
   // ---------- Parse the TE Factors into an array -----------------
-  if (!Energy.ParseIndices()) return 1;
+  //if (!Energy.ParseIndices()) return 1;
   
   // ---------- Parse the Force indices into an array ---------------
-  if (!Force.ParseIndices(argv[1])) return 1;
+  //if (!Force.ParseIndices(argv[1])) return 1;
 
   // ---------- Parse the KeySets into an array ---------------
Index: src/atom.cpp
===================================================================
--- src/atom.cpp	(revision 260934983a76c25770b6a1f5706e2fc4eb1eecec)
+++ src/atom.cpp	(revision db942effae4d25e57224c90078ac288434a4a783)
@@ -26,4 +26,5 @@
   SeparationVertex = false;
   LowpointNr = -1;
+  AdaptiveOrder = 0;
 };
 
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision 260934983a76c25770b6a1f5706e2fc4eb1eecec)
+++ src/builder.cpp	(revision db942effae4d25e57224c90078ac288434a4a783)
@@ -559,8 +559,5 @@
 void FragmentAtoms(molecule *mol, config *configuration)
 {
-  enum BondOrderScheme Scheme = NoScheme;
-  enum CutCyclicBond CutCyclic;
-  char schema;
-  int Order1, Order2;
+  int Order1;
   clock_t start, end;
   
@@ -568,42 +565,7 @@
   cout << Verbose(0) << "What's the desired bond order: ";
   cin >> Order1;
-  cout << Verbose(0) << "What's the desired bond scheme [(B)ottomUp/(T)opDown/(A)NOVA/(C)ombined]: ";
-  cin >> schema;
-  CutCyclic = SaturateBond;
-  switch (schema) {
-    default:
-        Scheme = NoScheme;
-      break;
-    case 'B':
-        CutCyclic = KeepBond;
-    case 'b':
-        Scheme = BottomUp;
-        Order2 = 0;
-      break;
-    case 'T':
-        CutCyclic = KeepBond;
-    case 't':
-        Scheme = TopDown;
-        Order2 = Order1;
-        Order1 = 0;
-      break;
-    case 'A':
-        CutCyclic = KeepBond;
-    case 'a':
-        Scheme = ANOVA;
-        Order2 = 0;
-      break;
-    case 'C':
-        CutCyclic = KeepBond;
-    case 'c':
-        cout << Verbose(0) << "Combined first selects subgraphs by BottomUp of bond order " << Order1 << " then TopDown to fragment these." << endl;
-        cout << Verbose(0) << "What's the desired bond order for TopDown: ";
-        cin >> Order2;
-        Scheme = Combined;
-      break;
-  };
   if (mol->first->next != mol->last) {  // there are bonds
     start = clock();
-    mol->FragmentMolecule((ofstream *)&cout, Order1, Order2, Scheme, configuration, CutCyclic);
+    mol->FragmentMolecule((ofstream *)&cout, Order1, configuration);
     end = clock();
     cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
@@ -759,6 +721,4 @@
   int j, axis, count, faktor;
   int MinimumRingSize = -1;
-  enum BondOrderScheme Scheme = NoScheme;
-  enum CutCyclicBond CutCyclic;
   enum ConfigStatus config_present = absent;
   MoleculeLeafClass *Subgraphs = NULL;
@@ -781,9 +741,10 @@
           case '?':
             cout << "MoleCuilder suite" << endl << "==================" << endl << endl;
-            cout << "Usage: " << argv[0] << "[config file] [-{acepsthH?vfrp}] [further arguments]" << endl; 
+            cout << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl; 
             cout << "or simply " << argv[0] << " without arguments for interactive session." << endl;
             cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;
             cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
             cout << "\t-e <file>\tSets the element database to be parsed from this file (default: elements.db in same dir as " << argv[0] << ")." << endl;
+            cout << "\t-e <dist> <order>\tFragments the molecule in BOSSANOVA manner and stores config files in same dir as config." << endl;
             cout << "\t-h/-H/-?\tGive this help screen." << endl;
             cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
@@ -954,59 +915,19 @@
                 break;
               case 'f':
-                int i,j;
                 if (flag ==0) flag = 2;  // only set if not already by other command line switch
-                cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << " and type " << argv[argptr+2] << "." << endl;
-                if (argc >= argptr+3) {
+                cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;
+                if (argc >= argptr+2) {
                   cout << Verbose(0) << "Creating connection matrix..." << endl;
                   start = clock();
                   mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr++]));
                   cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
-                  cout << Verbose(0) << "What's the desired bond scheme [(B)ottomUp/(T)opDown/(A)NOVA/(C)ombined]: ";
-                  CutCyclic = SaturateBond;
-                  switch (argv[argptr+1][0]) {
-                    default:
-                        Scheme = NoScheme;
-                      break;
-                    case 'B':
-                      CutCyclic = KeepBond;
-                    case 'b':
-                        Scheme = BottomUp;
-                        i = atoi(argv[argptr]);
-                        j = 0;
-                      break;
-                    case 'T':
-                      CutCyclic = KeepBond;
-                    case 't':
-                        Scheme = TopDown;
-                        i = 0;
-                        j = atoi(argv[argptr]);
-                      break;
-                    case 'A':
-                      CutCyclic = KeepBond;
-                    case 'a':
-                        Scheme = ANOVA;
-                        i = atoi(argv[argptr]);
-                        j = 0;
-                      break;
-    //                case 'C':
-    //                  CutCyclic = KeepBond;
-    //                case 'c':
-    //                    if (argc > 5) {
-    //                      i = atoi(argv[4]);
-    //                      j = atoi(argv[6]);
-    //                      Scheme = Combined;
-    //                    }
-    //                    else
-    //                      cerr << "Missing second bond order for TopDown fragmentation in combined approach." << endl;
-    //                  break;
-                  };
                   if (mol->first->next != mol->last) {
-                    mol->FragmentMolecule((ofstream *)&cout, i, j, Scheme, &configuration, CutCyclic);
+                    mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);
                   }
                   end = clock();
                   cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
-                  argptr+=2;
+                  argptr+=1;
                 } else {
-                  cerr << "Not enough arguments for fragmentation: -f <max. bond distance> <bond order> <type>" << endl;
+                  cerr << "Not enough arguments for fragmentation: -f <max. bond distance> <bond order>" << endl;
                 }
                 break;
Index: src/defs.hpp
===================================================================
--- src/defs.hpp	(revision 260934983a76c25770b6a1f5706e2fc4eb1eecec)
+++ src/defs.hpp	(revision db942effae4d25e57224c90078ac288434a4a783)
@@ -18,5 +18,5 @@
 #define VERSIONSTRING "v1.0"
 
-enum BondOrderScheme { NoScheme, BottomUp, TopDown, ANOVA, Combined }; //!< Fragmentation scheme used in BOSS
+//enum BondOrderScheme { NoScheme, BottomUp, TopDown, ANOVA, Combined }; //!< Fragmentation scheme used in BOSS
 
 enum EdgeType { Undetermined, TreeEdge, BackEdge }; //!< edge type in a graph after Depth-First-Search analysis.
@@ -24,5 +24,5 @@
 enum Shading { white, lightgray, darkgray, black };  //!< color in Breadth-First-Search analysis
 
-enum CutCyclicBond { KeepBond,  SaturateBond }; //!< Saturation scheme either atom- or bondwise 
+//enum CutCyclicBond { KeepBond,  SaturateBond }; //!< Saturation scheme either atom- or bondwise 
 
 // Specifting whether a value in the parameter file must be specified or is optional
@@ -44,6 +44,5 @@
 #define KEYSETFILE "KeySets.dat"
 #define ADJACENCYFILE "Adjacency.dat"
-#define TEFACTORSFILE "TE-Factors.dat"
-#define FORCESFILE "Forces-Factors.dat"
+#define ORDERATSITEFILE "OrderAtSite.dat" 
 #define FRAGMENTPREFIX "BondFragment"
 #define STANDARDCONFIG "unknown.conf"
@@ -54,4 +53,5 @@
 #define STANDARDHBONDANGLEDB "Hbondangle.db"
 
+#define UPDATECOUNT 10  //!< update ten sites per BOSSANOVA interval
 
 #endif /*DEFS_HPP_*/
Index: src/joiner.cpp
===================================================================
--- src/joiner.cpp	(revision 260934983a76c25770b6a1f5706e2fc4eb1eecec)
+++ src/joiner.cpp	(revision db942effae4d25e57224c90078ac288434a4a783)
@@ -46,8 +46,8 @@
 
   // ---------- Parse the TE Factors into an array -----------------
-  if (!Energy.ParseIndices()) return 1;
+  //if (!Energy.ParseIndices()) return 1;
   
   // ---------- Parse the Force indices into an array ---------------
-  if (!Force.ParseIndices(argv[1])) return 1;
+  //if (!Force.ParseIndices(argv[1])) return 1;
   //if (!ParseForceIndices(argv[1], data.ForcesIndices, data.AtomCounter, data.FragmentCounter)) return 1;
 
Index: src/moleculelist.cpp
===================================================================
--- src/moleculelist.cpp	(revision 260934983a76c25770b6a1f5706e2fc4eb1eecec)
+++ src/moleculelist.cpp	(revision db942effae4d25e57224c90078ac288434a4a783)
@@ -21,10 +21,7 @@
 MoleculeListClass::MoleculeListClass(int NumMolecules, int NumAtoms)
 {
-  TEList = (double*) Malloc(sizeof(double)*NumMolecules, "MoleculeListClass:MoleculeListClass: *TEList");
   ListOfMolecules = (molecule **) Malloc(sizeof(molecule *)*NumMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
-  for (int i=0;i<NumMolecules;i++) {
-    TEList[i] = 0.;
+  for (int i=0;i<NumMolecules;i++)
     ListOfMolecules[i] = NULL;
-  }
   NumberOfMolecules = NumMolecules;
   NumberOfTopAtoms = NumAtoms;
@@ -36,5 +33,5 @@
 MoleculeListClass::~MoleculeListClass()
 {
-  cout << Verbose(3) << this << ": Freeing ListOfMolcules and TEList." << endl;
+  cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
   for (int i=0;i<NumberOfMolecules;i++) {
     if (ListOfMolecules[i] != NULL) { // if NULL don't free
@@ -46,8 +43,11 @@
   cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
   Free((void **)&ListOfMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
-  cout << Verbose(4) << "Freeing TEList." << endl;
-  Free((void **)&TEList, "MoleculeListClass:MoleculeListClass: *TEList");
-};
-
+};
+
+/** Compare whether two molecules are equal.
+ * \param *a molecule one
+ * \param *n molecule two
+ * \return lexical value (-1, 0, +1)
+ */
 int MolCompare(const void *a, const void *b)
 {
@@ -122,165 +122,4 @@
 };
 
-/** Returns an allocated index list which fragment should remain and which are doubly appearing.
- * Again, we use linked cell however not with \a threshold which in most cases is too small and would
- * generate too many cells, but a too be specified \a celldistance which should e.g. be chosen in
- * the magnitude of the typical bond distance. In order to avoid so far the need of a vector list,
- * we simply use the molecule structure and the vector sitting in the atom structure.
- * \param *out output stream for debugging
- * \param **MapList empty but allocated(MoleculeListClass::NumberOfMolecules), containing on return mapping of for which equal molecule which atoms corresponds to which one (needed for ForceFactors)
- * \param threshold threshold value for molecule::IsEqualToWithinThreshold()
- * \param cell_size 6 entries specifying the unit cell vectors (matrix is symmetric)
- * \param celldistance needed for linked-cell technique to find possible equals
- * \return allocated index list with \a Num entries.
- */
-int * MoleculeListClass::GetMappingToUniqueFragments(ofstream *out, double threshold, double *cell_size, double celldistance)
-{
-  int j, UniqueIndex, HighestNumber;
-  //int NumberCells, divisor[NDIM], n[NDIM], N[NDIM], index, Index;
-  int Counter;
-  //molecule **CellList;
-  int *MapList = NULL;
-  size_t *SortMap = NULL;
-  //atom *Walker = NULL, *OtherWalker = NULL;
-  //molecule *LinkedCellContainer = new molecule(NULL);  // container for all the center of gravities stored as atoms
-  
-  /// Allocated MapList with range NumberOfMolecules. 
-  MapList = (int *) Malloc(sizeof(int)*NumberOfMolecules, "MoleculeListClass::GetMappingToUniqueFragments: *MapList");
-  SortMap = (size_t *) Malloc(sizeof(size_t)*NumberOfMolecules, "MoleculeListClass::GetMappingToUniqueFragments: *MapList");
-  for (j=0;j<NumberOfMolecules;j++) {
-    if (ListOfMolecules[j] ==  NULL)
-      cerr << "ERROR: Molecule " << j << " is NULL!" << endl;
-    //else 
-      //cerr << "ERROR: Molecule " << j << " is " << ListOfMolecules[j] << " with count " << ListOfMolecules[j]->AtomCount << "." << endl;
-    MapList[j] = -1;
-    SortMap[j] = 0;
-  }
-  *out << Verbose(1) << "Get mapping to unique fragments with " << NumberOfMolecules << " fragments total." << endl;
-  
-  // sort the list with heapsort according to sort function
-  gsl_heapsort_index(SortMap, ListOfMolecules, NumberOfMolecules, sizeof(molecule *), MolCompare);
-  // check next neighbours and remap
-  Counter = 0;
-  MapList[SortMap[0]] = Counter++;
-  for(int i=1;i<NumberOfMolecules;i++) {
-    if (MolCompare(&ListOfMolecules[SortMap[i]], &ListOfMolecules[SortMap[i-1]]) == 0)
-      MapList[SortMap[i]] = MapList[SortMap[i-1]];
-    else
-      MapList[SortMap[i]] = Counter++; 
-  } 
-  Free((void **)&SortMap, "MoleculeListClass::GetMappingToUniqueFragments: *SortMap");
-
-  *out << Verbose(2) << "Map: ";
-  for(int i=0;i<NumberOfMolecules;i++)
-    *out << MapList[i] << " ";
-  *out << endl;
-
-  // bring MapList indices into an ascending order
-  HighestNumber = Counter;
-  *out << Verbose(3) << "HighestNumber: " << HighestNumber << "." << endl;
-  SortMap = (size_t *) Malloc(sizeof(size_t)*NumberOfMolecules, "MoleculeListClass::GetMappingToUniqueFragments: *SortMap");
-  for(int i=0;i<NumberOfMolecules;i++)
-    SortMap[i] = HighestNumber;   // is the first number that will not appear anymore in list   
-  UniqueIndex = 0;
-  for(int i=0;i<NumberOfMolecules;i++)
-    if (MapList[i] != -1) {
-      if ((unsigned int)SortMap[MapList[i]] == (unsigned int)HighestNumber)
-        SortMap[MapList[i]] = UniqueIndex++;
-      MapList[i] = SortMap[MapList[i]];
-    } else {
-      MapList[i] = -1;
-    }
-  *out << Verbose(2) << "Ascending Map: ";
-  for(int i=0;i<NumberOfMolecules;i++)
-    *out << MapList[i] << " ";
-  *out << endl;
-  Free((void **)&SortMap, "MoleculeListClass::GetMappingToUniqueFragments: *SortMap");
-
-  /// Return the constructed list.
-  return MapList;
-};
-
-
-/** Takes a list of molecules and returns the list with doubles removed.
- * Checks by using molecule::IsEqualToWithinThreshold() whether fragments within the list
- * are actually doubles. If not, the pointer are copied to a new list, if so,
- * they are dropped. The new list is then copied back to the reallocated \a **FragmentList
- * and the new number of elements therein returned.
- * \param *out output stream for debugging
- * \param *Map mapping of which index goes on which (from molecule::GetMappingToUniqueFragments())
- * \return number of molecules in the new realloacted **FragmentList
- */
-int MoleculeListClass::ReduceFragmentToUniqueOnes(ofstream *out, int *Map)
-{
-  *out << Verbose(2) << "Begin of ReduceFragmentToUniqueOnes" << endl;
-  /// Allocate temporary lists of size \a Num.
-  molecule **NewMolList = (molecule **) Malloc(sizeof(molecule *)*NumberOfMolecules, "MoleculeListClass::ReduceFragmentToUniqueOnes: NewList");
-  double *NewTEList = (double *) Malloc(sizeof(double)*NumberOfMolecules, "MoleculeListClass::ReduceFragmentToUniqueOnes: *NewTEList");
-  for(int i=0;i<NumberOfMolecules;i++) {
-    NewTEList[i] = 0;
-    NewMolList[i] = NULL;
-  }
-  int ReducedNum = 0, j;
-  *out << Verbose(3) << "Reducing TE and Forces lists." << endl;
-  for(int i=0;i<NumberOfMolecules;i++) {                       /// Count new reduced number of molecules ...
-    //*out << Verbose(4) << i << "th TE Value " << " goes to " << Map[i] << ": " << TEList[i] << "." << endl;
-    NewTEList[Map[i]] += TEList[i];
-    if (Map[i] == i) {                           /// ... by checking if Map[i] != i, then Fragment[i] is double of Fragment[Map[i]].
-      *out << Verbose(4) << "Copying molecule " << i << " as it is unique so far ... " << ListOfMolecules[i]<< endl;
-      NewMolList[ReducedNum++] = ListOfMolecules[i];   /// ..., whilst copying each value to temporary list, ...
-    } else {  // else free molecule cause it appears doubly
-      *out << Verbose(4) << "Removing molecule " << i << "." << endl;
-      delete(ListOfMolecules[i]);
-      ListOfMolecules[i] = NULL;
-    }
-  }
-  /// Reallocate \a **FragmentList ...
-  *out << Verbose(3) << "Reallocating to Unique Fragments:" << endl;
-  if ((ReducedNum != 0) && (TEList != 0)) {
-    // factor list
-    *out << Verbose(4) << "Reallocating TEList [" << TEList << "] to " << ReducedNum << endl;
-    TEList = (double *) ReAlloc(TEList, sizeof(double)*ReducedNum, "MoleculeListClass::ReduceFragmentToUniqueOnes: *TEList");
-
-    *out << Verbose(4) << "Copying values to new list: " << endl;
-    j = 0;  // is ReducedNum index
-    for(int i=0;i<NumberOfMolecules;i++) {/// copy value from temporary to return list
-      if (Map[i] == i) {
-      //if (fabs(NewTEList[i]) > MYEPSILON) { // this is a good check also, if TEList[i] = 0 then fragment cancels and may be dropped!
-        TEList[j] = NewTEList[i];
-        *out << Verbose(5) << "TE [i" << i << "<->j" << j << "]:" << NewTEList[i] << "->" << TEList[j];
-        j++;
-      }
-      *out << endl;
-    }
-    if (j != ReducedNum)
-      *out << "Panic:  j " << j << " != ReducedNum " << ReducedNum << "." << endl;
-      
-
-    // molecule list
-    *out << Verbose(4) << "Reallocating ListOfMolecules ... " << endl;
-    ListOfMolecules = (molecule **) ReAlloc(ListOfMolecules, sizeof(molecule *)*ReducedNum, "MoleculeListClass::ReduceFragmentToUniqueOnes: ListOfMolecules");
-    NumberOfMolecules = ReducedNum;
-    /// ... and copy the reduced number back
-    *out << Verbose(5) << "Transfering to ListOfMolecules ... ";
-    for(int i=0;i<ReducedNum;i++) {
-      ListOfMolecules[i] = NewMolList[i];
-      *out << NewMolList[i] << " -> " << ListOfMolecules[i] << "\t";
-    }
-    *out << endl;
-
-
-  } else {
-    Free((void **)&ListOfMolecules, "MoleculeListClass::ReduceFragmentToUniqueOnes: ListOfMolecules");
-    NumberOfMolecules = ReducedNum;
-    *out << Verbose(3) << "ReducedNum is " << ReducedNum << ", Number of Molecules is " << NumberOfMolecules << "." << endl;
-  }
-  // free all the lists again
-  *out << Verbose(3) << "Freeing temporary lists." << endl;
-  Free((void **)&NewTEList, "MoleculeListClass::ReduceFragmentToUniqueOnes: *NewTEList");      /// free temporary list again
-  Free((void **)&NewMolList, "MoleculeListClass::ReduceFragmentToUniqueOnes: *NewMolList");
-  *out << Verbose(2) << "End of ReduceFragmentToUniqueOnes" << endl;
-  return(ReducedNum);
-};
-
 /** Simple output of the pointers in ListOfMolecules.
  * \param *out output stream
@@ -295,11 +134,10 @@
 
 /** Writes a config file for each molecule in the given \a **FragmentList.
- * Note that molecules with a TEList factor of 0 are not stored!
- * \param *prefix prefix for the file name
+ * \param *out output stream for debugging
  * \param *configuration standard configuration to attach atoms in fragment molecule to.
  * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
  * \return true - success (each file was written), false - something went wrong.
  */
-bool MoleculeListClass::OutputConfigForListOfFragments(char *prefix, config *configuration, int *SortIndex)
+bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex)
 {
   ofstream outputFragment;
@@ -310,352 +148,102 @@
   atom *Walker = NULL;
   vector BoxDimension;
-  char TEFilename[MAXSTRINGSIZE];
   char *FragmentNumber;
   int FragmentCounter = 0;
   ofstream output;
-  element *runner = NULL;
   
   // store the fragments as config and as xyz
   for(int i=0;i<NumberOfMolecules;i++) {
-    //if (TEList[i] != 0) {
-      strcpy(PathBackup, configuration->configpath);
-
-      // scan all atoms in the current molecule for their fathers and write each vertex index to forces file
-
-      // correct periodic
-      ListOfMolecules[i]->ScanForPeriodicCorrection((ofstream *)&cout);
-
-      // output xyz file
-      FragmentNumber = FixedDigitNumber(NumberOfMolecules, FragmentCounter++);
-      sprintf(FragmentName, "%s/%s%s.conf.xyz", PathBackup, prefix, FragmentNumber);
-      outputFragment.open(FragmentName, ios::out);
-      cout << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << " as XYZ ...";
-      if (intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment))
-        cout << " done." << endl;
-      else
-        cout << " failed." << endl;
-      result = result && intermediateResult;
-      outputFragment.close();
-      outputFragment.clear();
-  
-      cout << Verbose(2) << "Contained atoms: ";
-      Walker = ListOfMolecules[i]->start;
-      while (Walker->next != ListOfMolecules[i]->end) {
-        Walker = Walker->next;
-        cout << Walker->Name << " ";
-      }
-      cout << endl;
-      // prepare output of config file
-      sprintf(FragmentName, "%s/%s%s.conf", PathBackup, prefix, FragmentNumber);
-      outputFragment.open(FragmentName, ios::out);
-      strcpy(PathBackup, configuration->configpath);
-      sprintf(FragmentName, "%s/%s%s/", PathBackup, prefix, FragmentNumber);
-      
-      // center on edge
-      ListOfMolecules[i]->CenterEdge((ofstream *)&cout, &BoxDimension);
-      ListOfMolecules[i]->SetBoxDimension(&BoxDimension);  // update Box of atoms by boundary
-      int j = -1;
-      for (int k=0;k<3;k++) {
-        j += k+1;
-        BoxDimension.x[k] = 5.;
-        ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
-      }
-      ListOfMolecules[i]->Translate(&BoxDimension);
-
-      // also calculate necessary orbitals
-      ListOfMolecules[i]->CountElements();  // this is a bugfix, atoms should should actually be added correctly to this fragment 
-      ListOfMolecules[i]->CalculateOrbitals(*configuration);
-      
-      // change path in config
-      configuration->SetDefaultPath(FragmentName);
-      cout << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << " as config ...";
-      if (intermediateResult = configuration->Save(&outputFragment, ListOfMolecules[i]->elemente, ListOfMolecules[i]))
-        cout << " done." << endl;
-      else
-        cout << " failed." << endl;
-      // restore old config
-      configuration->SetDefaultPath(PathBackup);
-      
-      result = result && intermediateResult;
-      outputFragment.close();
-      outputFragment.clear();
-      Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
-    }
-
-  strcpy(PathBackup, configuration->configpath);
-  // open file for the total energy factor
-  sprintf(TEFilename, "%s/%s%s", PathBackup, prefix, TEFACTORSFILE);
-  output.open(TEFilename, ios::out);
-  // output TEList (later to file AtomicTEList.dat)
-  cout << Verbose(2) << "Saving " << prefix << " energy factors ...";
-  //cout << Verbose(1) << "Final ATEList: ";
-  //less output << prefix << "TE\t";
-  for(int i=0;i<NumberOfMolecules;i++) {
-    //cout << TEList[i] << "\t";
-    //if (TEList[i] != 0)
-      output << TEList[i] << "\t";
-  }
-  //cout << endl;
-  output << endl;
+    // save default path as it is changed for each fragment
+    strcpy(PathBackup, configuration->GetDefaultPath());
+
+    // correct periodic
+    ListOfMolecules[i]->ScanForPeriodicCorrection(out);
+
+    // output xyz file
+    FragmentNumber = FixedDigitNumber(NumberOfMolecules, FragmentCounter++);
+    sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    outputFragment.open(FragmentName, ios::out);
+    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << " as XYZ ...";
+    if (intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment))
+      *out << " done." << endl;
+    else
+      *out << " failed." << endl;
+    result = result && intermediateResult;
+    outputFragment.close();
+    outputFragment.clear();
+
+    *out << Verbose(2) << "Contained atoms: ";
+    Walker = ListOfMolecules[i]->start;
+    while (Walker->next != ListOfMolecules[i]->end) {
+      Walker = Walker->next;
+      *out << Walker->Name << " ";
+    }
+    *out << endl;
+    // prepare output of config file
+    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    outputFragment.open(FragmentName, ios::out);
+    strcpy(PathBackup, configuration->configpath);
+    sprintf(FragmentName, "%s/%s%s/", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    
+    // center on edge
+    ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
+    ListOfMolecules[i]->SetBoxDimension(&BoxDimension);  // update Box of atoms by boundary
+    int j = -1;
+    for (int k=0;k<3;k++) {
+      j += k+1;
+      BoxDimension.x[k] = 5.;
+      ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
+    }
+    ListOfMolecules[i]->Translate(&BoxDimension);
+
+    // also calculate necessary orbitals
+    ListOfMolecules[i]->CountElements();  // this is a bugfix, atoms should should actually be added correctly to this fragment 
+    ListOfMolecules[i]->CalculateOrbitals(*configuration);
+    
+    // change path in config
+    configuration->SetDefaultPath(FragmentName);
+    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << " as config ...";
+    if (intermediateResult = configuration->Save(&outputFragment, ListOfMolecules[i]->elemente, ListOfMolecules[i]))
+      *out << " done." << endl;
+    else
+      *out << " failed." << endl;
+    // restore old config
+    configuration->SetDefaultPath(PathBackup);
+    
+    result = result && intermediateResult;
+    outputFragment.close();
+    outputFragment.clear();
+    Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
+  }
+
+  // open KeySet file
+  sprintf(FragmentName, "%s/%s%s", configuration->configpath,  FRAGMENTPREFIX , KEYSETFILE);
+  output.open(FragmentName, ios::out);
+  *out << Verbose(2) << "Saving " << FRAGMENTPREFIX << " key sets of each subgraph ...";
+  for(int j=0;j<NumberOfMolecules;j++) {
+    Walker = ListOfMolecules[j]->start;
+    while(Walker->next != ListOfMolecules[j]->end) {
+      Walker = Walker->next;        
+#ifdef ADDHYDROGEN
+      if ((Walker->GetTrueFather() != NULL) && (Walker->type->Z != 1)) // if there is a real father, prints its index
+#else
+      if ((Walker->GetTrueFather() != NULL)) // if there is a real father, prints its index
+#endif 
+        output <<  Walker->GetTrueFather()->nr << "\t";
+#ifdef ADDHYDROGEN
+      else  // otherwise a -1 to indicate an added saturation hydrogen
+        output << "-1\t";
+#endif 
+    }
+    output << endl;
+  }
   output.close();
-  cout << " done." << endl;
-  
-  // open file for the force factors
-  sprintf(TEFilename, "%s/%s%s", PathBackup, prefix, FORCESFILE);
-  output.open(TEFilename, ios::out);
-  //cout << Verbose(1) << "Final AtomicForcesList: ";
-  cout << Verbose(2) << "Saving " << prefix << " force factors ...";
-  //output << prefix << "Forces" << endl;
-  for(int j=0;j<NumberOfMolecules;j++) {
-    //if (TEList[j] != 0) {
-    runner = ListOfMolecules[j]->elemente->start;
-    while (runner->next != ListOfMolecules[j]->elemente->end) { // go through every element
-      runner = runner->next;
-      if (ListOfMolecules[j]->ElementsInMolecule[runner->Z]) { // if this element got atoms
-        Walker = ListOfMolecules[j]->start;
-        while (Walker->next != ListOfMolecules[j]->end) { // go through every atom of this element
-          Walker = Walker->next;
-          if (Walker->type->Z == runner->Z) {
-            if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a real father, prints its index
-              //cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", its number " << Walker->GetTrueFather()->nr << " and index " << SortIndex[Walker->GetTrueFather()->nr] << "." << endl;  
-              output <<  SortIndex[Walker->GetTrueFather()->nr] << "\t";
-            } else  // otherwise a -1 to indicate an added saturation hydrogen
-              output << "-1\t";
-          }
-        }
-      }
-    }
-    //cout << endl << endl;
-    output << endl;
-  }
-  output.close();
-  cout << " done." << endl;
-
-  // open KeySet file
-  sprintf(TEFilename, "%s/%s%s", PathBackup, prefix, "KeySets.dat");
-  output.open(TEFilename, ios::out);
-  cout << Verbose(2) << "Saving " << prefix << " key sets of each subgraph ...";
-  for(int j=0;j<NumberOfMolecules;j++) {
-    //if (TEList[j] != 0) {
-      Walker = ListOfMolecules[j]->start;
-      while(Walker->next != ListOfMolecules[j]->end) {
-        Walker = Walker->next;        
-#ifdef ADDHYDROGEN
-        if ((Walker->GetTrueFather() != NULL) && (Walker->type->Z != 1)) // if there is a real father, prints its index
-#else
-        if ((Walker->GetTrueFather() != NULL)) // if there is a real father, prints its index
-#endif 
-          output <<  Walker->GetTrueFather()->nr << "\t";
-#ifdef ADDHYDROGEN
-        else  // otherwise a -1 to indicate an added saturation hydrogen
-          output << "-1\t";
-#endif 
-      }
-      //cout << endl << endl;
-      output << endl;
-    }
-  output.close();
-  cout << " done." << endl;
+  *out << " done." << endl;
   
   // printing final number
-  cout << "Final number of fragments: " << FragmentCounter << "." << endl;
+  *out << "Final number of fragments: " << FragmentCounter << "." << endl;
       
   return result;
 };
-
-/** Reduce the list of molecules to unique pieces.
- * molecule::IsEqualToWithinThreshold() is used by GetMappingToUniqueFragments() to provide a mapping for
- * ReduceFragmentToUniqueOnes(). 
- * \param *out output stream for debugging
- * \param cell_size 6 entries specifying the unit cell vectors (matrix is symmetric) for getting the map
- * \param celldistance needed for linked-cell technique when getting the map
- */
-void MoleculeListClass::ReduceToUniqueList(ofstream *out, double *cell_size, double celldistance)
-{
-  int *ReduceMapping = NULL;
-  
-  *out << Verbose(0) << "Begin of ReduceToUniqueList." << endl; 
-  // again, check whether there are equal fragments by ReduceToUniqueOnes
-  *out << Verbose(1) << "Get reduce mapping." << endl;
-  ReduceMapping = GetMappingToUniqueFragments(out, 0.05, cell_size, celldistance);
-  *out << Verbose(1) << "MapList: ";
-  for(int i=0;i<NumberOfMolecules;i++)
-    *out << ReduceMapping[i] << " ";
-  *out << endl << endl;
-  *out << Verbose(1) << "Apply reduce mapping." << endl;
-  ReduceFragmentToUniqueOnes(out, ReduceMapping);
-  Free((void **)&ReduceMapping, "MoleculeListClass::FragmentTopDown: *ReduceMapping");  // free map after reduce
-  *out << Verbose(1) << "Number of Fragments after Reducing is " << NumberOfMolecules << "." << endl;
-  *out << Verbose(0) << "End of ReduceToUniqueList." << endl; 
-};
-
-/** Fragments a molecule and resulting pieces recursively until the number of bonds
- * Here the idea is similar to the FragmentBottomUp() approach, only that we split up the molecule
- * bond per bond into left and right fragments and this recursively also with each yielded
- * fragment until there is no fragment with continuous number of bonds greater than \a BondDegree.
- * are below the given \a BondDegree.
- * \param *out output stream for debugging
- * \param BondDegree maximum number of continuous bonds in a fragment
- * \param bonddistance typical bond distance
- * \param *configuration configuration for writing config files for each fragment
- * \param CutCyclic cut cyclic bonds (saturate with hydrogen) or add 
- * \return pointer to MoleculeListClass with all the fragments or NULL if something went wrong.
- * \todo so far we volontarily mix up the above BondDegree definition and molecule::NoNonBonds
- * 
- * \bug  FragmentTopDown does not work right now due to NULL being given as cell_size to ReduceToUniqueOnes()
- */
-MoleculeListClass * MoleculeListClass::FragmentTopDown(ofstream *out, int BondDegree, double bonddistance, config *configuration, enum CutCyclicBond CutCyclic)
-{
-  int Num, No, Cyclic, NoBonds;
-  MoleculeListClass *ReturnList = NULL;
-  MoleculeListClass **FragmentsList = NULL;
-  molecule *CurrentMolecule = NULL;
-  bond *Binder = NULL;
- 
-  *out << Verbose(0) << "Begin of FragmentTopDown:" << this << "." << endl; 
-  Num = 0;
-  FragmentsList = (MoleculeListClass **) Malloc(sizeof(MoleculeListClass *)*NumberOfMolecules, "MoleculeListClass::FragmentTopDown: **FragmentsList");
-  // fragment each molecule into its own MoleculeListClass
-  for(int i=0;i<NumberOfMolecules;i++) {
-    CurrentMolecule = ListOfMolecules[i];
-    CurrentMolecule->CountAtoms(out);
-    Cyclic = CurrentMolecule->CountCyclicBonds(out);
-    NoBonds = 
-#ifdef ADDHYDROGEN
-CurrentMolecule->NoNonBonds
-#else 
-CurrentMolecule->BondCount
-#endif 
-                - ((CutCyclic == KeepBond) ? Cyclic : 0);
-    // check if NoNonBonds < BondDegree, then don't split any further: Here, correct for greatest continuous bond length!
-    *out << Verbose(0) << "Checking on Number of true Bonds " << NoBonds << " (i.e. no -H, if chosen no cyclic) greater than " << BondDegree << "." << endl;
-    if (NoBonds > BondDegree) {
-      //cerr << Verbose(1) << "TopDown Level of Molecule " << CurrentMolecule << ": " << (CurrentMolecule->NoNonBonds - BondDegree) << "." << endl; 
-      *out << Verbose(1) << "Breaking up molecule " << i << " in fragment list." << endl;
-      CurrentMolecule->CreateListOfBondsPerAtom(out);
-      // Get atomic fragments for the list
-      *out << Verbose(1) << "Getting Atomic fragments for molecule " << i << "." << endl;
-      MoleculeListClass *AtomicFragments = CurrentMolecule->GetAtomicFragments(out, NumberOfTopAtoms, configuration->GetIsAngstroem(), TEList[i]/(1. - NoBonds), CutCyclic);
-      
-      // build the atomic part of the list of molecules
-      *out << Verbose(1) << "Putting atomic fragment into list of molecules." << endl;
-      
-      // cutting cyclic bonds yields only a single fragment, not two as non-cyclic!
-      No = (CutCyclic == KeepBond) ? ((
-#ifdef ADDHYDROGEN
-CurrentMolecule->NoNonBonds - Cyclic)*2) : CurrentMolecule->NoNonBonds*2 - Cyclic; 
-#else 
-CurrentMolecule->BondCount - Cyclic)*2) : CurrentMolecule->BondCount*2 - Cyclic; 
-#endif 
-      
-      if (AtomicFragments != NULL) {
-        ReturnList = new MoleculeListClass(No+AtomicFragments->NumberOfMolecules, NumberOfTopAtoms);
-        for (int j=0;j<AtomicFragments->NumberOfMolecules;j++) { // shift all Atomic Fragments into this one here
-          ReturnList->ListOfMolecules[j+No] = AtomicFragments->ListOfMolecules[j];
-          AtomicFragments->ListOfMolecules[j] = NULL;
-          ReturnList->ListOfMolecules[j+No]->NoNonBonds = 0;
-          ReturnList->ListOfMolecules[j+No]->NoNonHydrogen = 1;
-          ReturnList->TEList[j+No] = AtomicFragments->TEList[j];
-        }
-        *out << Verbose(3) << "Freeing Atomic memory" << endl;
-        delete(AtomicFragments);
-        AtomicFragments = NULL;
-      } else {
-        *out << Verbose(2) << "AtomicFragments is NULL, filling list with whole molecule only." << endl; 
-        ReturnList = new MoleculeListClass(No, NumberOfTopAtoms);
-      }
-
-      // build the side pieces part of the list of molecules
-      *out << Verbose(1) << "Building side piece fragments and putting into list of molecules." << endl;
-      No = 0;
-      Binder = CurrentMolecule->first;
-      while(Binder->next != CurrentMolecule->last) {
-        Binder = Binder->next;
-#ifdef ADDHYDROGEN
-        if (Binder->HydrogenBond == 0)
-#endif        
-          if ((CutCyclic == SaturateBond) || (!Binder->Cyclic)) {
-            // split each bond into left and right side piece
-            if (Binder->Cyclic) {
-              molecule *dummy = NULL; // we lazily hand FragmentMoleculeByBond() a dummy as second fragment and ... 
-              *out << Verbose(1) << "Breaking up " << Binder->nr << "th and " << No << "th non-hydrogen, CYCLIC bond " << *Binder << " of molecule " << i << " [" << CurrentMolecule << "] into a single pieces." << endl; 
-              CurrentMolecule->FragmentMoleculeByBond(out, Binder, &(ReturnList->ListOfMolecules[No]), &(dummy), configuration->GetIsAngstroem(), CutCyclic);
-              delete(dummy);  // ... free the dummy which is due to the cyclic bond the exacy same fragment
-              *out << Verbose(1) << "Single fragment is " << ReturnList->ListOfMolecules[No] << "." << endl;
-      
-              // write down the necessary TE-summation in order to regain TE of whole molecule
-              *out << Verbose(2) << "Creating TEList for Fragment " << i << " of Bond " << No << "." << endl;
-#ifdef ADDHYDROGEN
-              ReturnList->TEList[No] = -TEList[i]/(1. -CurrentMolecule->NoNonBonds);
-#else              
-              ReturnList->TEList[No] = -TEList[i]/(1. -CurrentMolecule->BondCount);
-#endif
-            } else {
-              *out << Verbose(1) << "Breaking up " << Binder->nr << "th and " << No << "th non-hydrogen, non-cyclic bond " << *Binder << " of molecule " << i << " [" << CurrentMolecule << "] into left and right side pieces." << endl; 
-              CurrentMolecule->FragmentMoleculeByBond(out, Binder, &(ReturnList->ListOfMolecules[No]), &(ReturnList->ListOfMolecules[No+1]), configuration->GetIsAngstroem(), CutCyclic);
-              *out << Verbose(1) << "Left Fragment is " << ReturnList->ListOfMolecules[No] << ", right Fragment is " << ReturnList->ListOfMolecules[No+1] << "." << endl;
-      
-              // write down the necessary TE-summation in order to regain TE of whole molecule
-              *out << Verbose(2) << "Creating TEList for Fragment " << i << " of Bond " << No << "." << endl;
-              ReturnList->TEList[No] = -TEList[i]/(1. - NoBonds);
-              ReturnList->TEList[No+1] = -TEList[i]/(1. - NoBonds);
-            }
-            No += ((Binder->Cyclic) ? 1 : 2);
-          }
-      }
-
-      // Reduce this list
-      ReturnList->ReduceToUniqueList(out, NULL, bonddistance);
-
-      // recurse and receive List
-      *out << Verbose(1) << "Calling TopDown " << ReturnList<< " with filled list: ";
-      ReturnList->Output(out);
-      FragmentsList[i] = ReturnList->FragmentTopDown(out, BondDegree, bonddistance, configuration, CutCyclic);
-      *out << Verbose(1) << "Returning from TopDown " << ReturnList<< " with filled list: ";
-      ReturnList->Output(out);
-      
-      // remove list after we're done
-      *out << Verbose(2) << "Removing caller list." << endl;
-      delete(ReturnList);
-      ReturnList = NULL;
-      
-    } else {  // create a list with only a single molecule inside, transfer everything from "this" list to return list
-      *out << Verbose(1) << "Not breaking up molecule " << i << " in fragment list." << endl;
-      FragmentsList[i] = new MoleculeListClass(1, NumberOfTopAtoms);
-      FragmentsList[i]->ListOfMolecules[0] = ListOfMolecules[i];
-      ListOfMolecules[i] = NULL;
-      FragmentsList[i]->TEList[0] = TEList[i]; 
-    }
-    Num += FragmentsList[i]->NumberOfMolecules;
-  }
-  
-  // now, we have a list of MoleculeListClasses: combine all fragments list into one single list
-  *out << Verbose(2) << "Combining to a list of " << Num << " molecules and Memory cleanup of old list of lists." << endl; 
-  ReturnList = new MoleculeListClass(Num, NumberOfTopAtoms);
-  No = 0;
-  for(int i=0;i<NumberOfMolecules;i++) {
-    for(int j=0;j<FragmentsList[i]->NumberOfMolecules;j++) {
-      // transfer molecule
-      ReturnList->ListOfMolecules[No] = FragmentsList[i]->ListOfMolecules[j];
-      FragmentsList[i]->ListOfMolecules[j] = NULL;
-      // transfer TE factor
-      ReturnList->TEList[No] = FragmentsList[i]->TEList[j];
-      No++;
-    }
-    delete(FragmentsList[i]);
-    FragmentsList[i] = NULL;
-  }
-  Free((void **)&FragmentsList, "MoleculeListClass::FragmentTopDown: **FragmentsList");
-  
-  // Reduce the list to unique fragments
-  ReturnList->ReduceToUniqueList(out, NULL, bonddistance);
-  
-  // return pointer to list
-  *out << Verbose(0) << "Return with filled list: ";
-  ReturnList->Output(out);
-  
-  *out << Verbose(0) << "End of FragmentTopDown:" << this << "." << endl; 
-  return (ReturnList);
-}; 
 
 /******************************************* Class MoleculeLeafClass ************************************************/
Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision 260934983a76c25770b6a1f5706e2fc4eb1eecec)
+++ src/molecules.cpp	(revision db942effae4d25e57224c90078ac288434a4a783)
@@ -1898,5 +1898,5 @@
     return false;
   }
-  cout << Verbose(1) << "Parsing the KeySet file ... ";
+  cout << Verbose(1) << "Parsing the KeySet file ... " << endl;
   // open file and read
   sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, KEYSETFILE);
@@ -1942,5 +1942,5 @@
  * \return true - file written successfully, false - writing failed
  */
-bool molecule::StoreAdjacencyToFile(ofstream *out, char *path, int bondorder)
+bool molecule::StoreAdjacencyToFile(ofstream *out, char *path)
 {
   ofstream AdjacencyFile;
@@ -1953,5 +1953,4 @@
   cout << Verbose(1) << "Saving adjacency list ... ";
   if (AdjacencyFile != NULL) {
-    AdjacencyFile << "Order\t" << bondorder << endl;
     Walker = start;
     while(Walker->next != end) {
@@ -1977,8 +1976,7 @@
  * \param *path path to file
  * \param **ListOfAtoms allocated (molecule::AtomCount) and filled lookup table for ids (Atom::nr) to *Atom
- * \param bondorder check whether files matches desired bond order
  * \return true - structure is equal, false - not equivalence
  */ 
-bool molecule::CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms, int bondorder)
+bool molecule::CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms)
 {
   char *filename = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule - filename");
@@ -1995,53 +1993,46 @@
     int CurrentBondsOfAtom;
 
-    // determine order from file
-    File.getline(filename, MAXSTRINGSIZE);    
-    if (bondorder != atoi(&filename[5])) {
-      *out << "Bond order desired is " << bondorder << " and does not match one in file " << filename[6] << "." << endl;
-      status = false;
-    } else {
-      // Parse the file line by line and count the bonds
-      while (!File.eof()) {
-        File.getline(filename, MAXSTRINGSIZE);
-        stringstream line;
-        line.str(filename);
-        int AtomNr = -1;
-        line >> AtomNr;
-        CurrentBondsOfAtom = -1; // we count one too far due to line end
-        // parse into structure
-        if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
-          while (!line.eof())
-            line >> CurrentBonds[ ++CurrentBondsOfAtom ];
-          // compare against present bonds
-          //cout << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
-          if (CurrentBondsOfAtom == NumberOfBondsPerAtom[AtomNr]) {
-            for(int i=0;i<NumberOfBondsPerAtom[AtomNr];i++) {
-              int id = ListOfBondsPerAtom[AtomNr][i]->GetOtherAtom(ListOfAtoms[AtomNr])->nr;
-              int j = 0;
-              for (;(j<CurrentBondsOfAtom) && (CurrentBonds[j++] != id);); // check against all parsed bonds
-              if (CurrentBonds[j-1] != id) { // no match ? Then mark in ListOfAtoms
-                ListOfAtoms[AtomNr] = NULL;
-                NonMatchNumber++;
-                status = false;
-                //out << "[" << id << "]\t";
-              } else {
-                //out << id << "\t";
-              }
+    // Parse the file line by line and count the bonds
+    while (!File.eof()) {
+      File.getline(filename, MAXSTRINGSIZE);
+      stringstream line;
+      line.str(filename);
+      int AtomNr = -1;
+      line >> AtomNr;
+      CurrentBondsOfAtom = -1; // we count one too far due to line end
+      // parse into structure
+      if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
+        while (!line.eof())
+          line >> CurrentBonds[ ++CurrentBondsOfAtom ];
+        // compare against present bonds
+        //cout << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
+        if (CurrentBondsOfAtom == NumberOfBondsPerAtom[AtomNr]) {
+          for(int i=0;i<NumberOfBondsPerAtom[AtomNr];i++) {
+            int id = ListOfBondsPerAtom[AtomNr][i]->GetOtherAtom(ListOfAtoms[AtomNr])->nr;
+            int j = 0;
+            for (;(j<CurrentBondsOfAtom) && (CurrentBonds[j++] != id);); // check against all parsed bonds
+            if (CurrentBonds[j-1] != id) { // no match ? Then mark in ListOfAtoms
+              ListOfAtoms[AtomNr] = NULL;
+              NonMatchNumber++;
+              status = false;
+              //out << "[" << id << "]\t";
+            } else {
+              //out << id << "\t";
             }
-            //out << endl;
-          } else {
-            *out << "Number of bonds for Atom " << *ListOfAtoms[AtomNr] << " does not match, parsed " << CurrentBondsOfAtom << " against " << NumberOfBondsPerAtom[AtomNr] << "." << endl;
-            status = false;
           }
+          //out << endl;
+        } else {
+          *out << "Number of bonds for Atom " << *ListOfAtoms[AtomNr] << " does not match, parsed " << CurrentBondsOfAtom << " against " << NumberOfBondsPerAtom[AtomNr] << "." << endl;
+          status = false;
         }
       }
-      File.close();
-      File.clear();
-      if (status) { // if equal we parse the KeySetFile
-        *out << " done: Equal." << endl;
-        status = true;
-      } else
-        *out << " done: Not equal by " << NonMatchNumber << " atoms." << endl;
-    }
+    }
+    File.close();
+    File.clear();
+    if (status) { // if equal we parse the KeySetFile
+      *out << " done: Equal." << endl;
+      status = true;
+    } else
+      *out << " done: Not equal by " << NonMatchNumber << " atoms." << endl;
     Free((void **)&CurrentBonds, "molecule::CheckAdjacencyFileAgainstMolecule - **CurrentBonds");
   } else {
@@ -2057,11 +2048,8 @@
  * Writes for each fragment a config file.
  * \param *out output stream for debugging
- * \param BottomUpOrder up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
- * \param TopDownOrder up to how many neighbouring bonds a fragment contains in BondOrderScheme::TopDown scheme
- * \param Scheme which BondOrderScheme to use for the fragmentation
+ * \param Order up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
  * \param *configuration configuration for writing config files for each fragment
- * \param CutCyclic whether to add cut cyclic bond or to saturate
- */
-void molecule::FragmentMolecule(ofstream *out, int BottomUpOrder, int TopDownOrder, enum BondOrderScheme Scheme, config *configuration, enum CutCyclicBond CutCyclic)
+ */
+void molecule::FragmentMolecule(ofstream *out, int Order, config *configuration)
 {
 	MoleculeListClass **BondFragments = NULL;
@@ -2072,16 +2060,16 @@
   int AtomNo;
   int MinimumRingSize;
-  int TotalFragmentCounter = 0;
-  int FragmentCounter = 0;
+  int TotalFragmentCounter;
+  int FragmentCounter;
   MoleculeLeafClass *MolecularWalker = NULL;
   MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
   fstream File;
-  bool status = true;
-
-  cout << endl;
+  bool FragmentationToDo = true;
+
+  *out << endl;
 #ifdef ADDHYDROGEN
-  cout << Verbose(0) << "I will treat hydrogen special and saturate dangling bonds with it." << endl;
+  *out << Verbose(0) << "I will treat hydrogen special and saturate dangling bonds with it." << endl;
 #else
-  cout << Verbose(0) << "Hydrogen is treated just like the rest of the lot." << endl;
+  *out << Verbose(0) << "Hydrogen is treated just like the rest of the lot." << endl;
 #endif
 
@@ -2101,35 +2089,45 @@
   if (Walker->next != end) {  // everything went alright
     *out << " range of nuclear ids exceeded [0, AtomCount)." << endl;
-    status = false;
-  }
-  status = status && CheckAdjacencyFileAgainstMolecule(out, configuration->configpath, ListOfAtoms, BottomUpOrder); 
-  if (status) {  // NULL entries in ListOfAtoms contain NonMatches
-    status = status && ParseKeySetFile(out, configuration->configpath, ListOfAtoms, FragmentList, configuration->GetIsAngstroem());
-  }
+    FragmentationToDo = false;
+  }
+  FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(out, configuration->configpath, ListOfAtoms); 
+  if (FragmentationToDo)  // NULL entries in ListOfAtoms contain NonMatches
+    FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ListOfAtoms, FragmentList, configuration->GetIsAngstroem());
+  if (FragmentationToDo)  // parse the adaptive order per atom/site/vertex
+    FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(out, configuration->configpath);
   Free((void **)&ListOfAtoms, "molecule::FragmentMolecule - **ListOfAtoms");
 
   // =================================== Begin of FRAGMENTATION ===============================
-  if (!status) {
-    // === store Adjacency file ===
-    StoreAdjacencyToFile(out, configuration->configpath, BottomUpOrder);
-
+  if (FragmentationToDo) { // if we parsed Adjacancy, check whether OrderAtSite is above Order everywhere
+    Walker = start;
+    while (Walker->next != end) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron 
+      Walker = Walker->next;
+#ifdef ADDHYDROGEN
+      if (Walker->type->Z != 1) // skip hydrogen
+#endif
+        if (Walker->AdaptiveOrder < Order)
+          FragmentationToDo = false;
+    }
+  } 
+  
+  if (!FragmentationToDo) { // if we have still do something, FragmentationToDo will be false
     // === first perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs ===
-    Subgraphs = DepthFirstSearchAnalysis((ofstream *)&cout, false, MinimumRingSize);
+    Subgraphs = DepthFirstSearchAnalysis((ofstream *)&*out, false, MinimumRingSize);
     MolecularWalker = Subgraphs;
     // fill the bond structure of the individually stored subgraphs
     while (MolecularWalker->next != NULL) {
       MolecularWalker = MolecularWalker->next;
-      cout << Verbose(1) << "Creating adjacency list for subgraph " << MolecularWalker << "." << endl;
-      MolecularWalker->Leaf->CreateAdjacencyList((ofstream *)&cout, BondDistance);
-      MolecularWalker->Leaf->CreateListOfBondsPerAtom((ofstream *)&cout);
+      *out << Verbose(1) << "Creating adjacency list for subgraph " << MolecularWalker << "." << endl;
+      MolecularWalker->Leaf->CreateAdjacencyList(out, BondDistance);
+      MolecularWalker->Leaf->CreateListOfBondsPerAtom(out);
     }
     
     // === fragment all subgraphs ===
-    if ((MinimumRingSize != -1) && ((BottomUpOrder >= MinimumRingSize) || (TopDownOrder >= MinimumRingSize))) {
-      cout << Verbose(0) << "Bond order greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
+    if ((MinimumRingSize != -1) && (Order >= MinimumRingSize)) {
+      *out << Verbose(0) << "Bond order greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
     } else {
       FragmentCounter = 0;
       MolecularWalker = Subgraphs;
-      // count subgraphs
+      // count subgraphs and allocate fragments
       while (MolecularWalker->next != NULL) {
         MolecularWalker = MolecularWalker->next;
@@ -2137,10 +2135,35 @@
       }
       BondFragments = (MoleculeListClass **) Malloc(sizeof(MoleculeListClass *)*FragmentCounter, "molecule::FragmentMolecule - **BondFragments");
+      
+      // create RootStack for each Subgraph 
+      // NOTE: (keep this extern of following while loop, as lateron we may here look for which site to add to which subgraph)
+      KeyStack RootStack[FragmentCounter];
+
+      FragmentCounter = 0;
+      MolecularWalker = Subgraphs;
+      // count subgraphs and allocate fragments
+      while (MolecularWalker->next != NULL) {
+        MolecularWalker = MolecularWalker->next;
+        RootStack[FragmentCounter].clear();  
+        // find first root candidates  
+        Walker = MolecularWalker->Leaf->start;
+        while (Walker->next != MolecularWalker->Leaf->end) { // go through all (non-hydrogen) atoms
+          Walker = Walker->next;
+      #ifdef ADDHYDROGEN
+          if (Walker->type->Z != 1) // skip hydrogen
+      #endif
+            if (Walker->GetTrueFather()->AdaptiveOrder < Order) // only if Order is still greater
+              RootStack[FragmentCounter].push_front(Walker->nr);
+        }
+        FragmentCounter++;
+      }
+      
       // fill the bond fragment list
       FragmentCounter = 0;
+      TotalFragmentCounter = 0;
       MolecularWalker = Subgraphs;
       while (MolecularWalker->next != NULL) {
         MolecularWalker = MolecularWalker->next;
-        cout << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
+        *out << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
         if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
             // output ListOfBondsPerAtom for debugging
@@ -2165,38 +2188,22 @@
           *out << Verbose(0) << "Begin of bond fragmentation." << endl;
           BondFragments[FragmentCounter] = NULL;
-          if (Scheme == ANOVA) {
-            Graph *FragmentList = MolecularWalker->Leaf->FragmentBOSSANOVA(out,BottomUpOrder,configuration);
-
-            // allocate memory for the pointer array and transmorph graphs into full molecular fragments
-            int TotalNumberOfMolecules = 0;
-            for(Graph::iterator runner = FragmentList->begin(); runner != FragmentList->end(); runner++)
-              TotalNumberOfMolecules++;
-            BondFragments[FragmentCounter] = new MoleculeListClass(TotalNumberOfMolecules, AtomCount);
-            int k=0;
-            for(Graph::iterator runner = FragmentList->begin(); runner != FragmentList->end(); runner++) {
-              KeySet test = (*runner).first;
-              BondFragments[FragmentCounter]->ListOfMolecules[k] = StoreFragmentFromKeySet(out, test, configuration);
-              BondFragments[FragmentCounter]->TEList[k] = ((*runner).second).second;
-              k++;
-            }
-          }  
-          if ((Scheme == BottomUp) || (Scheme == Combined)) { // get overlapping subgraphs
-            BondFragments[FragmentCounter] = FragmentList = MolecularWalker->Leaf->FragmentBottomUp(out, BottomUpOrder, configuration, CutCyclic);
+          // call BOSSANOVA method
+          Graph *FragmentList = MolecularWalker->Leaf->FragmentBOSSANOVA(out, RootStack[FragmentCounter]);
+
+          // allocate memory for the pointer array and transmorph graphs into full molecular fragments
+          int TotalNumberOfMolecules = 0;
+          for(Graph::iterator runner = FragmentList->begin(); runner != FragmentList->end(); runner++)
+            TotalNumberOfMolecules++;
+          BondFragments[FragmentCounter] = new MoleculeListClass(TotalNumberOfMolecules, AtomCount);
+          int k=0;
+          for(Graph::iterator runner = FragmentList->begin(); runner != FragmentList->end(); runner++) {
+            KeySet test = (*runner).first;
+            *out << "Fragment No." << (*runner).second.first << " was created " << (int)(*runner).second.second << " time(s)." << endl;
+            BondFragments[FragmentCounter]->ListOfMolecules[k] = MolecularWalker->Leaf->StoreFragmentFromKeySet(out, test, configuration);
+            k++;
           }
-          if (Scheme == TopDown) { // initialise top level with whole molecule
-            *out << Verbose(2) << "Initial memory allocating and initialising for whole molecule." << endl;
-            FragmentList = new MoleculeListClass(1, MolecularWalker->Leaf->AtomCount);
-            FragmentList->ListOfMolecules[0] = MolecularWalker->Leaf->CopyMolecule();
-            FragmentList->TEList[0] = 1.;
-          }
-          if ((Scheme == Combined) || (Scheme == TopDown)) {
-            *out << Verbose(1) << "Calling TopDown." << endl;
-            BondFragments[FragmentCounter] = FragmentList->FragmentTopDown(out, TopDownOrder, BondDistance, configuration, CutCyclic);
-            // remove this molecule from list again and free wrapper list
-            delete(FragmentList);
-            FragmentList = NULL;
-          }
+          *out << k << "/" << BondFragments[FragmentCounter]->NumberOfMolecules << " fragments generated from the keysets." << endl;
         } else {
-          cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
+          *out << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
         }
         TotalFragmentCounter += BondFragments[FragmentCounter]->NumberOfMolecules;
@@ -2212,5 +2219,5 @@
         FragmentList->ListOfMolecules[TotalFragmentCounter] =  BondFragments[i]->ListOfMolecules[j];
         BondFragments[i]->ListOfMolecules[j] = NULL;
-        FragmentList->TEList[TotalFragmentCounter++] = BondFragments[i]->TEList[j];
+        TotalFragmentCounter++;
       }
       delete(BondFragments[i]);
@@ -2218,5 +2225,5 @@
     Free((void **)&BondFragments, "molecule::FragmentMolecule - **BondFragments");
   } else
-    cout << Verbose(0) << "Using fragments reconstructed from the KeySetFile." << endl;
+    *out << Verbose(0) << "Nothing to do, using only fragments reconstructed from the KeySetFile." << endl;
   // ==================================== End of FRAGMENTATION ================================
   
@@ -2241,9 +2248,15 @@
     }
     *out << Verbose(1) << "Writing " << FragmentList->NumberOfMolecules << " possible bond fragmentation configs" << endl;
-    if (FragmentList->OutputConfigForListOfFragments("BondFragment", configuration, SortIndex))
+    if (FragmentList->OutputConfigForListOfFragments(out, configuration, SortIndex))
       *out << Verbose(1) << "All configs written." << endl;
     else
       *out << Verbose(1) << "Some configs' writing failed." << endl;
     Free((void **)&SortIndex, "molecule::FragmentMolecule: *SortIndex"); 
+    
+    // === store Adjacency file ===
+    StoreAdjacencyToFile(out, configuration->configpath);
+
+    // Store adaptive orders into file
+    StoreOrderAtSiteFile(out, configuration->configpath);
     
     // restore orbital and Stop values
@@ -2256,4 +2269,5 @@
   } else
     *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
+  
   // free subgraph memory again
   if (Subgraphs != NULL) {
@@ -2266,4 +2280,79 @@
 
   *out << Verbose(0) << "End of bond fragmentation." << endl;
+};
+
+/** Stores pairs (Atom::nr, Atom::AdaptiveOrder) into file.
+ * Atoms not present in the file get "-1".
+ * \param *out output stream for debugging
+ * \param *path path to file ORDERATSITEFILE
+ * \return true - file writable, false - not writable
+ */
+bool molecule::StoreOrderAtSiteFile(ofstream *out, char *path)
+{
+  stringstream line;
+  ofstream file;
+
+  line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
+  file.open(line.str().c_str());
+  *out << Verbose(0) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl;
+  if (file != NULL) {
+    atom *Walker = start;
+    while (Walker->next != end) {
+      Walker = Walker->next;
+      file << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << endl;
+      *out << Verbose(2) << "Storing: " << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << "." << endl;
+    }
+    file.close();
+    return true;
+  } else {
+    return false;
+  }
+};
+
+/** Parses pairs(Atom::nr, Atom::AdaptiveOrder) from file and stores in molecule's Atom's.
+ * Atoms not present in the file get "0".
+ * \param *out output stream for debugging
+ * \param *path path to file ORDERATSITEFILEe
+ * \return true - file found and scanned, false - file not found
+ * \sa ParseKeySetFile() and CheckAdjacencyFileAgainstMolecule() as this is meant to be used in conjunction with the two
+ */
+bool molecule::ParseOrderAtSiteFromFile(ofstream *out, char *path)
+{
+  unsigned char *OrderArray = (unsigned char *) Malloc(sizeof(unsigned char)*AtomCount, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
+  bool status;
+  int AtomNr;
+  stringstream line;
+  ifstream file;
+  int Order;
+
+  *out << Verbose(1) << "Begin of ParseOrderAtSiteFromFile" << endl;
+  for(int i=0;i<AtomCount;i++)
+    OrderArray[i] = 0;
+  line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
+  file.open(line.str().c_str());
+  if (file != NULL) {
+    for (int i=0;i<AtomCount;i++) // initialise with 0
+      OrderArray[i] = 0;
+    while (!file.eof()) { // parse from file
+      file >> AtomNr;
+      file >> Order;
+      OrderArray[AtomNr] = (unsigned char) Order;
+      //*out << Verbose(2) << "AtomNr " << AtomNr << " with order " << (int)OrderArray[AtomNr] << "." << endl;
+    }
+    atom *Walker = start;
+    while (Walker->next != end) { // fill into atom classes
+      Walker = Walker->next;
+      Walker->AdaptiveOrder = OrderArray[Walker->nr];
+      *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << "." << endl; 
+    }
+    file.close();
+    status = true;
+  } else {
+    status = false;
+  }
+  Free((void **)&OrderArray, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
+  
+  *out << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl;
+  return status;
 };
 
@@ -2333,263 +2422,4 @@
 };
 
-/** Splits up a molecule into atomic, non-hydrogen, hydrogen-saturated fragments.
- * \param *out output stream for debugging
- * \param NumberOfTopAtoms number to initialise second parameter of MoleculeListClass with
- * \param IsAngstroem whether atomic coordination is in Angstroem (true) or atomic length/bohrradous (false)
- * \param factor additional factor TE and forces factors are multiplied with
- * \param CutCyclic whether to add cut cyclic bond or to saturate
- * \return MoleculelistClass of pointer to each atomic fragment.
- */
-MoleculeListClass * molecule::GetAtomicFragments(ofstream *out, int NumberOfTopAtoms, bool IsAngstroem, double factor, enum CutCyclicBond CutCyclic)
-{
-  atom *TopAtom = NULL, *BottomAtom = NULL; // Top = this, Bottom = AtomicFragment->ListOfMolecules[No]
-  atom *Walker = NULL;
-  MoleculeListClass *AtomicFragments = NULL;
-  atom **AtomList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::GetAtomicFragments: **AtomList");
-  for (int i=0;i<AtomCount;i++)
-    AtomList[i] = NULL;
-  bond **BondList = (bond **) Malloc(sizeof(bond *)*BondCount, "molecule::GetAtomicFragments: **AtomList");
-  for (int i=0;i<BondCount;i++)
-    BondList[i] = NULL;
-  int No = 0, Cyclic;
-
-  *out << Verbose(0) << "Begin of GetAtomicFragments." << endl;
-
-  *out << Verbose(1) << "Atoms in Molecule: ";
-  Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    *out << Walker << "\t";
-  } 
-  *out << endl;
-#ifdef ADDHYDROGEN
-  if (NoNonHydrogen != 0) {
-    AtomicFragments = new MoleculeListClass(NoNonHydrogen, NumberOfTopAtoms);
-  } else {
-    *out << Verbose(1) << "NoNonHydrogen is " << NoNonHydrogen << ", can't allocated MoleculeListClass." << endl;
-#else
-  if (AtomCount != 0) {
-    AtomicFragments = new MoleculeListClass(AtomCount, NumberOfTopAtoms);
-  } else {
-    *out << Verbose(1) << "AtomCount is " << AtomCount << ", can't allocated MoleculeListClass." << endl;
-#endif
-    return (AtomicFragments);
-  }
-  
-  TopAtom = start;
-  while (TopAtom->next != end) {
-    TopAtom = TopAtom->next;
-  //for(int i=0;i<AtomCount;i++) {
-#ifdef ADDHYDROGEN
-    if (TopAtom->type->Z != 1) {   // regard only non-hydrogen
-#endif
-      //TopAtom = AtomsInMolecule[i];
-      *out << Verbose(1) << "Current non-Hydrogen Atom: " << TopAtom->Name << endl;
-
-      // go through all bonds to check if cyclic
-      Cyclic = 0;
-      for(int i=0;i<NumberOfBondsPerAtom[TopAtom->nr];i++)
-        Cyclic += (ListOfBondsPerAtom[TopAtom->nr][i]->Cyclic) ? 1 : 0;
-
-#ifdef ADDHYDROGEN
-      if (No > NoNonHydrogen) {
-#else
-      if (No > AtomCount) {
-#endif
-        *out << Verbose(1) << "Access on created AtomicFragmentsListOfMolecules[" << No << "] beyond NumberOfMolecules " << AtomicFragments->NumberOfMolecules << "." << endl;
-        break;
-      }
-      if (AtomList[TopAtom->nr] == NULL) {
-        // create new molecule
-        AtomicFragments->ListOfMolecules[No] = new molecule(elemente); // allocate memory
-        AtomicFragments->TEList[No] = factor;
-        AtomicFragments->ListOfMolecules[No]->BondDistance = BondDistance;
-        
-        // add central atom
-        BottomAtom = AtomicFragments->ListOfMolecules[No]->AddCopyAtom(TopAtom); // add this central atom to molecule
-        AtomList[TopAtom->nr] = BottomAtom; // mark down in list
-
-        // create fragment          
-        *out << Verbose(1) << "New fragment around atom: " << TopAtom->Name << endl;
-        BreadthFirstSearchAdd(out,AtomicFragments->ListOfMolecules[No], AtomList, BondList, TopAtom, NULL, 0, IsAngstroem, (Cyclic == 0) ? SaturateBond : CutCyclic);
-        AtomicFragments->ListOfMolecules[No]->CountAtoms(out);
-        // actually we now have to reset both arrays to NULL again, but BFS is overkill anyway for getting the atomic fragments
-        // thus we do it in O(1) and avoid the O(n) which would make this routine O(N^2)!  
-        AtomList[TopAtom->nr] = NULL;  // remove this fragment's central atom again from the list
-        
-        *out << Verbose(1) << "Atoms in Fragment " << TopAtom->nr << ": ";
-        Walker = AtomicFragments->ListOfMolecules[No]->start;
-        while (Walker->next != AtomicFragments->ListOfMolecules[No]->end) {
-          Walker = Walker->next; 
-        //for(int k=0;k<AtomicFragments->ListOfMolecules[No]->AtomCount;k++)
-          *out << Walker << "(" << Walker->father << ")\t";
-        }
-        *out << endl;
-        No++;
-      }
-#ifdef ADDHYDROGEN
-    } else
-      *out << Verbose(1) << "Current Hydrogen Atom: " << TopAtom->Name << endl;
-#endif
-  }
-
-  // output of full list before reduction
-  if (AtomicFragments != NULL) {
-    *out << "AtomicFragments: ";
-    AtomicFragments->Output(out);
-    *out << endl;
-  } else 
-    *out << Verbose(1) << "AtomicFragments is zero on return, splitting failed." << endl;
-
-  // Reducing the atomic fragments
-  if (AtomicFragments != NULL) {
-    // check whether there are equal fragments by GetMappingToUniqueFragments
-    AtomicFragments->ReduceToUniqueList(out, &cell_size[0], BondDistance);
-  } else
-    *out << Verbose(1) << "AtomicFragments is zero, reducing failed." << endl;
-  Free((void **)&BondList, "molecule::GetAtomicFragments: **BondList");
-  Free((void **)&AtomList, "molecule::GetAtomicFragments: **AtomList");
-  *out << Verbose(0) << "End of GetAtomicFragments." << endl;
-  return (AtomicFragments);
-};
-
-/** Splits up the bond in a molecule into a left and a right fragment.
- * \param *out output stream for debugging
- * \param *Bond bond to broken up into getting allocated ...
- * \param *LeftFragment ... left fragment molecule ... (ptr to the memory cell wherein the adress for the Fragment is to be stored)
- * \param *RightFragment ... and right fragment molecule to be returned.(ptr to the memory cell wherein the adress for the Fragment is to be stored)
- * \param IsAngstroem whether atomic coordination is in Angstroem (true) or atomic length/bohrradous (false)
- * \param CutCyclic whether to add cut cyclic bond or not
- * \sa FragmentTopDown()
- */
-void molecule::FragmentMoleculeByBond(ofstream *out, bond *Bond, molecule **LeftFragment, molecule **RightFragment, bool IsAngstroem, enum CutCyclicBond CutCyclic)
-{
-  *out << Verbose(0) << "Begin of FragmentMoleculeByBond." << endl; 
-#ifdef ADDHYDROGEN
-  if ((Bond->leftatom->type->Z != 1) && (Bond->rightatom->type->Z != 1)) { // if both bond partners aren't hydrogen ...
-#endif    
-    *out << Verbose(1) << "Current Non-Hydrogen Bond: " << Bond->leftatom->Name << " and " << Bond->rightatom->Name << endl;
-    *LeftFragment = new molecule(elemente);
-    *RightFragment = new molecule(elemente);
-    // initialise marker list for all atoms
-    atom **AddedAtomListLeft = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::FragmentMoleculeByBond: **AddedAtomListLeft");
-    atom **AddedAtomListRight = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::FragmentMoleculeByBond: **AddedAtomListRight");
-    for (int i=0;i<AtomCount;i++) {
-      AddedAtomListLeft[i] = NULL;
-      AddedAtomListRight[i] = NULL;
-    }
-    bond **AddedBondListLeft = (bond **) Malloc(sizeof(bond *)*BondCount, "molecule::FragmentMoleculeByBond: **AddedBondListLeft");
-    bond **AddedBondListRight = (bond **) Malloc(sizeof(bond *)*BondCount, "molecule::FragmentMoleculeByBond: **AddedBondListRight");
-    for (int i=0;i<BondCount;i++) {
-      AddedBondListLeft[i] = NULL;
-      AddedBondListRight[i] = NULL;
-    }
-  
-    // tag and add all atoms that have to be included
-    *out << Verbose(1) << "Adding BFS-wise on left hand side with Bond Order " << NoNonBonds-1 << "." << endl;
-    AddedAtomListLeft[Bond->leftatom->nr] = (*LeftFragment)->AddCopyAtom(Bond->leftatom);
-    BreadthFirstSearchAdd(out, *LeftFragment, AddedAtomListLeft, AddedBondListLeft, Bond->leftatom, Bond, 
-#ifdef ADDHYDROGEN
-    NoNonBonds
-#else
-    BondCount
-#endif    
-                                                                                                            , IsAngstroem, CutCyclic);
-    *out << Verbose(1) << "Adding BFS-wise on right hand side with Bond Order " << NoNonBonds-1 << "." << endl;
-    AddedAtomListRight[Bond->rightatom->nr] = (*RightFragment)->AddCopyAtom(Bond->rightatom);
-    BreadthFirstSearchAdd(out, *RightFragment, AddedAtomListRight, AddedBondListRight, Bond->rightatom, Bond, 
-#ifdef ADDHYDROGEN
-    NoNonBonds
-#else
-    BondCount
-#endif    
-                                                                                                            , IsAngstroem, CutCyclic);  
- 
-    // count atoms
-    (*LeftFragment)->CountAtoms(out);
-    (*RightFragment)->CountAtoms(out);
-    // free all and exit
-    Free((void **)&AddedAtomListLeft, "molecule::FragmentMoleculeByBond: **AddedAtomListLeft");
-    Free((void **)&AddedAtomListRight, "molecule::FragmentMoleculeByBond: **AddedAtomListRight");
-    Free((void **)&AddedBondListLeft, "molecule::FragmentMoleculeByBond: **AddedBondListLeft");
-    Free((void **)&AddedBondListRight, "molecule::FragmentMoleculeByBond: **AddedBondListRight");
-#ifdef ADDHYDROGEN
-  }
-#endif
-  *out << Verbose(0) << "End of FragmentMoleculeByBond." << endl;
-};
-
-/** Returns the given \a **FragmentList filled with molecules around each bond including up to \a BondDegree neighbours.
- * \param *out output stream for debugging
- * \param BondOrder neighbours on each side to be ...
- * \param IsAngstroem whether atomic coordination is in Angstroem (true) or atomic length/bohrradous (false)
- * \param CutCyclic whether to add cut cyclic bond or to saturate
- * \sa FragmentBottomUp(), molecule::AddBondOrdertoMolecule()
- */
-MoleculeListClass * molecule::GetEachBondFragmentOfOrder(ofstream *out, int BondOrder, bool IsAngstroem, enum CutCyclicBond CutCyclic)
-{
-  /// Allocate memory for Bond list and go through each bond and fragment molecule up to bond order and fill the list to be returned.
-  MoleculeListClass *FragmentList = NULL;
-  atom **AddedAtomList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::GetBondFragmentOfOrder: **AddedAtomList");
-  bond **AddedBondList = (bond **) Malloc(sizeof(bond *)*BondCount, "molecule::GetBondFragmentOfOrder: **AddedBondList");
-  bond *Binder = NULL;
-
-  *out << Verbose(0) << "Begin of GetEachBondFragmentOfOrder." << endl;
-#ifdef ADDHYDROGEN
-  if (NoNonBonds != 0) {
-    FragmentList = new MoleculeListClass(NoNonBonds, AtomCount);
-  } else {
-    *out << Verbose(1) << "NoNonBonds is " << NoNonBonds << ", can't allocate list." << endl;
-#else
-  if (BondCount != 0) {
-    FragmentList = new MoleculeListClass(BondCount, AtomCount);
-  } else {
-    *out << Verbose(1) << "BondCount is " << BondCount << ", can't allocate list." << endl;
-#endif
-  }
-  int No = 0;
-  Binder = first;
-  while (Binder->next != last) { // get each bond, NULL is returned if it is a H-H bond
-    Binder = Binder->next;
-#ifdef ADDHYDROGEN
-    if ((Binder->leftatom->type->Z != 1) && (Binder->rightatom->type->Z != 1))  // if both bond partners aren't hydrogen ...
-#endif
-      if ((CutCyclic == SaturateBond) || (!Binder->Cyclic)) {
-        *out << Verbose(1) << "Getting Fragment for Non-Hydrogen Bond: " << Binder->leftatom->Name << " and " << Binder->rightatom->Name << endl;
-        FragmentList->ListOfMolecules[No] = new molecule(elemente);
-        // initialise marker list for all atoms
-        for (int i=0;i<AtomCount;i++)
-          AddedAtomList[i] = NULL;
-        for (int i=0;i<BondCount;i++)
-          AddedBondList[i] = NULL;
-      
-        // add root bond as first bond (this is needed later on fragmenting)
-        *out << Verbose(1) << "Adding Root Bond " << *Binder << " and its atom." << endl;
-        AddedAtomList[Binder->leftatom->nr] = FragmentList->ListOfMolecules[No]->AddCopyAtom(Binder->leftatom);
-        AddedAtomList[Binder->rightatom->nr] = FragmentList->ListOfMolecules[No]->AddCopyAtom(Binder->rightatom);
-        AddedBondList[Binder->nr] = FragmentList->ListOfMolecules[No]->AddBond(AddedAtomList[Binder->leftatom->nr], AddedAtomList[Binder->rightatom->nr], Binder->BondDegree);
-    
-        // tag and add all atoms that have to be included
-        *out << Verbose(1) << "Adding BFS-wise on left hand side." << endl;
-        BreadthFirstSearchAdd(out, FragmentList->ListOfMolecules[No], AddedAtomList, AddedBondList, Binder->leftatom, NULL, BondOrder, IsAngstroem, CutCyclic);  
-        *out << Verbose(1) << "Adding BFS-wise on right hand side." << endl;
-        BreadthFirstSearchAdd(out, FragmentList->ListOfMolecules[No], AddedAtomList, AddedBondList, Binder->rightatom, NULL, BondOrder, IsAngstroem, CutCyclic);  
-        
-        // count atoms
-        FragmentList->ListOfMolecules[No]->CountAtoms(out);
-        FragmentList->TEList[No] = 1.;
-        *out << Verbose(1) << "GetBondFragmentOfOrder: " << Binder->nr << "th Fragment: " << FragmentList->ListOfMolecules[No] << "." << endl;
-        No++;
-      }
-  }
-  // free all lists
-  Free((void **)&AddedAtomList, "molecule::GetBondFragmentOfOrder: **AddedAtomList");
-  Free((void **)&AddedBondList, "molecule::GetBondFragmentOfOrder: **AddedBondList");
-  // output and exit
-  FragmentList->Output(out);
-  *out << Verbose(0) << "End of GetEachBondFragmentOfOrder." << endl;
-  return (FragmentList);
-}; 
-
 /** Adds atoms up to \a BondCount distance from \a *Root and notes them down in \a **AddedAtomList.
  * Gray vertices are always enqueued in an AtomStackClass FIFO queue, the rest is usual BFS with adding vertices found was
@@ -2603,7 +2433,6 @@
  * \param BondOrder maximum distance for vertices to add
  * \param IsAngstroem lengths are in angstroem or bohrradii
- * \param CutCyclic whether to cut cyclic bonds (means saturate on need with hydrogen) or to always add
  */  
-void molecule::BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem, enum CutCyclicBond CutCyclic)
+void molecule::BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem)
 {
   atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::BreadthFirstSearchAdd: **PredecessorList");
@@ -2650,5 +2479,5 @@
           ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
           *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " " << ((ColorList[OtherAtom->nr] == white) ? "white" : "lightgray") << ", its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl; 
-          if ((((ShortestPathList[OtherAtom->nr] < BondOrder) && (Binder != Bond)) || (Binder->Cyclic && (CutCyclic == KeepBond))) ) { // Check for maximum distance
+          if ((((ShortestPathList[OtherAtom->nr] < BondOrder) && (Binder != Bond))) ) { // Check for maximum distance
             *out << Verbose(3);
             if (AddedAtomList[OtherAtom->nr] == NULL) { // add if it's not been so far
@@ -2678,10 +2507,8 @@
             else if (ShortestPathList[OtherAtom->nr] >= BondOrder)
               *out << Verbose(3) << "Not Queueing, is out of Bond Count of " << BondOrder;
-            if ((Binder->Cyclic && (CutCyclic != KeepBond)))
-              *out << ", is part of a cyclic bond yet we don't keep them, saturating bond with Hydrogen." << endl;
             if (!Binder->Cyclic)
               *out << ", is not part of a cyclic bond, saturating bond with Hydrogen." << endl;
             if (AddedBondList[Binder->nr] == NULL) {
-              if ((AddedAtomList[OtherAtom->nr] != NULL) && (CutCyclic == KeepBond)) { // .. whether we add or saturate
+              if ((AddedAtomList[OtherAtom->nr] != NULL)) { // .. whether we add or saturate
                 AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
                 AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
@@ -2698,5 +2525,5 @@
           // This has to be a cyclic bond, check whether it's present ...
           if (AddedBondList[Binder->nr] == NULL) {
-            if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder) || (CutCyclic == KeepBond))) { 
+            if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) { 
               AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
               AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
@@ -3106,5 +2933,4 @@
   bool **UsedList;
   bond **BondsPerSPList;
-  double TEFactor;
   int *BondsPerSPCount;
 };
@@ -3165,5 +2991,5 @@
           *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
           TouchedList[TouchedIndex++] = OtherWalker->nr;  // note as added
-          FragmentSearch->FragmentSet->insert(OtherWalker->GetTrueFather()->nr);
+          FragmentSearch->FragmentSet->insert(OtherWalker->nr);
             //FragmentSearch->UsedList[OtherWalker->nr][i] = true;
           //}
@@ -3205,10 +3031,9 @@
         *out << Verbose(1+verbosity) << "Enough items on stack for a fragment!" << endl;
         // store fragment as a KeySet
-        *out << Verbose(2) << "Found a new fragment[" << FragmentSearch->FragmentCounter << "], indices are: ";
-        for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++) {
-          *out << (*runner)+1 << " "; 
-        }
+        *out << Verbose(2) << "Found a new fragment[" << FragmentSearch->FragmentCounter << "], local nr.s are: ";
+        for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
+          *out << (*runner) << " "; 
         InsertFragmentIntoGraph(out, FragmentSearch);
-        Removal = LookForRemovalCandidate(out, FragmentSearch->FragmentSet, FragmentSearch->ShortestPathList);
+        //Removal = LookForRemovalCandidate(out, FragmentSearch->FragmentSet, FragmentSearch->ShortestPathList);
         //Removal = StoreFragmentFromStack(out, FragmentSearch->Root, FragmentSearch->Leaflet, FragmentSearch->FragmentStack, FragmentSearch->ShortestPathList,FragmentSearch->Labels, &FragmentSearch->FragmentCounter, FragmentSearch->configuration);
       }
@@ -3218,8 +3043,12 @@
       for(int j=0;j<TouchedIndex;j++) {
         Removal = TouchedList[j];
-        *out << Verbose(2+verbosity) << "Removing item nr. " << Removal+1 << " from snake stack." << endl;
+        *out << Verbose(2+verbosity) << "Removing item nr. " << Removal << " from snake stack." << endl;
         FragmentSearch->FragmentSet->erase(Removal);
         TouchedList[j] = -1;
       }
+      *out << Verbose(2) << "Remaining local nr.s on snake stack are: ";
+      for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
+        *out << (*runner) << " "; 
+      *out << endl;
       TouchedIndex = 0; // set Index to 0 for list of atoms added on this level
     } else {
@@ -3228,42 +3057,40 @@
   }
   Free((void **)&TouchedList, "molecule::SPFragmentGenerator: *TouchedList");  
-  *out << Verbose(1+verbosity) << "End of SPFragmentGenerator " << RootDistance << " away from Root " << *FragmentSearch->Root << " and SubOrder is " << SubOrder << "." << endl;
-};
-
-/** Creates a list of all unique fragments of certain vertex size from a given graph \a Fragment in the context of \a this molecule.
+  *out << Verbose(1+verbosity) << "End of SPFragmentGenerator, " << RootDistance << " away from Root " << *FragmentSearch->Root << " and SubOrder is " << SubOrder << "." << endl;
+};
+
+/** Creates a list of all unique fragments of certain vertex size from a given graph \a Fragment for a given root vertex in the context of \a this molecule.
  * 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 number of vertices
- * \param *ListOfGraph Graph structure to insert found fragments into
- * \param Fragment Restricted vertex set to use in context of molecule
- * \param TEFactor TEFactor to store in graphlist in the end
- * \param *configuration configuration needed for IsAngstroem
+ * \param Order bond order (limits BFS exploration and "number of digits" in power set generation 
+ * \param *ReturnKeySets Graph structure to insert found keysets/fragments into
+ * \param RestrictedKeySet Restricted vertex set to use in context of molecule
+ * \param RootKeyNr Atom::nr of the atom acting as current fragment root 
  * \return number of inserted fragments
  * \note ShortestPathList in FragmentSearch structure is probably due to NumberOfAtomsSPLevel and SP not needed anymore
  */
-int molecule::CreateListOfUniqueFragmentsOfOrder(ofstream *out, int Order, Graph *ListOfGraph, KeySet Fragment, double TEFactor, config *configuration) 
-{
-  int SP, UniqueIndex, RootKeyNr, AtomKeyNr;
-  int *NumberOfAtomsSPLevel = (int *) Malloc(sizeof(int)*Order, "molecule::CreateListOfUniqueFragmentsOfOrder: *SPLevelCount");
+int molecule::PowerSetGenerator(ofstream *out, int Order, Graph *ReturnKeySets, KeySet RestrictedKeySet, int RootKeyNr) 
+{
+  int SP, UniqueIndex, AtomKeyNr;
+  int *NumberOfAtomsSPLevel = (int *) Malloc(sizeof(int)*Order, "molecule::PowerSetGenerator: *SPLevelCount");
   atom *Walker = NULL, *OtherWalker = NULL;
   bond *Binder = NULL;
   bond **BondsList = NULL;
-  KeyStack RootStack;
   KeyStack AtomStack;
-  atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
+  atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::PowerSetGenerator: **PredecessorList");
   KeySet::iterator runner;
+  //int Count = RestrictedKeySet.size();
 
   // initialise the fragments structure
   struct UniqueFragments FragmentSearch;
-  FragmentSearch.BondsPerSPList = (bond **) Malloc(sizeof(bond *)*Order*2, "molecule::CreateListOfUniqueFragmentsOfOrder: ***BondsPerSPList");
-  FragmentSearch.BondsPerSPCount = (int *) Malloc(sizeof(int)*Order, "molecule::CreateListOfUniqueFragmentsOfOrder: *BondsPerSPCount");
-  FragmentSearch.ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
-  FragmentSearch.Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
+  FragmentSearch.BondsPerSPList = (bond **) Malloc(sizeof(bond *)*Order*2, "molecule::PowerSetGenerator: ***BondsPerSPList");
+  FragmentSearch.BondsPerSPCount = (int *) Malloc(sizeof(int)*Order, "molecule::PowerSetGenerator: *BondsPerSPCount");
+  FragmentSearch.ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *ShortestPathList");
+  FragmentSearch.Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *Labels");
   FragmentSearch.FragmentCounter = 0;
   FragmentSearch.FragmentSet = new KeySet;
-  FragmentSearch.configuration = configuration;
-  FragmentSearch.TEFactor = TEFactor;
-  FragmentSearch.Leaflet = ListOfGraph;      // set to insertion graph
+  FragmentSearch.Leaflet = ReturnKeySets;      // set to insertion graph
+  FragmentSearch.Root = FindAtom(RootKeyNr);
   for (int i=0;i<AtomCount;i++) {
     FragmentSearch.Labels[i] = -1;
@@ -3280,237 +3107,186 @@
   
   *out << endl;
-  *out << Verbose(0) << "Begin of CreateListOfUniqueFragmentsOfOrder with order " << Order << "." << endl;
-
-  RootStack.clear();  
-  // find first root candidates  
-  runner = Fragment.begin();
-  Walker = NULL;
-  while ((Walker == NULL) && (runner != Fragment.end())) { // search for first non-hydrogen atom
-    Walker = FindAtom((*runner));
-    *out << Verbose(4) << "Current Root candidate is " << Walker->Name << "." << endl;
+  *out << Verbose(0) << "Begin of PowerSetGenerator with order " << Order << " at Root " << *FragmentSearch.Root << "." << endl;
+
+  UniqueIndex = 0;
+  if (FragmentSearch.Labels[RootKeyNr] == -1)
+    FragmentSearch.Labels[RootKeyNr] = UniqueIndex++;
+  FragmentSearch.ShortestPathList[RootKeyNr] = 0;
+  // prepare the atom stack counters (number of atoms with certain SP on stack)
+  for (int i=0;i<Order;i++)
+    NumberOfAtomsSPLevel[i] = 0;
+  NumberOfAtomsSPLevel[0] = 1;  // for root
+  SP = -1;
+  *out << endl;
+  *out << Verbose(0) << "Starting BFS analysis ..." << endl;
+  // push as first on atom stack and goooo ...
+  AtomStack.clear();  
+  AtomStack.push_back(RootKeyNr);
+  *out << Verbose(0) << "Cleared AtomStack and pushed root as first item onto it." << endl;
+  // do a BFS search to fill the SP lists and label the found vertices
+  while (!AtomStack.empty()) {
+    // pop next atom
+    AtomKeyNr = AtomStack.front();
+    AtomStack.pop_front();
+    if (SP != -1)
+      NumberOfAtomsSPLevel[SP]--;
+    if ((SP == -1) || (NumberOfAtomsSPLevel[SP] == -1)) {
+      SP++;
+      NumberOfAtomsSPLevel[SP]--; // carry over "-1" to next level
+      *out << Verbose(1) << "New SP level reached: " << SP << ", creating new SP list with 0 item(s)";
+      if (SP > 0)
+        *out << ", old level closed with " << FragmentSearch.BondsPerSPCount[SP-1] << " item(s)." << endl;
+      else
+        *out << "." << endl;
+      FragmentSearch.BondsPerSPCount[SP] = 0;
+    } else {
+      *out << Verbose(1) << "Still " << NumberOfAtomsSPLevel[SP]+1 << " on this SP (" << SP << ") level, currently having " << FragmentSearch.BondsPerSPCount[SP] << " item(s)." << endl;
+    }
+    Walker = FindAtom(AtomKeyNr);
+    *out << Verbose(0) << "Current Walker is: " << *Walker << " with nr " << Walker->nr << " and label " << FragmentSearch.Labels[AtomKeyNr] << " and SP of " << SP << "." << endl;
+    // check for new sp level
+    // go through all its bonds
+    *out << Verbose(1) << "Going through all bonds of Walker." << endl;
+    for (int i=0;i<NumberOfBondsPerAtom[AtomKeyNr];i++) {
+      Binder = ListOfBondsPerAtom[AtomKeyNr][i];
+      OtherWalker = Binder->GetOtherAtom(Walker);
+      if ((RestrictedKeySet.find(OtherWalker->nr) != RestrictedKeySet.end())
 #ifdef ADDHYDROGEN
-    if (Walker->type->Z == 1) // skip hydrogen
-      Walker = NULL;
+       && (OtherWalker->type->Z != 1)
 #endif
-    runner++;
-  }
-  if (Walker != NULL)
-    RootStack.push_back(Walker->nr);
-  else
-    *out << Verbose(0) << "ERROR: Could not find an appropriate Root atom!" << endl;
-  
-  UniqueIndex = 0;
-  while (!RootStack.empty()) {
-    // Get Root and prepare
-    RootKeyNr = RootStack.front();
-    RootStack.pop_front();
-    FragmentSearch.Root = FindAtom(RootKeyNr);
-    if (FragmentSearch.Labels[RootKeyNr] == -1)
-      FragmentSearch.Labels[RootKeyNr] = UniqueIndex++;
-    FragmentSearch.ShortestPathList[RootKeyNr] = 0;
-    // prepare the atom stack counters (number of atoms with certain SP on stack)
-    for (int i=0;i<Order;i++)
-      NumberOfAtomsSPLevel[i] = 0;
-    NumberOfAtomsSPLevel[0] = 1;  // for root
-    SP = -1;
-    *out << endl;
-    *out << Verbose(0) << "Starting BFS analysis with current Root: " << *FragmentSearch.Root << "." << endl;
-    // push as first on atom stack and goooo ...
-    AtomStack.clear();  
-    AtomStack.push_back(RootKeyNr);
-    *out << Verbose(0) << "Cleared AtomStack and pushed root as first item onto it." << endl;
-    // do a BFS search to fill the SP lists and label the found vertices
-    while (!AtomStack.empty()) {
-      // pop next atom
-      AtomKeyNr = AtomStack.front();
-      AtomStack.pop_front();
-      if (SP != -1)
-        NumberOfAtomsSPLevel[SP]--;
-      if ((SP == -1) || (NumberOfAtomsSPLevel[SP] == -1)) {
-      ////if (SP < FragmentSearch.ShortestPathList[AtomKeyNr]) { // bfs has reached new SP level, hence allocate for new list
-        SP++;
-        NumberOfAtomsSPLevel[SP]--; // carry over "-1" to next level
-        ////SP = FragmentSearch.ShortestPathList[AtomKeyNr];
-        *out << Verbose(1) << "New SP level reached: " << SP << ", creating new SP list with 0 item(s)";
-        if (SP > 0)
-          *out << ", old level closed with " << FragmentSearch.BondsPerSPCount[SP-1] << " item(s)." << endl;
-        else
-          *out << "." << endl;
-        FragmentSearch.BondsPerSPCount[SP] = 0;
-      } else {
-        *out << Verbose(1) << "Still " << NumberOfAtomsSPLevel[SP]+1 << " on this SP (" << SP << ") level, currently having " << FragmentSearch.BondsPerSPCount[SP] << " item(s)." << endl;
-      }
-      Walker = FindAtom(AtomKeyNr);
-      *out << Verbose(0) << "Current Walker is: " << *Walker << " with nr " << Walker->nr << " and label " << FragmentSearch.Labels[AtomKeyNr] << " and SP of " << SP << "." << endl;
-      // check for new sp level
-      // go through all its bonds
-      *out << Verbose(1) << "Going through all bonds of Walker." << endl;
-      for (int i=0;i<NumberOfBondsPerAtom[AtomKeyNr];i++) {
-        Binder = ListOfBondsPerAtom[AtomKeyNr][i];
-        OtherWalker = Binder->GetOtherAtom(Walker);
-        if ((Fragment.find(OtherWalker->nr) != Fragment.end())
-#ifdef ADDHYDROGEN
-         && (OtherWalker->type->Z != 1)
-#endif
-                                                              ) {  // skip hydrogens and restrict to fragment
-          *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *Binder << "." << endl;
-          // set the label if not set (and push on root stack as well)
-          if (FragmentSearch.Labels[OtherWalker->nr] == -1) {
-            RootStack.push_back(OtherWalker->nr);
-            FragmentSearch.Labels[OtherWalker->nr] = UniqueIndex++;
-            *out << Verbose(3) << "Set label to " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
-          } else {
-            *out << Verbose(3) << "Label is already " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
-          }
-          if ((OtherWalker != PredecessorList[AtomKeyNr]) && (FragmentSearch.Labels[OtherWalker->nr] > FragmentSearch.Labels[RootKeyNr])) { // only pass through those with label bigger than Root's
-            // set shortest path if not set or longer
-            //if ((FragmentSearch.ShortestPathList[OtherWalker->nr] == -1) || (FragmentSearch.ShortestPathList[OtherWalker->nr] > FragmentSearch.ShortestPathList[AtomKeyNr]))  {
-              FragmentSearch.ShortestPathList[OtherWalker->nr] = SP+1;
-              *out << Verbose(3) << "Set Shortest Path to " << FragmentSearch.ShortestPathList[OtherWalker->nr] << "." << endl;
-            //} else {
-            //  *out << Verbose(3) << "Shortest Path is already " << FragmentSearch.ShortestPathList[OtherWalker->nr] << "." << endl;
-            //}
-            if (FragmentSearch.ShortestPathList[OtherWalker->nr] < Order) { // only pass through those within reach of Order
-              // push for exploration on stack (only if the SP of OtherWalker is longer than of Walker! (otherwise it has been added already!)
-              if (FragmentSearch.ShortestPathList[OtherWalker->nr] > SP) {
-                if (FragmentSearch.ShortestPathList[OtherWalker->nr] < (Order-1)) {
-                  *out << Verbose(3) << "Putting on atom stack for further exploration." << endl;
-                  PredecessorList[OtherWalker->nr] = Walker;    // note down the one further up the exploration tree
-                  AtomStack.push_back(OtherWalker->nr);
-                  NumberOfAtomsSPLevel[FragmentSearch.ShortestPathList[OtherWalker->nr]]++;
-                } else {
-                  *out << Verbose(3) << "Not putting on atom stack, is at end of reach." << 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]);
-                FragmentSearch.BondsPerSPCount[SP]++;
-                *out << Verbose(3) << "Added its bond to SP list, having now " << FragmentSearch.BondsPerSPCount[SP] << " item(s)." << endl;
+                                                            ) {  // skip hydrogens and restrict to fragment
+        *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *Binder << "." << endl;
+        // set the label if not set (and push on root stack as well)
+        if (FragmentSearch.Labels[OtherWalker->nr] == -1) {
+          FragmentSearch.Labels[OtherWalker->nr] = UniqueIndex++;
+          *out << Verbose(3) << "Set label to " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
+        } else {
+          *out << Verbose(3) << "Label is already " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
+        }
+        if ((OtherWalker != PredecessorList[AtomKeyNr]) && (FragmentSearch.Labels[OtherWalker->nr] > FragmentSearch.Labels[RootKeyNr])) { // only pass through those with label bigger than Root's
+          FragmentSearch.ShortestPathList[OtherWalker->nr] = SP+1;
+          *out << Verbose(3) << "Set Shortest Path to " << FragmentSearch.ShortestPathList[OtherWalker->nr] << "." << endl;
+          if (FragmentSearch.ShortestPathList[OtherWalker->nr] < Order) { // only pass through those within reach of Order
+            // push for exploration on stack (only if the SP of OtherWalker is longer than of Walker! (otherwise it has been added already!)
+            if (FragmentSearch.ShortestPathList[OtherWalker->nr] > SP) {
+              if (FragmentSearch.ShortestPathList[OtherWalker->nr] < (Order-1)) {
+                *out << Verbose(3) << "Putting on atom stack for further exploration." << endl;
+                PredecessorList[OtherWalker->nr] = Walker;    // note down the one further up the exploration tree
+                AtomStack.push_back(OtherWalker->nr);
+                NumberOfAtomsSPLevel[FragmentSearch.ShortestPathList[OtherWalker->nr]]++;
               } else {
-                *out << Verbose(3) << "Not putting on atom stack, nor adding bond, as " << *OtherWalker << "'s SP is shorter than Walker's." << endl;
+                *out << Verbose(3) << "Not putting on atom stack, is at end of reach." << endl;
               }
-            } else {
-              *out << Verbose(3) << "Not continuing, as " << *OtherWalker << " is out of reach of order " << Order << "." << endl;
-            }
-          } else {
-            *out << Verbose(3) << "Not passing on, as label of " << *OtherWalker << " " << FragmentSearch.Labels[OtherWalker->nr] << " is smaller than that of Root " << FragmentSearch.Labels[RootKeyNr] << " or this is my predecessor." << endl;
-          }
-        } else {
-          *out << Verbose(2) << "Is not in the Fragment or skipping hydrogen " << *OtherWalker << "." << endl;
-        }
-      }
-    }
-    // reset predecessor list
-    for(int i=0;i<Order;i++) {
-      Binder = FragmentSearch.BondsPerSPList[2*i];
-      *out << Verbose(1) << "Current SP level is " << i << "." << endl;
-      while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
+              // 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]);
+              FragmentSearch.BondsPerSPCount[SP]++;
+              *out << Verbose(3) << "Added its bond to SP list, having now " << FragmentSearch.BondsPerSPCount[SP] << " item(s)." << endl;
+            } else *out << Verbose(3) << "Not putting on atom stack, nor adding bond, as " << *OtherWalker << "'s SP is shorter than Walker's." << endl;
+          } else *out << Verbose(3) << "Not continuing, as " << *OtherWalker << " is out of reach of order " << Order << "." << endl;
+        } else *out << Verbose(3) << "Not passing on, as label of " << *OtherWalker << " " << FragmentSearch.Labels[OtherWalker->nr] << " is smaller than that of Root " << FragmentSearch.Labels[RootKeyNr] << " or this is my predecessor." << endl;
+      } else *out << Verbose(2) << "Is not in the retstricted keyset or skipping hydrogen " << *OtherWalker << "." << endl;
+    }
+  }
+  // reset predecessor list
+  for(int i=0;i<Order;i++) {
+    Binder = FragmentSearch.BondsPerSPList[2*i];
+    *out << Verbose(1) << "Current SP level is " << i << "." << endl;
+    while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
+      Binder = Binder->next;
+      PredecessorList[Binder->rightatom->nr] = NULL;  // By construction "OtherAtom" is always Bond::rightatom!
+    } 
+  }
+  *out << endl;
+  
+  // outputting all list for debugging
+  *out << Verbose(0) << "Printing all found lists." << endl;
+  for(int i=0;i<Order;i++) {
+    Binder = FragmentSearch.BondsPerSPList[2*i];
+    *out << Verbose(1) << "Current SP level is " << i << "." << endl;
+    while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
+      Binder = Binder->next;
+      *out << Verbose(2) << *Binder << endl;
+    } 
+  }
+  
+  // creating fragments with the found edge sets
+  SP = 0;
+  for(int i=0;i<Order;i++) { // sum up all found edges
+    Binder = FragmentSearch.BondsPerSPList[2*i];
+    while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
+      Binder = Binder->next;
+      SP ++; 
+    }
+  }
+  *out << Verbose(0) << "Total number of edges is " << SP << "." << endl;
+  if (SP >= (Order-1)) {
+    // start with root (push on fragment stack)
+    *out << Verbose(0) << "Starting fragment generation with " << *FragmentSearch.Root << ", local nr is " << FragmentSearch.Root->nr << "." << endl;
+    FragmentSearch.FragmentSet->clear();  
+    FragmentSearch.FragmentSet->insert(FragmentSearch.Root->nr);
+    
+    if (FragmentSearch.FragmentSet->size() == (unsigned int) Order) {
+      *out << Verbose(0) << "Enough items on stack already for a fragment!" << endl;
+      // store fragment as a KeySet
+      *out << Verbose(2) << "Found a new fragment[" << FragmentSearch.FragmentCounter << "], local nr.s are: ";
+      for(KeySet::iterator runner = FragmentSearch.FragmentSet->begin(); runner != FragmentSearch.FragmentSet->end(); runner++) {
+        *out << (*runner) << " "; 
+      }
+      *out << endl;
+      InsertFragmentIntoGraph(out, &FragmentSearch);
+    } else {
+      *out << Verbose(0) << "Preparing subset for this root and calling generator." << endl;
+      // prepare the subset and call the generator
+      BondsList = (bond **) Malloc(sizeof(bond *)*FragmentSearch.BondsPerSPCount[0], "molecule::PowerSetGenerator: **BondsList");
+      Binder = FragmentSearch.BondsPerSPList[0];
+      for(int i=0;i<FragmentSearch.BondsPerSPCount[0];i++) {
         Binder = Binder->next;
-        PredecessorList[Binder->rightatom->nr] = NULL;  // By construction "OtherAtom" is always Bond::rightatom!
-      } 
-    }
-    *out << endl;
-    *out << Verbose(0) << "Printing all found lists." << endl;
-    // outputting all list for debugging
-    for(int i=0;i<Order;i++) {
-      Binder = FragmentSearch.BondsPerSPList[2*i];
-      *out << Verbose(1) << "Current SP level is " << i << "." << endl;
-      while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
-        Binder = Binder->next;
-        *out << Verbose(2) << *Binder << endl;
-      } 
-    }
-    
-    // creating fragments with the found edge sets
-    SP = 0;
-    for(int i=0;i<Order;i++) { // sum up all found edges
-      Binder = FragmentSearch.BondsPerSPList[2*i];
-      while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
-        Binder = Binder->next;
-        SP ++; 
-      }
-    }
-    *out << Verbose(0) << "Total number of edges is " << SP << "." << endl;
-    if (SP >= (Order-1)) {
-      // start with root (push on fragment stack)
-      *out << Verbose(0) << "Starting fragment generation with " << *FragmentSearch.Root << "." << endl;
-      FragmentSearch.FragmentSet->clear();  
-      FragmentSearch.FragmentSet->insert(FragmentSearch.Root->GetTrueFather()->nr);
-      
-      if (FragmentSearch.FragmentSet->size() == (unsigned int) Order) {
-        *out << Verbose(0) << "Enough items on stack already for a fragment!" << endl;
-        // store fragment as a KeySet
-        *out << Verbose(2) << "Found a new fragment[" << FragmentSearch.FragmentCounter << "], indices are: ";
-        for(KeySet::iterator runner = FragmentSearch.FragmentSet->begin(); runner != FragmentSearch.FragmentSet->end(); runner++) {
-          *out << (*runner)+1 << " "; 
-        }
-        *out << endl;
-        InsertFragmentIntoGraph(out, &FragmentSearch);
-        //StoreFragmentFromStack(out, FragmentSearch.Root, FragmentSearch.Leaflet, FragmentSearch.FragmentStack, FragmentSearch.ShortestPathList,FragmentSearch.Labels, &FragmentSearch.FragmentCounter, FragmentSearch.configuration);
-      } else {
-        *out << Verbose(0) << "Preparing subset for this root and calling generator." << endl;
-        // prepare the subset and call the generator
-        BondsList = (bond **) Malloc(sizeof(bond *)*FragmentSearch.BondsPerSPCount[0], "molecule::CreateListOfUniqueFragmentsOfOrder: **BondsList");
-        Binder = FragmentSearch.BondsPerSPList[0];
-        for(int i=0;i<FragmentSearch.BondsPerSPCount[0];i++) {
-          Binder = Binder->next;
-          BondsList[i] = Binder;
-        }
-        SPFragmentGenerator(out, &FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order-1);  
-        Free((void **)&BondsList, "molecule::CreateListOfUniqueFragmentsOfOrder: **BondsList");
-      }
-    } else {
-      *out << Verbose(0) << "Not enough total number of edges to build " << Order << "-body fragments." << endl;
-    }
-    
-    // remove root from stack
-    *out << Verbose(0) << "Removing root again from stack." << endl;
-    FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);    
-
-    // free'ing the bonds lists
-    *out << Verbose(0) << "Free'ing all found lists. and resetting index lists" << endl;
-    for(int i=0;i<Order;i++) {
-      *out << Verbose(1) << "Current SP level is " << i << ": ";
-      Binder = FragmentSearch.BondsPerSPList[2*i];
-      while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
-        Binder = Binder->next;
-        FragmentSearch.ShortestPathList[Binder->leftatom->nr] = -1;
-        FragmentSearch.ShortestPathList[Binder->rightatom->nr] = -1;
-      }
-      // delete added bonds
-      cleanup(FragmentSearch.BondsPerSPList[2*i], FragmentSearch.BondsPerSPList[2*i+1]);
-      // also start and end node
-      *out << "cleaned." << endl;
-    }
-  }
-  
+        BondsList[i] = Binder;
+      }
+      SPFragmentGenerator(out, &FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order-1);  
+      Free((void **)&BondsList, "molecule::PowerSetGenerator: **BondsList");
+    }
+  } else {
+    *out << Verbose(0) << "Not enough total number of edges to build " << Order << "-body fragments." << endl;
+  }
+
+/*  // as FragmentSearch structure is used only once, we don't have to clean it anymore
+  // remove root from stack
+  *out << Verbose(0) << "Removing root again from stack." << endl;
+  FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);    
+
+  // free'ing the bonds lists
+  *out << Verbose(0) << "Free'ing all found lists. and resetting index lists" << endl;
+  for(int i=0;i<Order;i++) {
+    *out << Verbose(1) << "Current SP level is " << i << ": ";
+    Binder = FragmentSearch.BondsPerSPList[2*i];
+    while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
+      Binder = Binder->next;
+      // *out << "Removing atom " << Binder->leftatom->nr << " and " << Binder->rightatom->nr << "." << endl; // make sure numbers are local
+      FragmentSearch.ShortestPathList[Binder->leftatom->nr] = -1;
+      FragmentSearch.ShortestPathList[Binder->rightatom->nr] = -1;
+    }
+    // delete added bonds
+    cleanup(FragmentSearch.BondsPerSPList[2*i], FragmentSearch.BondsPerSPList[2*i+1]);
+    // also start and end node
+    *out << "cleaned." << endl;
+  }
+*/
   // free allocated memory
-  Free((void **)&NumberOfAtomsSPLevel, "molecule::CreateListOfUniqueFragmentsOfOrder: *SPLevelCount");
-  Free((void **)&PredecessorList, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
+  Free((void **)&NumberOfAtomsSPLevel, "molecule::PowerSetGenerator: *SPLevelCount");
+  Free((void **)&PredecessorList, "molecule::PowerSetGenerator: **PredecessorList");
   for(int i=0;i<Order;i++) { // delete start and end of each list
     delete(FragmentSearch.BondsPerSPList[2*i]);
     delete(FragmentSearch.BondsPerSPList[2*i+1]);
   }
-  Free((void **)&FragmentSearch.BondsPerSPList, "molecule::CreateListOfUniqueFragmentsOfOrder: ***BondsPerSPList");
-  Free((void **)&FragmentSearch.BondsPerSPCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *BondsPerSPCount");
-  Free((void **)&FragmentSearch.ShortestPathList, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
-  Free((void **)&FragmentSearch.Labels, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
+  Free((void **)&FragmentSearch.BondsPerSPList, "molecule::PowerSetGenerator: ***BondsPerSPList");
+  Free((void **)&FragmentSearch.BondsPerSPCount, "molecule::PowerSetGenerator: *BondsPerSPCount");
+  Free((void **)&FragmentSearch.ShortestPathList, "molecule::PowerSetGenerator: *ShortestPathList");
+  Free((void **)&FragmentSearch.Labels, "molecule::PowerSetGenerator: *Labels");
   delete(FragmentSearch.FragmentSet);
   
-//  // gather all the leaves together
-//  *out << Verbose(0) << "Copying all fragments into MoleculeList structure." << endl; 
-//  FragmentList = new Graph;
-//  UniqueIndex = 0;
-//  while ((FragmentSearch.Leaflet != NULL) && (UniqueIndex < FragmentSearch.FragmentCounter))  {
-//    FragmentList->insert();
-//    FragmentList->ListOfMolecules[UniqueIndex++] = FragmentSearch.Leaflet->Leaf;
-//    FragmentSearch.Leaflet->Leaf = NULL; // prevent molecule from being removed
-//    TempLeaf = FragmentSearch.Leaflet;
-//    FragmentSearch.Leaflet = FragmentSearch.Leaflet->previous;
-//    delete(TempLeaf);
-//  };
-
   // return list  
-  *out << Verbose(0) << "End of CreateListOfUniqueFragmentsOfOrder." << endl;
+  *out << Verbose(0) << "End of PowerSetGenerator." << endl;
   return FragmentSearch.FragmentCounter;
 };
@@ -3664,5 +3440,5 @@
   GraphTestPair testGraphInsert;
 
-  testGraphInsert = Fragment->Leaflet->insert(GraphPair (*Fragment->FragmentSet,pair<int,double>(Fragment->FragmentCounter,Fragment->TEFactor)));  // store fragment number and current factor
+  testGraphInsert = Fragment->Leaflet->insert(GraphPair (*Fragment->FragmentSet,pair<int,double>(Fragment->FragmentCounter,1)));  // store fragment number and current factor
   if (testGraphInsert.second) {
     *out << Verbose(2) << "KeySet " << Fragment->FragmentCounter << " successfully inserted." << endl;
@@ -3670,5 +3446,5 @@
   } else {
     *out << Verbose(2) << "KeySet " << Fragment->FragmentCounter << " failed to insert, present fragment is " << ((*(testGraphInsert.first)).second).first << endl;
-    ((*(testGraphInsert.first)).second).second += Fragment->TEFactor;
+    ((*(testGraphInsert.first)).second).second ++;  // increase the "created" counter
     *out << Verbose(2) << "New factor is " << ((*(testGraphInsert.first)).second).second << "." << endl;
   }
@@ -3687,4 +3463,5 @@
  * \param graph1 first (dest) graph
  * \param graph2 second (source) graph
+ * \param *counter keyset counter that gets increased 
  */
 inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter)
@@ -3705,16 +3482,18 @@
 
 
-/** Creates truncated BOSSANOVA expansion up to order \a k.
+/** Performs BOSSANOVA decomposition at selected sites, increasing the cutoff by one at these sites.
+ * 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 ANOVAOrder ANOVA expansion is truncated above this order
- * \param *configuration configuration for writing config files for each fragment
  * \return pointer to Graph list
  */
-Graph * molecule::FragmentBOSSANOVA(ofstream *out, int ANOVAOrder, config *configuration)
+Graph * molecule::FragmentBOSSANOVA(ofstream *out, KeyStack &RootStack)
 {
   Graph *FragmentList = NULL, ***FragmentLowerOrdersList = NULL;
-  //MoleculeListClass *FragmentMoleculeList = NULL;
   int Order, NumLevels, NumMolecules, TotalNumMolecules = 0, *NumMoleculesOfOrder = NULL;
   int counter = 0;
+  int UpgradeCount = RootStack.size();
+  KeyStack FragmentRootStack;
+  int RootKeyNr, RootNr;
   
   *out << Verbose(0) << "Begin of FragmentBOSSANOVA." << endl;
@@ -3722,8 +3501,8 @@
   // 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 = (int *) Malloc(sizeof(int)*ANOVAOrder, "molecule::FragmentBOSSANOVA: *NumMoleculesOfOrder");
-  FragmentLowerOrdersList = (Graph ***) Malloc(sizeof(Graph **)*ANOVAOrder, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");
-
-  // Construct the complete KeySet
+  NumMoleculesOfOrder = (int *) Malloc(sizeof(int)*UpgradeCount, "molecule::FragmentBOSSANOVA: *NumMoleculesOfOrder");
+  FragmentLowerOrdersList = (Graph ***) Malloc(sizeof(Graph **)*UpgradeCount, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");
+
+  // Construct the complete KeySet which we need for topmost level only (but for all Roots)
   atom *Walker = start;
   KeySet CompleteMolecule;
@@ -3733,54 +3512,73 @@
   } 
 
-  for (int BondOrder=1;BondOrder<=ANOVAOrder;BondOrder++) {
-    // 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)
-    NumLevels = 1 << (BondOrder-1); // (int)pow(2,BondOrder-1);
-    *out << Verbose(0) << "BondOrder is (" << BondOrder << "/" << ANOVAOrder << ") and NumLevels is " << NumLevels << "." << endl;
-  
+  // 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) {
+    RootKeyNr = RootStack.front();
+    RootStack.pop_front();
+    // increase adaptive order by one
+    Walker = FindAtom(RootKeyNr);
+    Walker->GetTrueFather()->AdaptiveOrder++;
+    Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;
+    
     // allocate memory for all lower level orders in this 1D-array of ptrs
-    FragmentLowerOrdersList[BondOrder-1] = (Graph **) Malloc(sizeof(Graph *)*NumLevels, "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
+    NumLevels = 1 << (Order); // (int)pow(2,Order);
+    FragmentLowerOrdersList[RootNr] = (Graph **) Malloc(sizeof(Graph *)*NumLevels, "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
   
     // create top order where nothing is reduced
     *out << Verbose(0) << "==============================================================================================================" << endl;
-    *out << Verbose(0) << "Creating list of unique fragments of Bond Order " << BondOrder << " itself." << endl;
+    *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", NumLevels is " << NumLevels << ", " << RootStack.size() << " Roots remaining." << endl;
+
     // Create list of Graphs of current Bond Order (i.e. F_{ij})
-    FragmentLowerOrdersList[BondOrder-1][0] =  new Graph;
-    NumMoleculesOfOrder[BondOrder-1] = CreateListOfUniqueFragmentsOfOrder(out, BondOrder, FragmentLowerOrdersList[BondOrder-1][0], CompleteMolecule, 1., configuration);
-    *out << Verbose(1) << "Number of resulting molecules is: " << NumMoleculesOfOrder[BondOrder-1] << "." << endl;
+    FragmentLowerOrdersList[RootNr][0] =  new Graph;
+    NumMoleculesOfOrder[RootNr] = PowerSetGenerator(out, Walker->AdaptiveOrder, FragmentLowerOrdersList[RootNr][0], CompleteMolecule, RootKeyNr);
+    *out << Verbose(1) << "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
     NumMolecules = 0;
     
-    // create lower order fragments
-    *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
-    Order = BondOrder;
-    for (int source=0;source<(NumLevels >> 1);source++) { // 1-terms don't need any more splitting, that's why only half is gone through (shift again)
-
-      // step down to next order at (virtual) boundary of powers of 2 in array
-      while (source >= (1 << (BondOrder-Order))) // (int)pow(2,BondOrder-Order))    
-        Order--;
-      *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
-      for (int SubOrder=Order;SubOrder>1;SubOrder--) {
-        int dest = source + (1 << (BondOrder-SubOrder));
-        *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
-        *out << Verbose(0) << "Current SubOrder is: " << SubOrder-1 << " with source " << source << " to destination " << dest << "." << endl;
-  
-        // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
-        //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[BondOrder-1][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
-        //NumMolecules = 0;
-        FragmentLowerOrdersList[BondOrder-1][dest] = new Graph;
-        for(Graph::iterator runner = (*FragmentLowerOrdersList[BondOrder-1][source]).begin();runner != (*FragmentLowerOrdersList[BondOrder-1][source]).end(); runner++) {
-          NumMolecules += CreateListOfUniqueFragmentsOfOrder(out,SubOrder-1, FragmentLowerOrdersList[BondOrder-1][dest], (*runner).first, -(*runner).second.second, configuration);
+    if ((NumLevels >> 1) > 0) {
+      // create lower order fragments
+      *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
+      Order = Walker->AdaptiveOrder;
+      for (int source=0;source<(NumLevels >> 1);source++) { // 1-terms don't need any more splitting, that's why only half is gone through (shift again)
+        // step down to next order at (virtual) boundary of powers of 2 in array
+        while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))    
+          Order--;
+        *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
+        for (int SubOrder=Order;SubOrder>1;SubOrder--) {
+          int dest = source + (1 << (Walker->AdaptiveOrder-SubOrder));
+          *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
+          *out << Verbose(0) << "Current SubOrder is: " << SubOrder-1 << " with source " << source << " to destination " << dest << "." << endl;
+    
+          // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
+          //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
+          //NumMolecules = 0;
+          FragmentLowerOrdersList[RootNr][dest] = new Graph;
+          for(Graph::iterator runner = (*FragmentLowerOrdersList[RootNr][source]).begin();runner != (*FragmentLowerOrdersList[RootNr][source]).end(); runner++) {
+            for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
+              FragmentList = new Graph();
+              PowerSetGenerator(out, SubOrder-1, FragmentList, (*runner).first, *sprinter);
+              // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
+              *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
+              InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], *FragmentList, &NumMolecules);
+              delete(FragmentList);
+            }
+          }
+          *out << Verbose(1) << "Number of resulting molecules for SubOrder " << SubOrder << " is: " << NumMolecules << "." << endl;
         }
-        *out << Verbose(1) << "Number of resulting molecules is: " << NumMolecules << "." << endl;
-      }
-    }
-    // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current BondOrder
-    //NumMoleculesOfOrder[BondOrder-1] = NumMolecules;
-    TotalNumMolecules += NumMoleculesOfOrder[BondOrder-1];
-    *out << Verbose(1) << "Number of resulting molecules for Order " << BondOrder << " is: " << NumMoleculesOfOrder[BondOrder-1] << "." << endl;    
-  }
+      }
+    }
+    // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
+    //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
+    TotalNumMolecules += NumMoleculesOfOrder[RootNr];
+    *out << Verbose(1) << "Number of resulting molecules for Order " << Walker->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
+    RootStack.push_back(RootKeyNr); // put back on stack
+    RootNr++;
+  }
+  *out << Verbose(0) << "==============================================================================================================" << endl;
   *out << Verbose(1) << "Total number of resulting molecules is: " << TotalNumMolecules << "." << endl;
+  *out << Verbose(0) << "==============================================================================================================" << endl;
   // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein) 
   // 5433222211111111
@@ -3789,15 +3587,20 @@
   // 21
   // 1
-  // Subsequently, we combine same orders into a single list (FragmentByOrderList) and reduce these by order
+  // Subsequently, we combine all into a single list (FragmentList)
 
   *out << Verbose(0) << "Combining the lists of all orders per order and finally into a single one." << endl;
   FragmentList = new Graph;
-  for (int BondOrder=1;BondOrder<=ANOVAOrder;BondOrder++) {
-    NumLevels = 1 << (BondOrder-1);
+  RootNr = 0;
+  while (!RootStack.empty()) {
+    RootKeyNr = RootStack.front();
+    RootStack.pop_front();
+    Walker = FindAtom(RootKeyNr); 
+    NumLevels = 1 << (Walker->AdaptiveOrder - 1);
     for(int i=0;i<NumLevels;i++) {
-      InsertGraphIntoGraph(out, *FragmentList, (*FragmentLowerOrdersList[BondOrder-1][i]), &counter);
-      delete(FragmentLowerOrdersList[BondOrder-1][i]);
-    }
-    Free((void **)&FragmentLowerOrdersList[BondOrder-1], "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
+      InsertGraphIntoGraph(out, *FragmentList, (*FragmentLowerOrdersList[RootNr][i]), &counter);
+      delete(FragmentLowerOrdersList[RootNr][i]);
+    }
+    Free((void **)&FragmentLowerOrdersList[RootNr], "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
+    RootNr++;
   }
   Free((void **)&FragmentLowerOrdersList, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");
@@ -3807,110 +3610,4 @@
   return FragmentList;
 };
-
-/** Fragments a molecule, taking \a BondDegree neighbours into accent.
- * First of all, we have to split up the molecule into bonds ranging out till \a BondDegree.
- * These fragments serve in the following as the basis the calculate the bond energy of the bond
- * they originated from. Thus, they are split up in a left and a right part, each calculated for
- * the total energy, including the fragment as a whole and then we get: 
- * E(fragment) - E(left) - E(right) = E(bond)
- * The splitting up is done via Breadth-First-Search, \sa BreadthFirstSearchAdd().
- * \param *out output stream for debugging
- * \param BondOrder up to how many neighbouring bonds a fragment contains
- * \param *configuration configuration for writing config files for each fragment
- * \param CutCyclic whether to add cut cyclic bond or to saturate
- * \return pointer to MoleculeListClass with all the fragments or NULL if something went wrong.
- */
-MoleculeListClass * molecule::FragmentBottomUp(ofstream *out, int BondOrder, config *configuration, enum CutCyclicBond CutCyclic)
-{
-  int Num;
-  MoleculeListClass *FragmentList =  NULL, **FragmentsList = NULL;
-  bond *Bond = NULL;
-
-  *out << Verbose(0) << "Begin of FragmentBottomUp." << endl; 
-  FragmentsList = (MoleculeListClass **) Malloc(sizeof(MoleculeListClass *)*2, "molecule::FragmentBottomUp: **FragmentsList");
-  *out << Verbose(0) << "Getting Atomic fragments." << endl;
-  FragmentsList[0] = GetAtomicFragments(out, AtomCount, configuration->GetIsAngstroem(), 1., CutCyclic);
-
-  // create the fragments including each one bond of the original molecule and up to \a BondDegree neighbours
-  *out << Verbose(0) << "Getting " << 
-#ifdef ADDHYDROGEN  
-  NoNonBonds 
-#else
-  BondCount
-#endif
-                                        << " Bond fragments." << endl;
-  FragmentList = GetEachBondFragmentOfOrder(out, BondOrder, configuration->GetIsAngstroem(), CutCyclic);
-
-  // check whether there are equal fragments by ReduceToUniqueOnes
-  FragmentList->ReduceToUniqueList(out, &cell_size[0], BondDistance);
-
-	*out << Verbose(0) << "Begin of Separating " << FragmentList->NumberOfMolecules << " Fragments into Bond pieces." << endl;
-	// as we have the list, we have to take each fragment split it relative to its originating
-	// bond into left and right and store these in a new list
-  *out << Verbose(2) << endl << "Allocating MoleculeListClass" << endl; 
-  FragmentsList[1] =  new MoleculeListClass(3*FragmentList->NumberOfMolecules, AtomCount); // for each molecule the whole and its left and right part 
-  *out << Verbose(2) << "Creating TEList." << endl; 
-  // and create TE summation for these bond energy approximations (bond = whole - left - right)
-  for(int i=0;i<FragmentList->NumberOfMolecules;i++) {
-    // make up factors to regain total energy of whole molecule
-    FragmentsList[1]->TEList[3*i] = FragmentList->TEList[i];       // bond energy is 1 * whole 
-    FragmentsList[1]->TEList[3*i+1] = -FragmentList->TEList[i];    // - 1. * left part
-    FragmentsList[1]->TEList[3*i+2] = -FragmentList->TEList[i];    // - 1. * right part
-
-    // shift the pointer on whole molecule to new list in order to avoid that this molecule is deleted on deconstructing FragmentList
-    FragmentsList[1]->ListOfMolecules[3*i] = FragmentList->ListOfMolecules[i]; 
-    *out << Verbose(2) << "shifting whole fragment pointer for fragment " << FragmentList->ListOfMolecules[i] << " -> " <<  FragmentsList[1]->ListOfMolecules[3*i] << "." << endl;
-    // create bond matrix and count bonds
-    *out << Verbose(2) << "Updating bond list for fragment " << i << " [" << FragmentList << "]: " << FragmentList->ListOfMolecules[i] << endl; 
-    // create list of bonds per atom for this fragment (atoms were counted above)
-    FragmentsList[1]->ListOfMolecules[3*i]->CreateListOfBondsPerAtom(out);
-     
-    *out << Verbose(0) << "Getting left & right fragments for fragment " << i << "." << endl;
-    // the bond around which the fragment has been setup is always the first by construction (bond partners are first added atoms)
-    Bond = FragmentsList[1]->ListOfMolecules[3*i]->first->next; // is the bond between atom 0 and another in the middle
-    FragmentsList[1]->ListOfMolecules[3*i]->FragmentMoleculeByBond(out, Bond, &(FragmentsList[1]->ListOfMolecules[3*i+1]), &(FragmentsList[1]->ListOfMolecules[3*i+2]), configuration->GetIsAngstroem(), CutCyclic);
-    if ((FragmentsList[1]->ListOfMolecules[3*i+1] == NULL) || (FragmentsList[1]->ListOfMolecules[3*i+2] == NULL)) {
-      *out << Verbose(2) << "Left and/or Right Fragment is NULL!" << endl;
-    } else {
-      *out << Verbose(3) << "Left Fragment is " <<  FragmentsList[1]->ListOfMolecules[3*i+1] << ": " << endl;
-      FragmentsList[1]->ListOfMolecules[3*i+1]->Output(out);
-      *out << Verbose(3) << "Right Fragment is " <<  FragmentsList[1]->ListOfMolecules[3*i+2] << ": " << endl;
-      FragmentsList[1]->ListOfMolecules[3*i+2]->Output(out);
-      *out << endl;
-    }
-    // remove in old list, so that memory for this molecule is not free'd on final delete of this list
-    FragmentList->ListOfMolecules[i] = NULL;
-  }
-	*out << Verbose(0) << "End of Separating Fragments into Bond pieces." << endl;
-  delete(FragmentList);
-  FragmentList = NULL;
-
-  // combine atomic and bond list
-  FragmentList = new MoleculeListClass(FragmentsList[0]->NumberOfMolecules + FragmentsList[1]->NumberOfMolecules, AtomCount);
-  Num = 0;
-  for(int i=0;i<2;i++) {
-    for(int j=0;j<FragmentsList[i]->NumberOfMolecules;j++) {
-      // transfer molecule
-      FragmentList->ListOfMolecules[Num] = FragmentsList[i]->ListOfMolecules[j];
-      FragmentsList[i]->ListOfMolecules[j] = NULL;
-      // transfer TE factor
-      FragmentList->TEList[Num] = FragmentsList[i]->TEList[j];
-      Num++;
-    }
-    delete(FragmentsList[i]);
-    FragmentsList[i] = NULL;
-  }
-  *out << Verbose(2) << "Memory cleanup and return with filled list." << endl; 
-  Free((void **)&FragmentsList, "molecule::FragmentBottomUp: **FragmentsList");
-
-  // reducing list
-  FragmentList->ReduceToUniqueList(out, &cell_size[0], BondDistance);
-	
-	// write configs for all fragements (are written in FragmentMolecule)
-  // free FragmentList
-  *out << Verbose(0) << "End of FragmentBottomUp." << endl; 
-  return FragmentList;
-};
-
 
 /** Comparision function for GSL heapsort on distances in two molecules.
Index: src/molecules.hpp
===================================================================
--- src/molecules.hpp	(revision 260934983a76c25770b6a1f5706e2fc4eb1eecec)
+++ src/molecules.hpp	(revision db942effae4d25e57224c90078ac288434a4a783)
@@ -356,4 +356,5 @@
     int LowpointNr; //!< needed in DepthFirstSearchAnalysis() to detect nonseparable components, is the lowest possible number of an atom to reach via tree edges only followed by at most one back edge.
     bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, given in DepthFirstSearchAnalysis()
+    unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
   
   atom();
@@ -519,19 +520,15 @@
   
   /// Fragment molecule by two different approaches:
-  void FragmentMolecule(ofstream *out, int BottomUpOrder, int TopDownOrder, enum BondOrderScheme Scheme, config *configuration, enum CutCyclicBond CutCyclic);
-  bool StoreAdjacencyToFile(ofstream *out, char *path, int bondorder);
-  bool CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms, int bondorder);
+  void FragmentMolecule(ofstream *out, int Order, config *configuration);
+  bool StoreAdjacencyToFile(ofstream *out, char *path);
+  bool CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms);
+  bool ParseOrderAtSiteFromFile(ofstream *out, char *path);
+  bool StoreOrderAtSiteFile(ofstream *out, char *path);
   bool ParseKeySetFile(ofstream *out, char *filename, atom **ListOfAtoms, MoleculeListClass *&FragmentList, bool IsAngstroem);
   bool ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet);
-  MoleculeListClass * GetAtomicFragments(ofstream *out, int NumberOfTopAtoms, bool IsAngstroem, double factor, enum CutCyclicBond CutCyclic);
-  void BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem, enum CutCyclicBond CutCyclic);
-  /// -# BottomUp  
-  MoleculeListClass * FragmentBottomUp(ofstream *out, int BondOrder, config *configuration, enum CutCyclicBond CutCyclic);
-  MoleculeListClass * GetEachBondFragmentOfOrder(ofstream *out, int BondOrder, bool IsAngstroem, enum CutCyclicBond CutCyclic);
-  /// -# TopDown
-  void FragmentMoleculeByBond(ofstream *out, bond *Bond, molecule **LeftFragment, molecule **RightFragment, bool IsAngstroem, enum CutCyclicBond CutCyclic);
+  void BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem);
   /// -# BOSSANOVA
-  Graph * FragmentBOSSANOVA(ofstream *out, int ANOVAOrder, config *configuration);
-	int CreateListOfUniqueFragmentsOfOrder(ofstream *out, int Order, Graph *ListOfGraph, KeySet Fragment, double TEFactor, config *configuration);
+  Graph * FragmentBOSSANOVA(ofstream *out, KeyStack &RootStack);
+	int PowerSetGenerator(ofstream *out, int Order, Graph *ReturnKeySets, KeySet RestrictedKeySet, int RootKeyNr);
   bool BuildInducedSubgraph(ofstream *out, const molecule *Father);
   molecule * StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem);
@@ -560,5 +557,4 @@
     int NumberOfMolecules;        //!< Number of entries in \a **FragmentList and of to be returned one.
     int NumberOfTopAtoms;         //!< Number of atoms in the molecule from which all fragments originate
-    double *TEList;               //!< List of factors when summing over total energies of all fragment
     
   MoleculeListClass();
@@ -566,13 +562,6 @@
   ~MoleculeListClass();
 
-  MoleculeListClass * FragmentTopDown(ofstream *out, int BondDegree, double bonddistance, config *configuration, enum CutCyclicBond Saturation);
-
-  /// Reduced list to unique molecules.
-  int * GetMappingToUniqueFragments(ofstream *out, double threshold, double *cell_size, double celldistance);
-  int ReduceFragmentToUniqueOnes(ofstream *out, int *Map);
-  void ReduceToUniqueList(ofstream *out, double *cell_size, double celldistance);
-
   /// Output configs.
-  bool OutputConfigForListOfFragments(char *prefix, config *configuration, int *SortIndex);
+  bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
   void Output(ofstream *out);
   
Index: src/parser.cpp
===================================================================
--- src/parser.cpp	(revision 260934983a76c25770b6a1f5706e2fc4eb1eecec)
+++ src/parser.cpp	(revision db942effae4d25e57224c90078ac288434a4a783)
@@ -43,5 +43,5 @@
   stringstream line;
 
-  line << argv[1] << FRAGMENTPREFIX << TEFACTORSFILE;
+  //line << argv[1] << FRAGMENTPREFIX << TEFACTORSFILE;
   return FilePresent(line.str().c_str());
 };
@@ -106,5 +106,5 @@
   // count the number of matrices
   MatrixCounter = -1; // we count one too much
-  line << name << FRAGMENTPREFIX << TEFACTORSFILE;
+  line << name << FRAGMENTPREFIX << KEYSETFILE;
   input.open(line.str().c_str(), ios::in);
   if (input == NULL) {
@@ -417,33 +417,10 @@
 bool ForceMatrix::ParseIndices(char *name) 
 {
-  ifstream input;
-  char *FragmentNumber = NULL;
-  char filename[1023];
-  stringstream line;
-  
   cout << "Parsing force indices." << endl;
   Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "ForceMatrix::ParseIndices: **Indices");
-  line << name << FRAGMENTPREFIX << FORCESFILE;
-  input.open(line.str().c_str(), ios::in);
-  //cout << "Opening " << line.str() << " ... "  << input << endl;
-  if (input == NULL) {
-    cout << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
-    return false;
-  }
-  for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
-    // get the number of atoms for this fragment
-    stringstream line;
-    input.getline(filename, 1023);
-    line.str(filename);
-    // parse the values
+  for (int i=0;i<MatrixCounter;i++) {
     Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "ForceMatrix::ParseIndices: *Indices[]");
-    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
-    cout << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
-    Free((void **)&FragmentNumber, "ForceMatrix::ParseIndices: *FragmentNumber");
-    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
-      line >> Indices[i][j];
-      cout << " " << Indices[i][j];
-    }
-    cout << endl;
+    for(int j=0;j<RowCounter[i];j++)
+      Indices[i][j] = j;
   }
   Indices[MatrixCounter] = (int *) Malloc(sizeof(int)*RowCounter[MatrixCounter], "ForceMatrix::ParseIndices: *Indices[]");
@@ -451,5 +428,4 @@
     Indices[MatrixCounter][j] = j;
   }
-  input.close();
   return true;
 };
Index: src/vector.cpp
===================================================================
--- src/vector.cpp	(revision 260934983a76c25770b6a1f5706e2fc4eb1eecec)
+++ src/vector.cpp	(revision db942effae4d25e57224c90078ac288434a4a783)
@@ -78,8 +78,8 @@
   //vector Shifted, TranslationVector;
   vector TestVector;
-  *out << Verbose(1) << "Begin of KeepPeriodic." << endl;
-  *out << Verbose(2) << "Vector is: ";
-  Output(out);
-  *out << endl;
+//  *out << Verbose(1) << "Begin of KeepPeriodic." << endl;
+//  *out << Verbose(2) << "Vector is: ";
+//  Output(out);
+//  *out << endl;
   TestVector.CopyVector(this);
   TestVector.InverseMatrixMultiplication(matrix);
@@ -93,8 +93,8 @@
   TestVector.MatrixMultiplication(matrix);
   CopyVector(&TestVector);
-  *out << Verbose(2) << "New corrected vector is: ";
-  Output(out);
-  *out << endl;
-  *out << Verbose(1) << "End of KeepPeriodic." << endl;
+//  *out << Verbose(2) << "New corrected vector is: ";
+//  Output(out);
+//  *out << endl;
+//  *out << Verbose(1) << "End of KeepPeriodic." << endl;
 };
 
