Index: src/analyzer.cpp
===================================================================
--- src/analyzer.cpp	(revision bd65791c80eaf19f1799a7c9599f95606c4a4796)
+++ src/analyzer.cpp	(revision 2459b1b8b885ddb84c9c4feab2f53d9d08d718c0)
@@ -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/defs.hpp
===================================================================
--- src/defs.hpp	(revision bd65791c80eaf19f1799a7c9599f95606c4a4796)
+++ src/defs.hpp	(revision 2459b1b8b885ddb84c9c4feab2f53d9d08d718c0)
@@ -44,4 +44,6 @@
 #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"
Index: src/joiner.cpp
===================================================================
--- src/joiner.cpp	(revision bd65791c80eaf19f1799a7c9599f95606c4a4796)
+++ src/joiner.cpp	(revision 2459b1b8b885ddb84c9c4feab2f53d9d08d718c0)
@@ -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 bd65791c80eaf19f1799a7c9599f95606c4a4796)
+++ src/moleculelist.cpp	(revision 2459b1b8b885ddb84c9c4feab2f53d9d08d718c0)
@@ -133,4 +133,58 @@
 };
 
+
+/** Store force indices, i.e. the connection between the nuclear index in the total molecule config and the respective atom in fragment config.
+ * \param *out output stream for debugging
+ * \param *path path to file
+ * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
+ * \return true - file written successfully, false - writing failed
+ */
+bool MoleculeListClass::StoreForcesFile(ofstream *out, char *path, int *SortIndex)
+{
+  bool status = true;
+  ofstream ForcesFile;
+  stringstream line;
+  atom *Walker = NULL;
+  element *runner = NULL;
+
+  // open file for the force factors
+  *out << Verbose(1) << "Saving  force factors ... ";
+  line << path << "/" << FRAGMENTPREFIX << FORCESFILE;
+  ForcesFile.open(line.str().c_str(), ios::out);
+  if (ForcesFile != NULL) {
+    //cout << Verbose(1) << "Final AtomicForcesList: ";
+    //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 rea
+                //cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
+                ForcesFile <<  SortIndex[Walker->GetTrueFather()->nr] << "\t";
+                } else  // otherwise a -1 to indicate an added saturation hydrogen
+                  ForcesFile << "-1\t";
+              }
+            }
+          }
+      }
+      ForcesFile << endl;
+    }
+    ForcesFile.close();
+    *out << Verbose(1) << "done." << endl;
+  } else {
+    status = false;
+    *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
+  }
+  ForcesFile.close();
+  
+  return status;
+};
+
 /** Writes a config file for each molecule in the given \a **FragmentList.
  * \param *out output stream for debugging
@@ -216,28 +270,5 @@
     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();
-  *out << " done." << endl;
+  cout << " done." << endl;
   
   // printing final number
Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision bd65791c80eaf19f1799a7c9599f95606c4a4796)
+++ src/molecules.cpp	(revision 2459b1b8b885ddb84c9c4feab2f53d9d08d718c0)
@@ -1249,5 +1249,4 @@
   atom *Walker = NULL, *OtherWalker = NULL;
   int No, NoBonds;
-  bond *Binder = NULL;
   int NumberCells, divisor[NDIM], n[NDIM], N[NDIM], index, Index, j;
   molecule **CellList;
@@ -1365,5 +1364,5 @@
 		if (BondCount != 0) {
       NoCyclicBonds = 0;
-		  *out << Verbose(1) << "correct Bond degree of each bond" << endl;
+		  *out << Verbose(1) << "Correcting Bond degree of each bond ... ";
 		  do {
 		    No = 0; // No acts as breakup flag (if 1 we still continue)
@@ -1377,5 +1376,5 @@
             for(j=0;j<NumberOfBondsPerAtom[Walker->nr];j++)
               NoBonds += ListOfBondsPerAtom[Walker->nr][j]->BondDegree;
-	          *out << Verbose(3) << "Walker: " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+	          //*out << Verbose(3) << "Walker: " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
 	          if ((int)(Walker->type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check NoBonds of other atom
 	            // count valence of second partner
@@ -1383,5 +1382,5 @@
               for(j=0;j<NumberOfBondsPerAtom[OtherWalker->nr];j++)
                 NoBonds += ListOfBondsPerAtom[OtherWalker->nr][j]->BondDegree;
-              *out << Verbose(3) << "OtherWalker: " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+              //*out << Verbose(3) << "OtherWalker: " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
 	            if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) // increase bond degree by one
                 ListOfBondsPerAtom[Walker->nr][i]->BondDegree++;
@@ -1390,17 +1389,17 @@
 		    }
 		  } while (No);
-		
+		*out << " done." << endl;
 		} else
 			*out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
 	  *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << "." << endl;
 		
-	  // output bonds for debugging (if bond chain list was correctly installed)
-	  *out << Verbose(1) << endl << "From contents of bond chain list:";
-    Binder = first;
-    while(Binder->next != last) {
-      Binder = Binder->next;
-			*out << *Binder << "\t" << endl;
-    }
-    *out << endl;
+//	  // output bonds for debugging (if bond chain list was correctly installed)
+//	  *out << Verbose(1) << endl << "From contents of bond chain list:";
+//	  bond *Binder = first;
+//    while(Binder->next != last) {
+//      Binder = Binder->next;
+//			*out << *Binder << "\t" << endl;
+//    }
+//    *out << endl;
   } else
   	*out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl; 
@@ -1612,7 +1611,4 @@
   }
 
-  // further analysis of the found cycles (print rings, get minimum cycle length)
-  CyclicStructureAnalysis(out, MinimumRingSize);
-  
   // free all and exit
   delete(AtomStack);
@@ -1880,56 +1876,136 @@
 
 /** Parses the KeySet file and fills \a *FragmentList from the known molecule structure.
+ * Does two-pass scanning:
+ * -# Scans the keyset file and initialises a temporary graph
+ * -# Scans TEFactors file and sets the TEFactor of each key set in the temporary graph accordingly
+ * Finally, the temporary graph is inserted into the given \a FragmentList for return.
  * \param *out output stream for debugging
  * \param *path path to file
- * \param **ListOfAtoms allocated (molecule::AtomCount) and filled lookup table for ids (Atom::nr) to *Atom
- * \param *FragmentList NULL, filled on return
+ * \param *FragmentList empty, filled on return
  * \param IsAngstroem whether we have Ansgtroem or bohrradius
  * \return true - parsing successfully, false - failure on parsing (FragmentList will be NULL)
  */
