Changeset 2459b1

Timestamp:

May 5, 2008, 4:19:05 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

3c4b08

Parents:

bd6579

Message:

ForcesFile and TEFactors are _needed_, reincorporated. UniqueFragments structure is now in BOSSANOVA

ForcesFile is again written, as we need it for the hydrogen not coming saturation (forces!)
TEFactors are back, as despite my notion they are needed in the evaluation
UniqueFragments structure is shifted from PowerSetGenerator to FragmentBOSSANOVA. Actually only for the labels - however, the if was changed to test the real indices (which are also always the same, which is better for adaptive runs) - but might be more useful there still.
analyzer and joiner again parse indices.

Location:

src

Files:

• analyzer.cpp (modified) (1 diff)
• defs.hpp (modified) (1 diff)
• joiner.cpp (modified) (1 diff)
• moleculelist.cpp (modified) (2 diffs)
• molecules.cpp (modified) (43 diffs)
• molecules.hpp (modified) (2 diffs)
• parser.cpp (modified) (2 diffs)

src/analyzer.cpp

@@ -58 +58 @@
 
   // ---------- 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 ---------------

src/defs.hpp

@@ -44 +44 @@
 #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"

src/joiner.cpp

@@ -46 +46 @@
 
   // ---------- 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;
 

src/moleculelist.cpp

@@ -133 +133 @@
 };
 
+
+/** 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 +270 @@
     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

src/molecules.cpp

@@ -1249 +1249 @@
   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 +1364 @@
                 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 +1376 @@
             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 +1382 @@
               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 +1389 @@
                     }
                   } 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 +1611 @@
   }
 
-  // further analysis of the found cycles (print rings, get minimum cycle length)
-  CyclicStructureAnalysis(out, MinimumRingSize);
-  
   // free all and exit
   delete(AtomStack);
@@ -1880 +1876 @@
 
 /** 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 +2022 @@
   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 +2038 @@
     }
     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 +2055 @@
 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 +2071 @@
     // 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 +2107 @@
     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 +2131 @@
  * -# 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 +2140 @@
 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 +2153 @@
   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 +2168 @@
 #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 +2187 @@
     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 +2260 @@
 #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 +2295 @@
       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 +2308 @@
         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 +2426 @@
       }
     }
-    *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 +2451 @@
     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 +2485 @@
   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 +2494 @@
     }
     file.close();
+    *out << Verbose(1) << "done." << endl;
     return true;
   } else {
+    *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
     return false;
   }
@@ -2355 +2539 @@
     }
     file.close();
+    *out << Verbose(1) << "done." << endl;
     status = true;
   } else {
+    *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
     status = false;
   }
@@ -2938 +3124 @@
   int CurrentIndex;
   int *Labels;
+  double TEFactor;
   int *ShortestPathList;
   bool **UsedList;
@@ -3083 +3270 @@
  * 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 +3286 @@
   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 +3351 @@
           *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 +3372 @@
             } 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 +3439 @@
   }
 
-/*  // 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 +3460 @@
     *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 +3618 @@
   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 +3624 @@
   } 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 +3672 @@
  * 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 +3684 @@
   KeyStack FragmentRootStack;
   int RootKeyNr, RootNr;
+  struct UniqueFragments FragmentSearch;
   
   *out << Verbose(0) << "Begin of FragmentBOSSANOVA." << endl;
@@ -3530 +3692 @@
   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 +3726 @@
     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 +3764 @@
           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 +3775 @@
           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 +3792 @@
     //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 +3819 @@
   // 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 +3846 @@
   
   *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

src/molecules.hpp

@@ -525 +525 @@
   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 +565 @@
 
   /// Output configs.
+  bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
   bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
   void Output(ofstream *out);

src/parser.cpp

@@ -43 +43 @@
   stringstream line;
 
-  //line << argv[1] << FRAGMENTPREFIX << TEFACTORSFILE;
+  line << argv[1] << FRAGMENTPREFIX << KEYSETFILE;
   return FilePresent(line.str().c_str());
 };
@@ -114 +114 @@
   }
   while (!input.eof()) {
-    input >> tmp;
+    input.getline(filename, 1023);
     MatrixCounter++;
   }