-bool molecule::ParseKeySetFile(ofstream *out, char *path, atom **ListOfAtoms, MoleculeListClass *&FragmentList, bool IsAngstroem)
+bool molecule::ParseKeySetFile(ofstream *out, char *path, Graph *&FragmentList, bool IsAngstroem)
 {
   bool status = true;
-  ifstream KeySetFile;
+  ifstream InputFile;
   stringstream line;
+  GraphTestPair testGraphInsert;
+  int NumberOfFragments = 0;
+  double TEFactor;
   char *filename = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::ParseKeySetFile - filename");
   
-  if (FragmentList != NULL) { // check list pointer
-    cerr << "Error: FragmentList was not NULL as supposed to be, already atoms present therein?" << endl;
-    return false;
-  }
-  cout << Verbose(1) << "Parsing the KeySet file ... " << endl;
-  // open file and read
+  if (FragmentList == NULL) { // check list pointer
+    FragmentList = new Graph;
+  }
+  
+  // 1st pass: open file and read
+  *out << Verbose(1) << "Parsing the KeySet file ... " << endl;
   sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, KEYSETFILE);
-  KeySetFile.open(filename);
-  if (KeySetFile != NULL) {
+  InputFile.open(filename);
+  if (InputFile != NULL) {
     // each line represents a new fragment
-    int NumberOfFragments = 0;
     char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::ParseKeySetFile - *buffer");
-    // 1. scan through file to know number of fragments
-    while (!KeySetFile.eof()) {
-      KeySetFile.getline(buffer, MAXSTRINGSIZE);
-      if (strlen(buffer) > 0) // there is at least on possible number on the parsed line
-        NumberOfFragments++;
-    }
-    // 2. allocate the MoleculeListClass accordingly
-    FragmentList = new MoleculeListClass(NumberOfFragments, AtomCount);
-    // 3. Clear File, go to beginning and parse again, now adding found ids to each keyset and converting into molecules
-    KeySetFile.clear();
-    KeySetFile.seekg(ios::beg);
+    // 1. parse keysets and insert into temp. graph
+    while (!InputFile.eof()) {
+      InputFile.getline(buffer, MAXSTRINGSIZE);
+      KeySet CurrentSet;
+      if ((strlen(buffer) > 0) && (ScanBufferIntoKeySet(out, buffer, CurrentSet))) {  // if at least one valid atom was added, write config
+        testGraphInsert = FragmentList->insert(GraphPair (CurrentSet,pair<int,double>(NumberOfFragments++,1)));  // store fragment number and current factor
+        if (!testGraphInsert.second) {
+          cerr << "KeySet file must be corrupt as there are two equal key sets therein!" << endl;
+        }
+        //FragmentList->ListOfMolecules[NumberOfFragments++] = StoreFragmentFromKeySet(out, CurrentSet, IsAngstroem);
+      }
+    }
+    // 2. Free and done
+    InputFile.close();
+    InputFile.clear();
+    Free((void **)&buffer, "molecule::ParseKeySetFile - *buffer");
+    *out << Verbose(1) << "done." << endl;
+  } else {
+    *out << Verbose(1) << "File " << filename << " not found." << endl;
+    status = false;
+  }
+  
+  // 2nd pass: open TEFactors file and read
+  *out << Verbose(1) << "Parsing the TEFactors file ... " << endl;
+  sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, TEFACTORSFILE);
+  InputFile.open(filename);
+  if (InputFile != NULL) {
+    // 3. add found TEFactors to each keyset
     NumberOfFragments = 0;
-    while ((!KeySetFile.eof()) && (FragmentList->NumberOfMolecules > NumberOfFragments)) {
-      KeySetFile.getline(buffer, MAXSTRINGSIZE);
-      KeySet CurrentSet;
-      if ((strlen(buffer) > 0) && (ScanBufferIntoKeySet(out, buffer, CurrentSet)))  // if at least one valid atom was added, write config
-        FragmentList->ListOfMolecules[NumberOfFragments++] = StoreFragmentFromKeySet(out, CurrentSet, IsAngstroem);
+    for(Graph::iterator runner = FragmentList->begin();runner != FragmentList->end(); runner++) {
+      if (!InputFile.eof()) {
+        InputFile >> TEFactor;
+        (*runner).second.second = TEFactor;
+        *out << Verbose(2) << "Setting " << ++NumberOfFragments << " fragment's TEFactor to " << (*runner).second.second << "." << endl; 
+      } else {
+        status = false;
+        break;
+      }
     }
     // 4. Free and done
-    Free((void **)&buffer, "molecule::ParseKeySetFile - *buffer");
-    cout << "done." << endl;
+    InputFile.close();
+    *out << Verbose(1) << "done." << endl;
   } else {
-    cout << "File not found." << endl;
+    *out << Verbose(1) << "File " << filename << " not found." << endl;
     status = false;
   }
+
   Free((void **)&filename, "molecule::ParseKeySetFile - filename");
+
+  return status;
+};
+
+/** Stores keysets and TEFactors to file.
+ * \param *out output stream for debugging
+ * \param KeySetList Graph with Keysets and factors
+ * \param *path path to file
+ * \return true - file written successfully, false - writing failed
+ */
+bool molecule::StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path)
+{
+  ofstream output;
+  bool status =  true;
+  string line;
+  string::iterator ende;
+
+  // open KeySet file
+  line = path;
+  line.append("/");
+  line += FRAGMENTPREFIX;
+  ende = line.end();
+  line += KEYSETFILE;
+  output.open(line.c_str(), ios::out);
+  *out << Verbose(1) << "Saving key sets of the total graph ... ";
+  if(output != NULL) {
+    for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++) {
+      for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++)
+        output << *sprinter << "\t";
+      output << endl;
+    }
+    *out << "done." << endl;
+  } else {
+    cerr << "Unable to open " << line << " for writing keysets!" << endl;
+    status = false;
+  }
+  output.close();
+  output.clear();
+
+  // open TEFactors file
+  line.erase(ende, line.end());
+  line += TEFACTORSFILE;
+  output.open(line.c_str(), ios::out);
+  *out << Verbose(1) << "Saving TEFactors of the total graph ... ";
+  if(output != NULL) {
+    for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++)
+      output << (*runner).second.second << endl;
+    *out << Verbose(1) << "done." << endl;
+  } else {
+    *out << Verbose(1) << "failed to open file" << line << "." << endl;
+    status = false;
+  }
+  output.close();
 
   return status;
@@ -1946,10 +2022,10 @@
   ofstream AdjacencyFile;
   atom *Walker = NULL;
-  char *filename = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::StoreAdjacencyToFile - filename");
+  stringstream line;
   bool status = true;
 
-  sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, ADJACENCYFILE);
-  AdjacencyFile.open(filename);
-  cout << Verbose(1) << "Saving adjacency list ... ";
+  line << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
+  AdjacencyFile.open(line.str().c_str(), ios::out);
+  *out << Verbose(1) << "Saving adjacency list ... ";
   if (AdjacencyFile != NULL) {
     Walker = start;
@@ -1962,10 +2038,9 @@
     }
     AdjacencyFile.close();
-    cout << "done." << endl;
+    *out << Verbose(1) << "done." << endl;
   } else {
-    cout << "failed." << endl;
+    *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
     status = false;
   }
-  Free((void **)&filename, "molecule::StoreAdjacencyToFile - filename");
   
   return status;
@@ -1980,11 +2055,12 @@
 bool molecule::CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms)
 {
-  char *filename = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule - filename");
   ifstream File;
+  stringstream line;
   bool status = true;
-  
-  sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, ADJACENCYFILE);
-  File.open(filename);
-  *out << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ...";
+  char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
+  
+  line << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
+  File.open(line.str().c_str(), ios::out);
+  *out << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... " << endl;
   if (File != NULL) {
     // allocate storage structure
@@ -1995,7 +2071,7 @@
     // Parse the file line by line and count the bonds
     while (!File.eof()) {
-      File.getline(filename, MAXSTRINGSIZE);
+      File.getline(buffer, MAXSTRINGSIZE);
       stringstream line;
-      line.str(filename);
+      line.str(buffer);
       int AtomNr = -1;
       line >> AtomNr;
@@ -2031,14 +2107,14 @@
     File.clear();
     if (status) { // if equal we parse the KeySetFile
-      *out << " done: Equal." << endl;
+      *out << Verbose(1) << "done: Equal." << endl;
       status = true;
     } else
-      *out << " done: Not equal by " << NonMatchNumber << " atoms." << endl;
+      *out << Verbose(1) << "done: Not equal by " << NonMatchNumber << " atoms." << endl;
     Free((void **)&CurrentBonds, "molecule::CheckAdjacencyFileAgainstMolecule - **CurrentBonds");
   } else {
-    *out << " Adjacency file not found." << endl;
+    *out << Verbose(1) << "Adjacency file not found." << endl;
     status = false;
   }
-  Free((void **)&filename, "molecule::CheckAdjacencyFileAgainstMolecule - filename");
+  Free((void **)&buffer, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
   
   return status;
@@ -2055,4 +2131,7 @@
  * -# combines the generated molecule lists from all subgraphs
  * -# saves to disk: fragment configs, adjacency, orderatsite, keyset files
+ * Note that as we split "this" molecule up into a list of subgraphs, i.e. a MoleculeListClass, we have two sets
+ * of vertex indices: Global always means the index in "this" molecule, whereas local refers to the molecule or
+ * subgraph in the MoleculeListClass.
  * \param *out output stream for debugging
  * \param Order up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
@@ -2061,12 +2140,12 @@
 void molecule::FragmentMolecule(ofstream *out, int Order, config *configuration)
 {
-	MoleculeListClass **BondFragments = NULL;
-  MoleculeListClass *FragmentList = NULL;
+	MoleculeListClass *BondFragments = NULL;
   atom *Walker = NULL;
+  atom *OtherWalker = NULL;
+  bond *Binder = NULL;
   int *SortIndex = NULL;
   element *runner = NULL;
   int AtomNo;
   int MinimumRingSize;
-  int TotalFragmentCounter;
   int FragmentCounter;
   MoleculeLeafClass *MolecularWalker = NULL;
@@ -2074,5 +2153,12 @@
   fstream File;
   bool FragmentationToDo = true;
-
+  Graph **FragmentList = NULL;
+  Graph *TempFragmentList = NULL;
+  Graph TotalGraph;     // graph with all keysets however local numbers
+  KeySet *TempSet = NULL;
+  int TotalNumberOfKeySets = 0;
+  int KeySetCounter = 0;
+  atom ***ListOfLocalAtoms = NULL;
+  
   *out << endl;
 #ifdef ADDHYDROGEN
@@ -2082,8 +2168,10 @@
 #endif
 
+  // ===== 1. Check whether bond structure is same as stored in files ====
+  
   // fill the adjacency list
   CreateListOfBondsPerAtom(out);
 
-  // === compare it with adjacency file ===
+  // create lookup table for Atom::nr
   atom **ListOfAtoms = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::FragmentMolecule - **ListOfAtoms");
   Walker = start;
@@ -2099,15 +2187,72 @@
     FragmentationToDo = false;
   }
+  // === compare it with adjacency file ===
   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());
+  Free((void **)&ListOfAtoms, "molecule::FragmentMolecule - **ListOfAtoms");
+
+  // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
+  Subgraphs = DepthFirstSearchAnalysis((ofstream *)&*out, false, MinimumRingSize);
+  
+  // fill index lookup list for each subgraph from global nr to local pointer
+  ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **)*FragmentCounter, "molecule::FragmentMolecule - **ListOfLocalAtoms");
+  for (int i=0;i<FragmentCounter;i++) {
+    ListOfLocalAtoms[i] = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::FragmentMolecule - *ListOfLocalAtoms[]");
+    for(int j=0;j<AtomCount;j++)
+      ListOfLocalAtoms[i][j] = NULL;
+  }
+  FragmentCounter = 0;
+  MolecularWalker = Subgraphs;
+  while (MolecularWalker->next != NULL) {
+    MolecularWalker = MolecularWalker->next;
+    Walker = MolecularWalker->Leaf->start;
+    while (Walker->next != MolecularWalker->Leaf->end) {
+      Walker = Walker->next;
+#ifdef ADDHYDROGEN
+      if (Walker->type->Z != 1) // skip hydrogen
+#endif
+        ListOfLocalAtoms[FragmentCounter][Walker->GetTrueFather()->nr] = Walker;  // set present global id index to the local pointer
+    }
+    FragmentCounter++;
+  }
+
+  // fill the bond structure of the individually stored subgraphs
+  FragmentCounter = 0;
+  MolecularWalker = Subgraphs;
+  while (MolecularWalker->next != NULL) {
+    MolecularWalker = MolecularWalker->next;
+    *out << Verbose(1) << "Creating adjacency list for subgraph " << MolecularWalker << "." << endl;
+    Walker = MolecularWalker->Leaf->start;
+    while (Walker->next != MolecularWalker->Leaf->end) {
+      Walker = Walker->next;
+      AtomNo = Walker->father->nr;  // global id of the current walker
+      for(int i=0;i<NumberOfBondsPerAtom[AtomNo];i++) { // go through father's bonds and copy them all
+        Binder = ListOfBondsPerAtom[AtomNo][i];
+        OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->father)->nr];    // local copy of current bond partner of walker
+        if ((OtherWalker != NULL) && (OtherWalker->nr > Walker->nr))
+          MolecularWalker->Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);
+      }
+    }
+    //MolecularWalker->Leaf->CreateAdjacencyList(out, BondDistance);
+    MolecularWalker->Leaf->CreateListOfBondsPerAtom(out);
+    FragmentCounter++;
+  }
+
+  // ===== 3. if structure still valid, parse key set file and others =====
+  if (FragmentationToDo) {
+    // parse key sets into new graph
+    TempFragmentList = new Graph;
+    FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, TempFragmentList, configuration->GetIsAngstroem());
+  } else
+    *out << Verbose(1) << "Creation of lookup table Atom::nr <-> Atom failed!" << endl;
   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 (FragmentationToDo) { // if we parsed Adjacancy, check whether OrderAtSite is above Order everywhere
+  else
+     *out << Verbose(1) << "Parsing keyset file failed" << endl;
+  
+  // ===== 4. check globally whether there's something to do actually (first adaptivity check)
+  if (FragmentationToDo) { // check whether OrderAtSite is above Order everywhere
+    FragmentationToDo = false;
     Walker = start;
-    while (Walker->next != end) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron 
+    while (Walker->next != end) {
       Walker = Walker->next;
 #ifdef ADDHYDROGEN
@@ -2115,24 +2260,34 @@
 #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 *)&*out, false, MinimumRingSize);
+          FragmentationToDo = true;
+    }
+  } else
+    *out << Verbose(1) << "Parsing order at site file failed" << endl;
+  
+  if (!FragmentationToDo)
+    *out << Verbose(0) << "Order at every site is already equal or above desired order " << Order << "." << endl;
+
+  // =================================== Begin of FRAGMENTATION =============================== 
+  // check cyclic lengths
+  if ((MinimumRingSize != -1) && (Order >= MinimumRingSize)) {
+    *out << Verbose(0) << "Bond order " << Order << " greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
+  } else {
+    FragmentCounter = 0;
     MolecularWalker = Subgraphs;
-    // fill the bond structure of the individually stored subgraphs
+    // count subgraphs and allocate fragments
     while (MolecularWalker->next != NULL) {
       MolecularWalker = MolecularWalker->next;
-      *out << Verbose(1) << "Creating adjacency list for subgraph " << MolecularWalker << "." << endl;
-      MolecularWalker->Leaf->CreateAdjacencyList(out, BondDistance);
-      MolecularWalker->Leaf->CreateListOfBondsPerAtom(out);
-    }
+      FragmentCounter++;
+    }
+    FragmentList = (Graph **) Malloc(sizeof(Graph *)*FragmentCounter, "molecule::FragmentMolecule - **BondFragments");
     
-    // === fragment all subgraphs ===
-    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 {
+    // ===== 6a. fill RootStack for each subgraph (second adaptivity check) ===== 
+    // NOTE: (keep this extern of following while loop, as lateron we may here look for which site to add to which subgraph)
+    // fill the root stack
+    KeyStack RootStack[FragmentCounter];
+
+    if (Order == -1) {  // means we want to increase order adaptively
+      cerr << "Adaptive chosing of sites not yet implemented!" << endl;
+    } else {  // means we just want to increase it at every site by one
       FragmentCounter = 0;
       MolecularWalker = Subgraphs;
@@ -2140,19 +2295,6 @@
       while (MolecularWalker->next != NULL) {
         MolecularWalker = MolecularWalker->next;
-        FragmentCounter++;
-      }
-      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;
+        // find first root candidates  
         RootStack[FragmentCounter].clear();  
-        // find first root candidates  
         Walker = MolecularWalker->Leaf->start;
         while (Walker->next != MolecularWalker->Leaf->end) { // go through all (non-hydrogen) atoms
@@ -2166,76 +2308,105 @@
         FragmentCounter++;
       }
-      
-      // fill the bond fragment list
+    }
+
+   // ===== 6b. assign each keyset to its respective subgraph ===== 
+    if ((TempFragmentList != NULL) && (TempFragmentList->size() != 0)) { // if there are some scanned keysets at all
+      // spread the keysets
       FragmentCounter = 0;
-      TotalFragmentCounter = 0;
       MolecularWalker = Subgraphs;
       while (MolecularWalker->next != NULL) {
         MolecularWalker = MolecularWalker->next;
-        *out << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
-        if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
-            // output ListOfBondsPerAtom for debugging
-          *out << Verbose(0) << endl << "From Contents of ListOfBondsPerAtom, all non-hydrogen atoms:" << endl;
-          Walker = MolecularWalker->Leaf->start;
-          while (Walker->next != MolecularWalker->Leaf->end) {
-            Walker = Walker->next; 
-  #ifdef ADDHYDROGEN
-            if (Walker->type->Z != 1) {   // regard only non-hydrogen
-  #endif
-              *out << Verbose(0) << "Atom " << Walker->Name << " has Bonds: "<<endl;
-              for(int j=0;j<MolecularWalker->Leaf->NumberOfBondsPerAtom[Walker->nr];j++) {
-                *out << Verbose(1) << *(MolecularWalker->Leaf->ListOfBondsPerAtom)[Walker->nr][j] << endl;
-              }
-  #ifdef ADDHYDROGEN
+        // assign scanned keysets
+        FragmentList[FragmentCounter] = new Graph;
+        TempSet = new KeySet;
+        KeySetCounter = 0;
+        for(Graph::iterator runner = TempFragmentList->begin();runner != TempFragmentList->end(); runner++) { // key sets contain global numbers!
+          if ( ListOfLocalAtoms[FragmentCounter][FindAtom(*((*runner).first.begin()))->nr]->nr != -1) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
+            // translate keyset to local numbers
+            for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
+              TempSet->insert(ListOfLocalAtoms[FragmentCounter][FindAtom(*sprinter)->nr]->nr);
+            // insert into FragmentList
+            FragmentList[FragmentCounter]->insert(GraphPair (*TempSet, pair<int,double>(KeySetCounter++, (*runner).second.second)));
+          }
+          TempSet->clear();
+        }
+        delete(TempSet);
+        if (KeySetCounter == 0)
+          delete(FragmentList[FragmentCounter]);
+        else
+          *out << Verbose(1) << KeySetCounter << " atoms were put into subgraph " << FragmentCounter << "." << endl;
+        FragmentCounter++;
+      }
+    } else  // otherwise make sure all lists are initialised to NULL
+      for(int i=0;i<FragmentCounter;i++)
+        FragmentList[i] = NULL;
+    
+
+   // ===== 7. fill the bond fragment list =====
+    FragmentCounter = 0;
+    TotalNumberOfKeySets = 0;
+    MolecularWalker = Subgraphs;
+    while (MolecularWalker->next != NULL) {
+      MolecularWalker = MolecularWalker->next;
+      *out << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
+      if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
+          // output ListOfBondsPerAtom for debugging
+        *out << Verbose(0) << endl << "From Contents of ListOfBondsPerAtom, all non-hydrogen atoms:" << endl;
+        Walker = MolecularWalker->Leaf->start;
+        while (Walker->next != MolecularWalker->Leaf->end) {
+          Walker = Walker->next; 
+#ifdef ADDHYDROGEN
+          if (Walker->type->Z != 1) {   // regard only non-hydrogen
+#endif
+            *out << Verbose(0) << "Atom " << Walker->Name << " has Bonds: "<<endl;
+            for(int j=0;j<MolecularWalker->Leaf->NumberOfBondsPerAtom[Walker->nr];j++) {
+              *out << Verbose(1) << *(MolecularWalker->Leaf->ListOfBondsPerAtom)[Walker->nr][j] << endl;
             }
-  #endif
+#ifdef ADDHYDROGEN
           }
-          *out << endl;
+#endif
+        }
+        *out << endl;
+      
+        *out << Verbose(0) << endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= " << endl;
+        *out << Verbose(0) << "Begin of bond fragmentation." << endl;
         
-          *out << Verbose(0) << endl << " ========== BOND ENERGY ========================= " << endl;
-          *out << Verbose(0) << "Begin of bond fragmentation." << endl;
-          BondFragments[FragmentCounter] = NULL;
-          // 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++;
-          }
-          *out << k << "/" << BondFragments[FragmentCounter]->NumberOfMolecules << " fragments generated from the keysets." << endl;
-        } else {
-          *out << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
+        // call BOSSANOVA method
+        MolecularWalker->Leaf->FragmentBOSSANOVA(out, FragmentList[FragmentCounter], RootStack[FragmentCounter]);
+        
+        // translate list into global numbers (i.e. valid in "this" molecule, not in MolecularWalker->Leaf)
+        KeySet *TempSet = new KeySet;
+        for(Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) {
+          for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
+            TempSet->insert((MolecularWalker->Leaf->FindAtom(*sprinter))->GetTrueFather()->nr);
+          TotalGraph.insert(GraphPair(*TempSet, pair<int,double>(TotalNumberOfKeySets++, (*runner).second.second)));
+          TempSet->clear();
         }
-        TotalFragmentCounter += BondFragments[FragmentCounter]->NumberOfMolecules;
-        FragmentCounter++;  // next fragment list
+        delete(TempSet);
+        delete(FragmentList[FragmentCounter]);
+      } else {
+        cerr << "Subgraph " << MolecularWalker << " has no atoms!" << endl;
       }
-    }
-  
-    // === combine bond fragments list into a single one ===
-    FragmentList = new MoleculeListClass(TotalFragmentCounter, AtomCount);
-    TotalFragmentCounter = 0;
-    for (int i=0;i<FragmentCounter;i++) {
-      for(int j=0;j<BondFragments[i]->NumberOfMolecules;j++) {
-        FragmentList->ListOfMolecules[TotalFragmentCounter] =  BondFragments[i]->ListOfMolecules[j];
-        BondFragments[i]->ListOfMolecules[j] = NULL;
-        TotalFragmentCounter++;
-      }
-      delete(BondFragments[i]);
-    }
-    Free((void **)&BondFragments, "molecule::FragmentMolecule - **BondFragments");
-  } else
-    *out << Verbose(0) << "Nothing to do, using only fragments reconstructed from the KeySetFile." << endl;
-  // ==================================== End of FRAGMENTATION ================================
-  
-  // === Save fragments' configuration to disk ===
-  if (FragmentList != NULL) {
+      FragmentCounter++;  // next fragment list
+    }
+  }
+  Free((void **)&FragmentList, "molecule::FragmentMolecule - **FragmentList");
+
+  // ===== 8. gather keyset lists (graphs) from all subgraphs and transform into MoleculeListClass =====
+  // allocate memory for the pointer array and transmorph graphs into full molecular fragments
+  BondFragments = new MoleculeListClass(TotalGraph.size(), AtomCount);
+  int k=0;
+  for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
+    KeySet test = (*runner).first;
+    *out << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl;
+    BondFragments->ListOfMolecules[k] = StoreFragmentFromKeySet(out, test, configuration);
+    k++;
+  }
+  *out << k << "/" << BondFragments->NumberOfMolecules << " fragments generated from the keysets." << endl;
+
+ // ==================================== End of FRAGMENTATION ================================
+  
+  // ===== 10. Save fragments' configuration and keyset files et al to disk ===
+  if (BondFragments->NumberOfMolecules != 0) {
     // create a SortIndex to map from BFS labels to the sequence in which the atoms are given in the config file
     SortIndex = (int *) Malloc(sizeof(int)*AtomCount, "molecule::FragmentMolecule: *SortIndex");
@@ -2255,15 +2426,20 @@
       }
     }
-    *out << Verbose(1) << "Writing " << FragmentList->NumberOfMolecules << " possible bond fragmentation configs" << endl;
-    if (FragmentList->OutputConfigForListOfFragments(out, configuration, SortIndex))
+    *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
+    if (BondFragments->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"); 
+      *out << Verbose(1) << "Some config writing failed." << endl;
+
+    // store force index reference file
+    BondFragments->StoreForcesFile(out, configuration->configpath, SortIndex);
     
-    // === store Adjacency file ===
+    // store keysets file 
+    StoreKeySetFile(out, TotalGraph, configuration->configpath);
+
+    // store Adjacency file 
     StoreAdjacencyToFile(out, configuration->configpath);
 
-    // Store adaptive orders into file
+    // store adaptive orders into file
     StoreOrderAtSiteFile(out, configuration->configpath);
     
@@ -2275,8 +2451,14 @@
     delete(FragmentList); // remove bond molecule from memory
     FragmentList = NULL;
+    Free((void **)&SortIndex, "molecule::FragmentMolecule: *SortIndex"); 
   } else
     *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
   
+  // free the index lookup list
+  for (int i=0;i<FragmentCounter;i++)
+    Free((void **)&ListOfLocalAtoms[i], "molecule::FragmentMolecule - *ListOfLocalAtoms[]");
+  Free((void **)&ListOfLocalAtoms, "molecule::FragmentMolecule - **ListOfLocalAtoms");
   // free subgraph memory again
+  delete(TempFragmentList);
   if (Subgraphs != NULL) {
     while (Subgraphs->next != NULL) {
@@ -2303,5 +2485,5 @@
   line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
   file.open(line.str().c_str());
-  *out << Verbose(0) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl;
+  *out << Verbose(1) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl;
   if (file != NULL) {
     atom *Walker = start;
@@ -2312,6 +2494,8 @@
     }
     file.close();
+    *out << Verbose(1) << "done." << endl;
     return true;
   } else {
+    *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
     return false;
   }
@@ -2355,6 +2539,8 @@
     }
     file.close();
+    *out << Verbose(1) << "done." << endl;
     status = true;
   } else {
+    *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
     status = false;
   }
@@ -2938,4 +3124,5 @@
   int CurrentIndex;
   int *Labels;
+  double TEFactor;
   int *ShortestPathList;
   bool **UsedList;
@@ -3083,12 +3270,11 @@
  * with SP of 2, then those with SP of 3, then those with SP of 4 and so on.
  * \param *out output stream for debugging
- * \param Order bond order (limits BFS exploration and "number of digits" in power set generation 
- * \param *ReturnKeySets Graph structure to insert found keysets/fragments into
+ * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
+ * \param FragmentSearch UniqueFragments structure containing TEFactor, root atom and so on 
  * \param RestrictedKeySet Restricted vertex set to use in context of molecule
- * \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::PowerSetGenerator(ofstream *out, int Order, Graph *ReturnKeySets, KeySet RestrictedKeySet, int RootKeyNr) 
+int molecule::PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet) 
 {
   int SP, UniqueIndex, AtomKeyNr;
@@ -3100,29 +3286,11 @@
   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::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.Leaflet = ReturnKeySets;      // set to insertion graph
-  FragmentSearch.Root = FindAtom(RootKeyNr);
+  int RootKeyNr = FragmentSearch.Root->nr;
+  int Counter = FragmentSearch.FragmentCounter;
+
   for (int i=0;i<AtomCount;i++) {
-    FragmentSearch.Labels[i] = -1;
-    FragmentSearch.ShortestPathList[i] = -1;
     PredecessorList[i] = NULL;
   }
-  for (int i=0;i<Order;i++) {
-    FragmentSearch.BondsPerSPList[2*i] = new bond();    // start node
-    FragmentSearch.BondsPerSPList[2*i+1] = new bond();  // end node
-    FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1];     // intertwine these two 
-    FragmentSearch.BondsPerSPList[2*i+1]->previous = FragmentSearch.BondsPerSPList[2*i];
-    FragmentSearch.BondsPerSPCount[i] = 0;
-  }  
-  
+
   *out << endl;
   *out << Verbose(0) << "Begin of PowerSetGenerator with order " << Order << " at Root " << *FragmentSearch.Root << "." << endl;
@@ -3183,5 +3351,5 @@
           *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
+        if ((OtherWalker != PredecessorList[AtomKeyNr]) && (OtherWalker->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
           FragmentSearch.ShortestPathList[OtherWalker->nr] = SP+1;
           *out << Verbose(3) << "Set Shortest Path to " << FragmentSearch.ShortestPathList[OtherWalker->nr] << "." << endl;
@@ -3204,6 +3372,6 @@
             } 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;
+        } else *out << Verbose(3) << "Not passing on, as index of " << *OtherWalker << " " << OtherWalker->nr << " is smaller than that of Root " << RootKeyNr << " or this is my predecessor." << endl;
+      } else *out << Verbose(2) << "Is not in the restricted keyset or skipping hydrogen " << *OtherWalker << "." << endl;
     }
   }
@@ -3271,5 +3439,5 @@
   }
 
-/*  // as FragmentSearch structure is used only once, we don't have to clean it anymore
+  // 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;
@@ -3292,21 +3460,12 @@
     *out << "cleaned." << endl;
   }
-*/
+
   // free allocated memory
   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::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);
   
   // return list  
   *out << Verbose(0) << "End of PowerSetGenerator." << endl;
-  return FragmentSearch.FragmentCounter;
+  return (FragmentSearch.FragmentCounter - Counter);
 };
 
@@ -3459,5 +3618,5 @@
   GraphTestPair testGraphInsert;
 
-  testGraphInsert = Fragment->Leaflet->insert(GraphPair (*Fragment->FragmentSet,pair<int,double>(Fragment->FragmentCounter,1)));  // store fragment number and current factor
+  testGraphInsert = Fragment->Leaflet->insert(GraphPair (*Fragment->FragmentSet,pair<int,double>(Fragment->FragmentCounter,Fragment->TEFactor)));  // store fragment number and current factor
   if (testGraphInsert.second) {
     *out << Verbose(2) << "KeySet " << Fragment->FragmentCounter << " successfully inserted." << endl;
@@ -3465,5 +3624,5 @@
   } else {
     *out << Verbose(2) << "KeySet " << Fragment->FragmentCounter << " failed to insert, present fragment is " << ((*(testGraphInsert.first)).second).first << endl;
-    ((*(testGraphInsert.first)).second).second ++;  // increase the "created" counter
+    ((*(testGraphInsert.first)).second).second += Fragment->TEFactor;  // increase the "created" counter
     *out << Verbose(2) << "New factor is " << ((*(testGraphInsert.first)).second).second << "." << endl;
   }
@@ -3513,9 +3672,11 @@
  * as these copies are filtered out via use of the hash table (KeySet).
  * \param *out output stream for debugging
+ * \param Fragment&*List list of already present keystacks (adaptive scheme) or empty list
+ * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
  * \return pointer to Graph list
  */
-Graph * molecule::FragmentBOSSANOVA(ofstream *out, KeyStack &RootStack)
-{
-  Graph *FragmentList = NULL, ***FragmentLowerOrdersList = NULL;
+void molecule::FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack)
+{
+  Graph ***FragmentLowerOrdersList = NULL;
   int Order, NumLevels, NumMolecules, TotalNumMolecules = 0, *NumMoleculesOfOrder = NULL;
   int counter = 0;
@@ -3523,4 +3684,5 @@
   KeyStack FragmentRootStack;
   int RootKeyNr, RootNr;
+  struct UniqueFragments FragmentSearch;
   
   *out << Verbose(0) << "Begin of FragmentBOSSANOVA." << endl;
@@ -3530,4 +3692,17 @@
   NumMoleculesOfOrder = (int *) Malloc(sizeof(int)*UpgradeCount, "molecule::FragmentBOSSANOVA: *NumMoleculesOfOrder");
   FragmentLowerOrdersList = (Graph ***) Malloc(sizeof(Graph **)*UpgradeCount, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");
+
+  // initialise the fragments structure
+  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.Root = FindAtom(RootKeyNr);
+  for (int i=0;i<AtomCount;i++) {
+    FragmentSearch.Labels[i] = -1;
+    FragmentSearch.ShortestPathList[i] = -1;
+  }
 
   // Construct the complete KeySet which we need for topmost level only (but for all Roots)
@@ -3551,16 +3726,30 @@
     Walker->GetTrueFather()->AdaptiveOrder++;
     Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;
-    
+
+    // initialise Order-dependent entries of UniqueFragments structure
+    FragmentSearch.BondsPerSPList = (bond **) Malloc(sizeof(bond *)*Order*2, "molecule::PowerSetGenerator: ***BondsPerSPList");
+    FragmentSearch.BondsPerSPCount = (int *) Malloc(sizeof(int)*Order, "molecule::PowerSetGenerator: *BondsPerSPCount");
+    for (int i=0;i<Order;i++) {
+      FragmentSearch.BondsPerSPList[2*i] = new bond();    // start node
+      FragmentSearch.BondsPerSPList[2*i+1] = new bond();  // end node
+      FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1];     // intertwine these two 
+      FragmentSearch.BondsPerSPList[2*i+1]->previous = FragmentSearch.BondsPerSPList[2*i];
+      FragmentSearch.BondsPerSPCount[i] = 0;
+    }  
+
     // allocate memory for all lower level orders in this 1D-array of ptrs
-    NumLevels = 1 << (Order); // (int)pow(2,Order);
+    NumLevels = 1 << (Order-1); // (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 KeySets of Bond Order " << Order << " for " << *Walker << ", NumLevels is " << NumLevels << ", " << RootStack.size() << " Roots remaining." << endl;
+    *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", NumLevels is " << NumLevels << ", " << (RootStack.size()-RootNr-1) << " Roots remaining." << endl;
 
     // Create list of Graphs of current Bond Order (i.e. F_{ij})
     FragmentLowerOrdersList[RootNr][0] =  new Graph;
-    NumMoleculesOfOrder[RootNr] = PowerSetGenerator(out, Walker->AdaptiveOrder, FragmentLowerOrdersList[RootNr][0], CompleteMolecule, RootKeyNr);
+    FragmentSearch.TEFactor = 1.;
+    FragmentSearch.Leaflet = FragmentLowerOrdersList[RootNr][0];      // set to insertion graph
+    FragmentSearch.Root = Walker;
+    NumMoleculesOfOrder[RootNr] = PowerSetGenerator(out, Walker->AdaptiveOrder, FragmentSearch, CompleteMolecule);
     *out << Verbose(1) << "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
     NumMolecules = 0;
@@ -3575,8 +3764,8 @@
           Order--;
         *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
-        for (int SubOrder=Order;SubOrder>1;SubOrder--) {
-          int dest = source + (1 << (Walker->AdaptiveOrder-SubOrder));
+        for (int SubOrder=Order-1;SubOrder>0;SubOrder--) {
+          int dest = source + (1 << (Walker->AdaptiveOrder-(SubOrder+1)));
           *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
-          *out << Verbose(0) << "Current SubOrder is: " << SubOrder-1 << " with source " << source << " to destination " << dest << "." << endl;
+          *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
     
           // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
@@ -3586,10 +3775,12 @@
           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);
+              Graph TempFragmentList;
+              FragmentSearch.TEFactor = -(*runner).second.second;
+              FragmentSearch.Leaflet = &TempFragmentList;      // set to insertion graph
+              FragmentSearch.Root = FindAtom(*sprinter);
+              NumMoleculesOfOrder[RootNr] += PowerSetGenerator(out, SubOrder, FragmentSearch, (*runner).first);
               // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
               *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
-              InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], *FragmentList, &NumMolecules);
-              delete(FragmentList);
+              InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], TempFragmentList, &NumMolecules);
             }
           }
@@ -3601,11 +3792,25 @@
     //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
     TotalNumMolecules += NumMoleculesOfOrder[RootNr];
-    *out << Verbose(1) << "Number of resulting molecules for Order " << Walker->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
+    *out << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
     RootStack.push_back(RootKeyNr); // put back on stack
     RootNr++;
+
+    // free Order-dependent entries of UniqueFragments structure for next loop cycle
+    Free((void **)&FragmentSearch.BondsPerSPCount, "molecule::PowerSetGenerator: *BondsPerSPCount");
+    for (int i=0;i<Order;i++) {
+      delete(FragmentSearch.BondsPerSPList[2*i]);
+      delete(FragmentSearch.BondsPerSPList[2*i+1]);
+    }  
+    Free((void **)&FragmentSearch.BondsPerSPList, "molecule::PowerSetGenerator: ***BondsPerSPList");
   }
   *out << Verbose(0) << "==============================================================================================================" << endl;
   *out << Verbose(1) << "Total number of resulting molecules is: " << TotalNumMolecules << "." << endl;
   *out << Verbose(0) << "==============================================================================================================" << endl;
+
+  // cleanup FragmentSearch structure
+  Free((void **)&FragmentSearch.ShortestPathList, "molecule::PowerSetGenerator: *ShortestPathList");
+  Free((void **)&FragmentSearch.Labels, "molecule::PowerSetGenerator: *Labels");
+  delete(FragmentSearch.FragmentSet);
+  
   // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein) 
   // 5433222211111111
@@ -3614,8 +3819,14 @@
   // 21
   // 1
+
   // Subsequently, we combine all into a single list (FragmentList)
 
   *out << Verbose(0) << "Combining the lists of all orders per order and finally into a single one." << endl;
-  FragmentList = new Graph;
+  if (FragmentList == NULL) {
+    FragmentList = new Graph;
+    counter = 0;
+  } else {
+    counter = FragmentList->size();
+  }
   RootNr = 0;
   while (!RootStack.empty()) {
@@ -3635,8 +3846,7 @@
   
   *out << Verbose(0) << "End of FragmentBOSSANOVA." << endl;
-  return FragmentList;
-};
-
-/** Comparision function for GSL heapsort on distances in two molecules.
+};
+
+/** Comparison function for GSL heapsort on distances in two molecules.
  * \param *a
  * \param *b
Index: src/molecules.hpp
===================================================================
--- src/molecules.hpp	(revision bd65791c80eaf19f1799a7c9599f95606c4a4796)
+++ src/molecules.hpp	(revision 2459b1b8b885ddb84c9c4feab2f53d9d08d718c0)
@@ -525,10 +525,12 @@
   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 ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList, bool IsAngstroem);
+  bool StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path);
+  bool StoreForcesFile(ofstream *out, MoleculeListClass *BondFragments, char *path, int *SortIndex);
   bool ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet);
   void BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem);
   /// -# BOSSANOVA
-  Graph * FragmentBOSSANOVA(ofstream *out, KeyStack &RootStack);
-	int PowerSetGenerator(ofstream *out, int Order, Graph *ReturnKeySets, KeySet RestrictedKeySet, int RootKeyNr);
+  void FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack);
+	int PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet);
   bool BuildInducedSubgraph(ofstream *out, const molecule *Father);
   molecule * StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem);
@@ -563,4 +565,5 @@
 
   /// Output configs.
+  bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
   bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
   void Output(ofstream *out);
Index: src/parser.cpp
===================================================================
--- src/parser.cpp	(revision bd65791c80eaf19f1799a7c9599f95606c4a4796)
+++ src/parser.cpp	(revision 2459b1b8b885ddb84c9c4feab2f53d9d08d718c0)
@@ -43,5 +43,5 @@
   stringstream line;
 
-  //line << argv[1] << FRAGMENTPREFIX << TEFACTORSFILE;
+  line << argv[1] << FRAGMENTPREFIX << KEYSETFILE;
   return FilePresent(line.str().c_str());
 };
@@ -114,5 +114,5 @@
   }
   while (!input.eof()) {
-    input >> tmp;
+    input.getline(filename, 1023);
     MatrixCounter++;
   }
