Changeset c75363 for molecuilder/src/moleculelist.cpp

Timestamp:

May 2, 2008, 1:25:48 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

30fda2

Parents:

e936b3

Message:

HUGE REWRITE to allow for adaptive increase of the bond order, first working commit

Lots of code was thrown out:
-BottomUp, TopDown and GetAtomicFragments
-TEFactors are now used as "CreationCounters", i.e. the count how often a fragment as been created (ideal would be only once)
-ReduceToUniqueOnes and stuff all thrown out, since they are out-dated since use of hash table
Other changes:
-CreateListofUniqueFragments renamed to PowerSetGenerator
-PowerSetGenerator goes not over all reachable roots, but one given by calling function FragmentBOSSANOVA
-FragmentBOSSANOVA loops over all possible root sites and hands this over to PowerSetGenerator
-by the virtue of the hash table it is not important anymore whether created keysets are unique or not, as this is recognized in log(n). Hence, the label < label is not important anymore (and not applicable in an adaptive scheme with old, parsed keysets and unknown labels) (THIS IS HOWEVER NOT DONE YET!)
The setup is then as follows:

1. FragmentMolecule
   1. parses adjacency, keysets and orderatsite files
   2. performs DFS to find connected subgraphs (to leave this in was a design decision: might be useful later)
   3. a RootStack is created for every subgraph (here, later we implement the "update 10 sites with highest energy contribution", and that's why this consciously not done in the following loop)
   4. in a loop over all subgraphs

d1. calls FragmentBOSSANOVA with this RootStack and within the subgraph molecule structure
d2. creates molecule (fragment)s from the returned keysets (StoreFragmentFromKeySet)

5. combines the generated molecule lists from all subgraphs
6. saves to disk: fragment configs, adjacency, orderatsite, keyset files

2. FragmentBOSSANOVA
   1. constructs a complete keyset of the molecule
   2. In a loop over all possible roots from the given rootstack

b1. increases order of root site
b2. calls PowerSetGenerator with this order, the complete keyset and the rootkeynr
b3. for all consecutive lower levels PowerSetGenerator is called with the suborder, the higher order keyset as the restricted one and each site in the set as the root)
b4. these are merged into a fragment list of keysets

3. All fragment lists (for all orders, i.e. from all destination fields) are merged into one list for return

3. PowerSetGenerator
   1. initialises UniqueFragments structure
   2. fills edge list via BFS
   3. creates the fragment by calling recursive function SPFragmentGenerator with UniqueFragments structure, 0 as root distance, the edge set, its dimension and the current suborder
   4. Free'ing structure
4. SPFragmentGenerator (nothing much has changed here)
   1. loops over every possible combination (2^{dimension of edge set)}

a1. inserts current set, if there's still space left

a11. yes: calls SPFragmentGenerator with structure, created new edge list and size respective to root distance+1
a12. no: stores fragment into keyset list by calling InsertFragmentIntoGraph

a2. removes all items added into the snake stack (in UniqueFragments structure) added during level (root distance) and current set

File:

• molecuilder/src/moleculelist.cpp (modified) (6 diffs)

molecuilder/src/moleculelist.cpp

@@ -21 +21 @@
 MoleculeListClass::MoleculeListClass(int NumMolecules, int NumAtoms)
 {
-  TEList = (double*) Malloc(sizeof(double)*NumMolecules, "MoleculeListClass:MoleculeListClass: *TEList");
   ListOfMolecules = (molecule **) Malloc(sizeof(molecule *)*NumMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
-  for (int i=0;i<NumMolecules;i++) {
-    TEList[i] = 0.;
+  for (int i=0;i<NumMolecules;i++)
     ListOfMolecules[i] = NULL;
-  }
   NumberOfMolecules = NumMolecules;
   NumberOfTopAtoms = NumAtoms;
@@ -36 +33 @@
 MoleculeListClass::~MoleculeListClass()
 {
-  cout << Verbose(3) << this << ": Freeing ListOfMolcules and TEList." << endl;
+  cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
   for (int i=0;i<NumberOfMolecules;i++) {
     if (ListOfMolecules[i] != NULL) { // if NULL don't free
@@ -46 +43 @@
   cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
   Free((void **)&ListOfMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
-  cout << Verbose(4) << "Freeing TEList." << endl;
-  Free((void **)&TEList, "MoleculeListClass:MoleculeListClass: *TEList");
-};
-
+};
+
+/** Compare whether two molecules are equal.
+ * \param *a molecule one
+ * \param *n molecule two
+ * \return lexical value (-1, 0, +1)
+ */
 int MolCompare(const void *a, const void *b)
 {
@@ -122 +122 @@
 };
 
-/** Returns an allocated index list which fragment should remain and which are doubly appearing.
- * Again, we use linked cell however not with \a threshold which in most cases is too small and would
- * generate too many cells, but a too be specified \a celldistance which should e.g. be chosen in
- * the magnitude of the typical bond distance. In order to avoid so far the need of a vector list,
- * we simply use the molecule structure and the vector sitting in the atom structure.
- * \param *out output stream for debugging
- * \param **MapList empty but allocated(MoleculeListClass::NumberOfMolecules), containing on return mapping of for which equal molecule which atoms corresponds to which one (needed for ForceFactors)
- * \param threshold threshold value for molecule::IsEqualToWithinThreshold()
- * \param cell_size 6 entries specifying the unit cell vectors (matrix is symmetric)
- * \param celldistance needed for linked-cell technique to find possible equals
- * \return allocated index list with \a Num entries.
- */
-int * MoleculeListClass::GetMappingToUniqueFragments(ofstream *out, double threshold, double *cell_size, double celldistance)
-{
-  int j, UniqueIndex, HighestNumber;
-  //int NumberCells, divisor[NDIM], n[NDIM], N[NDIM], index, Index;
-  int Counter;
-  //molecule **CellList;
-  int *MapList = NULL;
-  size_t *SortMap = NULL;
-  //atom *Walker = NULL, *OtherWalker = NULL;
-  //molecule *LinkedCellContainer = new molecule(NULL);  // container for all the center of gravities stored as atoms
-  
-  /// Allocated MapList with range NumberOfMolecules. 
-  MapList = (int *) Malloc(sizeof(int)*NumberOfMolecules, "MoleculeListClass::GetMappingToUniqueFragments: *MapList");
-  SortMap = (size_t *) Malloc(sizeof(size_t)*NumberOfMolecules, "MoleculeListClass::GetMappingToUniqueFragments: *MapList");
-  for (j=0;j<NumberOfMolecules;j++) {
-    if (ListOfMolecules[j] ==  NULL)
-      cerr << "ERROR: Molecule " << j << " is NULL!" << endl;
-    //else 
-      //cerr << "ERROR: Molecule " << j << " is " << ListOfMolecules[j] << " with count " << ListOfMolecules[j]->AtomCount << "." << endl;
-    MapList[j] = -1;
-    SortMap[j] = 0;
-  }
-  *out << Verbose(1) << "Get mapping to unique fragments with " << NumberOfMolecules << " fragments total." << endl;
-  
-  // sort the list with heapsort according to sort function
-  gsl_heapsort_index(SortMap, ListOfMolecules, NumberOfMolecules, sizeof(molecule *), MolCompare);
-  // check next neighbours and remap
-  Counter = 0;
-  MapList[SortMap[0]] = Counter++;
-  for(int i=1;i<NumberOfMolecules;i++) {
-    if (MolCompare(&ListOfMolecules[SortMap[i]], &ListOfMolecules[SortMap[i-1]]) == 0)
-      MapList[SortMap[i]] = MapList[SortMap[i-1]];
-    else
-      MapList[SortMap[i]] = Counter++; 
-  } 
-  Free((void **)&SortMap, "MoleculeListClass::GetMappingToUniqueFragments: *SortMap");
-
-  *out << Verbose(2) << "Map: ";
-  for(int i=0;i<NumberOfMolecules;i++)
-    *out << MapList[i] << " ";
-  *out << endl;
-
-  // bring MapList indices into an ascending order
-  HighestNumber = Counter;
-  *out << Verbose(3) << "HighestNumber: " << HighestNumber << "." << endl;
-  SortMap = (size_t *) Malloc(sizeof(size_t)*NumberOfMolecules, "MoleculeListClass::GetMappingToUniqueFragments: *SortMap");
-  for(int i=0;i<NumberOfMolecules;i++)
-    SortMap[i] = HighestNumber;   // is the first number that will not appear anymore in list   
-  UniqueIndex = 0;
-  for(int i=0;i<NumberOfMolecules;i++)
-    if (MapList[i] != -1) {
-      if ((unsigned int)SortMap[MapList[i]] == (unsigned int)HighestNumber)
-        SortMap[MapList[i]] = UniqueIndex++;
-      MapList[i] = SortMap[MapList[i]];
-    } else {
-      MapList[i] = -1;
-    }
-  *out << Verbose(2) << "Ascending Map: ";
-  for(int i=0;i<NumberOfMolecules;i++)
-    *out << MapList[i] << " ";
-  *out << endl;
-  Free((void **)&SortMap, "MoleculeListClass::GetMappingToUniqueFragments: *SortMap");
-
-  /// Return the constructed list.
-  return MapList;
-};
-
-
-/** Takes a list of molecules and returns the list with doubles removed.
- * Checks by using molecule::IsEqualToWithinThreshold() whether fragments within the list
- * are actually doubles. If not, the pointer are copied to a new list, if so,
- * they are dropped. The new list is then copied back to the reallocated \a **FragmentList
- * and the new number of elements therein returned.
- * \param *out output stream for debugging
- * \param *Map mapping of which index goes on which (from molecule::GetMappingToUniqueFragments())
- * \return number of molecules in the new realloacted **FragmentList
- */
-int MoleculeListClass::ReduceFragmentToUniqueOnes(ofstream *out, int *Map)
-{
-  *out << Verbose(2) << "Begin of ReduceFragmentToUniqueOnes" << endl;
-  /// Allocate temporary lists of size \a Num.
-  molecule **NewMolList = (molecule **) Malloc(sizeof(molecule *)*NumberOfMolecules, "MoleculeListClass::ReduceFragmentToUniqueOnes: NewList");
-  double *NewTEList = (double *) Malloc(sizeof(double)*NumberOfMolecules, "MoleculeListClass::ReduceFragmentToUniqueOnes: *NewTEList");
-  for(int i=0;i<NumberOfMolecules;i++) {
-    NewTEList[i] = 0;
-    NewMolList[i] = NULL;
-  }
-  int ReducedNum = 0, j;
-  *out << Verbose(3) << "Reducing TE and Forces lists." << endl;
-  for(int i=0;i<NumberOfMolecules;i++) {                       /// Count new reduced number of molecules ...
-    //*out << Verbose(4) << i << "th TE Value " << " goes to " << Map[i] << ": " << TEList[i] << "." << endl;
-    NewTEList[Map[i]] += TEList[i];
-    if (Map[i] == i) {                           /// ... by checking if Map[i] != i, then Fragment[i] is double of Fragment[Map[i]].
-      *out << Verbose(4) << "Copying molecule " << i << " as it is unique so far ... " << ListOfMolecules[i]<< endl;
-      NewMolList[ReducedNum++] = ListOfMolecules[i];   /// ..., whilst copying each value to temporary list, ...
-    } else {  // else free molecule cause it appears doubly
-      *out << Verbose(4) << "Removing molecule " << i << "." << endl;
-      delete(ListOfMolecules[i]);
-      ListOfMolecules[i] = NULL;
-    }
-  }
-  /// Reallocate \a **FragmentList ...
-  *out << Verbose(3) << "Reallocating to Unique Fragments:" << endl;
-  if ((ReducedNum != 0) && (TEList != 0)) {
-    // factor list
-    *out << Verbose(4) << "Reallocating TEList [" << TEList << "] to " << ReducedNum << endl;
-    TEList = (double *) ReAlloc(TEList, sizeof(double)*ReducedNum, "MoleculeListClass::ReduceFragmentToUniqueOnes: *TEList");
-
-    *out << Verbose(4) << "Copying values to new list: " << endl;
-    j = 0;  // is ReducedNum index
-    for(int i=0;i<NumberOfMolecules;i++) {/// copy value from temporary to return list
-      if (Map[i] == i) {
-      //if (fabs(NewTEList[i]) > MYEPSILON) { // this is a good check also, if TEList[i] = 0 then fragment cancels and may be dropped!
-        TEList[j] = NewTEList[i];
-        *out << Verbose(5) << "TE [i" << i << "<->j" << j << "]:" << NewTEList[i] << "->" << TEList[j];
-        j++;
-      }
-      *out << endl;
-    }
-    if (j != ReducedNum)
-      *out << "Panic:  j " << j << " != ReducedNum " << ReducedNum << "." << endl;
-      
-
-    // molecule list
-    *out << Verbose(4) << "Reallocating ListOfMolecules ... " << endl;
-    ListOfMolecules = (molecule **) ReAlloc(ListOfMolecules, sizeof(molecule *)*ReducedNum, "MoleculeListClass::ReduceFragmentToUniqueOnes: ListOfMolecules");
-    NumberOfMolecules = ReducedNum;
-    /// ... and copy the reduced number back
-    *out << Verbose(5) << "Transfering to ListOfMolecules ... ";
-    for(int i=0;i<ReducedNum;i++) {
-      ListOfMolecules[i] = NewMolList[i];
-      *out << NewMolList[i] << " -> " << ListOfMolecules[i] << "\t";
-    }
-    *out << endl;
-
-
-  } else {
-    Free((void **)&ListOfMolecules, "MoleculeListClass::ReduceFragmentToUniqueOnes: ListOfMolecules");
-    NumberOfMolecules = ReducedNum;
-    *out << Verbose(3) << "ReducedNum is " << ReducedNum << ", Number of Molecules is " << NumberOfMolecules << "." << endl;
-  }
-  // free all the lists again
-  *out << Verbose(3) << "Freeing temporary lists." << endl;
-  Free((void **)&NewTEList, "MoleculeListClass::ReduceFragmentToUniqueOnes: *NewTEList");      /// free temporary list again
-  Free((void **)&NewMolList, "MoleculeListClass::ReduceFragmentToUniqueOnes: *NewMolList");
-  *out << Verbose(2) << "End of ReduceFragmentToUniqueOnes" << endl;
-  return(ReducedNum);
-};
-
 /** Simple output of the pointers in ListOfMolecules.
  * \param *out output stream
@@ -295 +134 @@
 
 /** Writes a config file for each molecule in the given \a **FragmentList.
- * Note that molecules with a TEList factor of 0 are not stored!
- * \param *prefix prefix for the file name
+ * \param *out output stream for debugging
  * \param *configuration standard configuration to attach atoms in fragment molecule to.
  * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
  * \return true - success (each file was written), false - something went wrong.
  */
-bool MoleculeListClass::OutputConfigForListOfFragments(char *prefix, config *configuration, int *SortIndex)
+bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex)
 {
   ofstream outputFragment;
@@ -310 +148 @@
   atom *Walker = NULL;
   vector BoxDimension;
-  char TEFilename[MAXSTRINGSIZE];
   char *FragmentNumber;
   int FragmentCounter = 0;
   ofstream output;
-  element *runner = NULL;
   
   // store the fragments as config and as xyz
   for(int i=0;i<NumberOfMolecules;i++) {
-    //if (TEList[i] != 0) {
-      strcpy(PathBackup, configuration->configpath);
-
-      // scan all atoms in the current molecule for their fathers and write each vertex index to forces file
-
-      // correct periodic
-      ListOfMolecules[i]->ScanForPeriodicCorrection((ofstream *)&cout);
-
-      // output xyz file
-      FragmentNumber = FixedDigitNumber(NumberOfMolecules, FragmentCounter++);
-      sprintf(FragmentName, "%s/%s%s.conf.xyz", PathBackup, prefix, FragmentNumber);
-      outputFragment.open(FragmentName, ios::out);
-      cout << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << " as XYZ ...";
-      if (intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment))
-        cout << " done." << endl;
-      else
-        cout << " failed." << endl;
-      result = result && intermediateResult;
-      outputFragment.close();
-      outputFragment.clear();
-  
-      cout << Verbose(2) << "Contained atoms: ";
-      Walker = ListOfMolecules[i]->start;
-      while (Walker->next != ListOfMolecules[i]->end) {
-        Walker = Walker->next;
-        cout << Walker->Name << " ";
-      }
-      cout << endl;
-      // prepare output of config file
-      sprintf(FragmentName, "%s/%s%s.conf", PathBackup, prefix, FragmentNumber);
-      outputFragment.open(FragmentName, ios::out);
-      strcpy(PathBackup, configuration->configpath);
-      sprintf(FragmentName, "%s/%s%s/", PathBackup, prefix, FragmentNumber);
-      
-      // center on edge
-      ListOfMolecules[i]->CenterEdge((ofstream *)&cout, &BoxDimension);
-      ListOfMolecules[i]->SetBoxDimension(&BoxDimension);  // update Box of atoms by boundary
-      int j = -1;
-      for (int k=0;k<3;k++) {
-        j += k+1;
-        BoxDimension.x[k] = 5.;
-        ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
-      }
-      ListOfMolecules[i]->Translate(&BoxDimension);
-
-      // also calculate necessary orbitals
-      ListOfMolecules[i]->CountElements();  // this is a bugfix, atoms should should actually be added correctly to this fragment 
-      ListOfMolecules[i]->CalculateOrbitals(*configuration);
-      
-      // change path in config
-      configuration->SetDefaultPath(FragmentName);
-      cout << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << " as config ...";
-      if (intermediateResult = configuration->Save(&outputFragment, ListOfMolecules[i]->elemente, ListOfMolecules[i]))
-        cout << " done." << endl;
-      else
-        cout << " failed." << endl;
-      // restore old config
-      configuration->SetDefaultPath(PathBackup);
-      
-      result = result && intermediateResult;
-      outputFragment.close();
-      outputFragment.clear();
-      Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
-    }
-
-  strcpy(PathBackup, configuration->configpath);
-  // open file for the total energy factor
-  sprintf(TEFilename, "%s/%s%s", PathBackup, prefix, TEFACTORSFILE);
-  output.open(TEFilename, ios::out);
-  // output TEList (later to file AtomicTEList.dat)
-  cout << Verbose(2) << "Saving " << prefix << " energy factors ...";
-  //cout << Verbose(1) << "Final ATEList: ";
-  //less output << prefix << "TE\t";
-  for(int i=0;i<NumberOfMolecules;i++) {
-    //cout << TEList[i] << "\t";
-    //if (TEList[i] != 0)
-      output << TEList[i] << "\t";
-  }
-  //cout << endl;
-  output << endl;
+    // save default path as it is changed for each fragment
+    strcpy(PathBackup, configuration->GetDefaultPath());
+
+    // correct periodic
+    ListOfMolecules[i]->ScanForPeriodicCorrection(out);
+
+    // output xyz file
+    FragmentNumber = FixedDigitNumber(NumberOfMolecules, FragmentCounter++);
+    sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    outputFragment.open(FragmentName, ios::out);
+    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << " as XYZ ...";
+    if (intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment))
+      *out << " done." << endl;
+    else
+      *out << " failed." << endl;
+    result = result && intermediateResult;
+    outputFragment.close();
+    outputFragment.clear();
+
+    *out << Verbose(2) << "Contained atoms: ";
+    Walker = ListOfMolecules[i]->start;
+    while (Walker->next != ListOfMolecules[i]->end) {
+      Walker = Walker->next;
+      *out << Walker->Name << " ";
+    }
+    *out << endl;
+    // prepare output of config file
+    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    outputFragment.open(FragmentName, ios::out);
+    strcpy(PathBackup, configuration->configpath);
+    sprintf(FragmentName, "%s/%s%s/", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    
+    // center on edge
+    ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
+    ListOfMolecules[i]->SetBoxDimension(&BoxDimension);  // update Box of atoms by boundary
+    int j = -1;
+    for (int k=0;k<3;k++) {
+      j += k+1;
+      BoxDimension.x[k] = 5.;
+      ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
+    }
+    ListOfMolecules[i]->Translate(&BoxDimension);
+
+    // also calculate necessary orbitals
+    ListOfMolecules[i]->CountElements();  // this is a bugfix, atoms should should actually be added correctly to this fragment 
+    ListOfMolecules[i]->CalculateOrbitals(*configuration);
+    
+    // change path in config
+    configuration->SetDefaultPath(FragmentName);
+    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << " as config ...";
+    if (intermediateResult = configuration->Save(&outputFragment, ListOfMolecules[i]->elemente, ListOfMolecules[i]))
+      *out << " done." << endl;
+    else
+      *out << " failed." << endl;
+    // restore old config
+    configuration->SetDefaultPath(PathBackup);
+    
+    result = result && intermediateResult;
+    outputFragment.close();
+    outputFragment.clear();
+    Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
+  }
+
+  // open KeySet file
+  sprintf(FragmentName, "%s/%s%s", configuration->configpath,  FRAGMENTPREFIX , KEYSETFILE);
+  output.open(FragmentName, ios::out);
+  *out << Verbose(2) << "Saving " << FRAGMENTPREFIX << " key sets of each subgraph ...";
+  for(int j=0;j<NumberOfMolecules;j++) {
+    Walker = ListOfMolecules[j]->start;
+    while(Walker->next != ListOfMolecules[j]->end) {
+      Walker = Walker->next;        
+#ifdef ADDHYDROGEN
+      if ((Walker->GetTrueFather() != NULL) && (Walker->type->Z != 1)) // if there is a real father, prints its index
+#else
+      if ((Walker->GetTrueFather() != NULL)) // if there is a real father, prints its index
+#endif 
+        output <<  Walker->GetTrueFather()->nr << "\t";
+#ifdef ADDHYDROGEN
+      else  // otherwise a -1 to indicate an added saturation hydrogen
+        output << "-1\t";
+#endif 
+    }
+    output << endl;
+  }
   output.close();
-  cout << " done." << endl;
-  
-  // open file for the force factors
-  sprintf(TEFilename, "%s/%s%s", PathBackup, prefix, FORCESFILE);
-  output.open(TEFilename, ios::out);
-  //cout << Verbose(1) << "Final AtomicForcesList: ";
-  cout << Verbose(2) << "Saving " << prefix << " force factors ...";
-  //output << prefix << "Forces" << endl;
-  for(int j=0;j<NumberOfMolecules;j++) {
-    //if (TEList[j] != 0) {
-    runner = ListOfMolecules[j]->elemente->start;
-    while (runner->next != ListOfMolecules[j]->elemente->end) { // go through every element
-      runner = runner->next;
-      if (ListOfMolecules[j]->ElementsInMolecule[runner->Z]) { // if this element got atoms
-        Walker = ListOfMolecules[j]->start;
-        while (Walker->next != ListOfMolecules[j]->end) { // go through every atom of this element
-          Walker = Walker->next;
-          if (Walker->type->Z == runner->Z) {
-            if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a real father, prints its index
-              //cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", its number " << Walker->GetTrueFather()->nr << " and index " << SortIndex[Walker->GetTrueFather()->nr] << "." << endl;  
-              output <<  SortIndex[Walker->GetTrueFather()->nr] << "\t";
-            } else  // otherwise a -1 to indicate an added saturation hydrogen
-              output << "-1\t";
-          }
-        }
-      }
-    }
-    //cout << endl << endl;
-    output << endl;
-  }
-  output.close();
-  cout << " done." << endl;
-
-  // open KeySet file
-  sprintf(TEFilename, "%s/%s%s", PathBackup, prefix, "KeySets.dat");
-  output.open(TEFilename, ios::out);
-  cout << Verbose(2) << "Saving " << prefix << " key sets of each subgraph ...";
-  for(int j=0;j<NumberOfMolecules;j++) {
-    //if (TEList[j] != 0) {
-      Walker = ListOfMolecules[j]->start;
-      while(Walker->next != ListOfMolecules[j]->end) {
-        Walker = Walker->next;        
-#ifdef ADDHYDROGEN
-        if ((Walker->GetTrueFather() != NULL) && (Walker->type->Z != 1)) // if there is a real father, prints its index
-#else
-        if ((Walker->GetTrueFather() != NULL)) // if there is a real father, prints its index
-#endif 
-          output <<  Walker->GetTrueFather()->nr << "\t";
-#ifdef ADDHYDROGEN
-        else  // otherwise a -1 to indicate an added saturation hydrogen
-          output << "-1\t";
-#endif 
-      }
-      //cout << endl << endl;
-      output << endl;
-    }
-  output.close();
-  cout << " done." << endl;
+  *out << " done." << endl;
   
   // printing final number
-  cout << "Final number of fragments: " << FragmentCounter << "." << endl;
+  *out << "Final number of fragments: " << FragmentCounter << "." << endl;
       
   return result;
 };
-
-/** Reduce the list of molecules to unique pieces.
- * molecule::IsEqualToWithinThreshold() is used by GetMappingToUniqueFragments() to provide a mapping for
- * ReduceFragmentToUniqueOnes(). 
- * \param *out output stream for debugging
- * \param cell_size 6 entries specifying the unit cell vectors (matrix is symmetric) for getting the map
- * \param celldistance needed for linked-cell technique when getting the map
- */
-void MoleculeListClass::ReduceToUniqueList(ofstream *out, double *cell_size, double celldistance)
-{
-  int *ReduceMapping = NULL;
-  
-  *out << Verbose(0) << "Begin of ReduceToUniqueList." << endl; 
-  // again, check whether there are equal fragments by ReduceToUniqueOnes
-  *out << Verbose(1) << "Get reduce mapping." << endl;
-  ReduceMapping = GetMappingToUniqueFragments(out, 0.05, cell_size, celldistance);
-  *out << Verbose(1) << "MapList: ";
-  for(int i=0;i<NumberOfMolecules;i++)
-    *out << ReduceMapping[i] << " ";
-  *out << endl << endl;
-  *out << Verbose(1) << "Apply reduce mapping." << endl;
-  ReduceFragmentToUniqueOnes(out, ReduceMapping);
-  Free((void **)&ReduceMapping, "MoleculeListClass::FragmentTopDown: *ReduceMapping");  // free map after reduce
-  *out << Verbose(1) << "Number of Fragments after Reducing is " << NumberOfMolecules << "." << endl;
-  *out << Verbose(0) << "End of ReduceToUniqueList." << endl; 
-};
-
-/** Fragments a molecule and resulting pieces recursively until the number of bonds
- * Here the idea is similar to the FragmentBottomUp() approach, only that we split up the molecule
- * bond per bond into left and right fragments and this recursively also with each yielded
- * fragment until there is no fragment with continuous number of bonds greater than \a BondDegree.
- * are below the given \a BondDegree.
- * \param *out output stream for debugging
- * \param BondDegree maximum number of continuous bonds in a fragment
- * \param bonddistance typical bond distance
- * \param *configuration configuration for writing config files for each fragment
- * \param CutCyclic cut cyclic bonds (saturate with hydrogen) or add 
- * \return pointer to MoleculeListClass with all the fragments or NULL if something went wrong.
- * \todo so far we volontarily mix up the above BondDegree definition and molecule::NoNonBonds
- * 
- * \bug  FragmentTopDown does not work right now due to NULL being given as cell_size to ReduceToUniqueOnes()
- */
-MoleculeListClass * MoleculeListClass::FragmentTopDown(ofstream *out, int BondDegree, double bonddistance, config *configuration, enum CutCyclicBond CutCyclic)
-{
-  int Num, No, Cyclic, NoBonds;
-  MoleculeListClass *ReturnList = NULL;
-  MoleculeListClass **FragmentsList = NULL;
-  molecule *CurrentMolecule = NULL;
-  bond *Binder = NULL;
- 
-  *out << Verbose(0) << "Begin of FragmentTopDown:" << this << "." << endl; 
-  Num = 0;
-  FragmentsList = (MoleculeListClass **) Malloc(sizeof(MoleculeListClass *)*NumberOfMolecules, "MoleculeListClass::FragmentTopDown: **FragmentsList");
-  // fragment each molecule into its own MoleculeListClass
-  for(int i=0;i<NumberOfMolecules;i++) {
-    CurrentMolecule = ListOfMolecules[i];
-    CurrentMolecule->CountAtoms(out);
-    Cyclic = CurrentMolecule->CountCyclicBonds(out);
-    NoBonds = 
-#ifdef ADDHYDROGEN
-CurrentMolecule->NoNonBonds
-#else 
-CurrentMolecule->BondCount
-#endif 
-                - ((CutCyclic == KeepBond) ? Cyclic : 0);
-    // check if NoNonBonds < BondDegree, then don't split any further: Here, correct for greatest continuous bond length!
-    *out << Verbose(0) << "Checking on Number of true Bonds " << NoBonds << " (i.e. no -H, if chosen no cyclic) greater than " << BondDegree << "." << endl;
-    if (NoBonds > BondDegree) {
-      //cerr << Verbose(1) << "TopDown Level of Molecule " << CurrentMolecule << ": " << (CurrentMolecule->NoNonBonds - BondDegree) << "." << endl; 
-      *out << Verbose(1) << "Breaking up molecule " << i << " in fragment list." << endl;
-      CurrentMolecule->CreateListOfBondsPerAtom(out);
-      // Get atomic fragments for the list
-      *out << Verbose(1) << "Getting Atomic fragments for molecule " << i << "." << endl;
-      MoleculeListClass *AtomicFragments = CurrentMolecule->GetAtomicFragments(out, NumberOfTopAtoms, configuration->GetIsAngstroem(), TEList[i]/(1. - NoBonds), CutCyclic);
-      
-      // build the atomic part of the list of molecules
-      *out << Verbose(1) << "Putting atomic fragment into list of molecules." << endl;
-      
-      // cutting cyclic bonds yields only a single fragment, not two as non-cyclic!
-      No = (CutCyclic == KeepBond) ? ((
-#ifdef ADDHYDROGEN
-CurrentMolecule->NoNonBonds - Cyclic)*2) : CurrentMolecule->NoNonBonds*2 - Cyclic; 
-#else 
-CurrentMolecule->BondCount - Cyclic)*2) : CurrentMolecule->BondCount*2 - Cyclic; 
-#endif 
-      
-      if (AtomicFragments != NULL) {
-        ReturnList = new MoleculeListClass(No+AtomicFragments->NumberOfMolecules, NumberOfTopAtoms);
-        for (int j=0;j<AtomicFragments->NumberOfMolecules;j++) { // shift all Atomic Fragments into this one here
-          ReturnList->ListOfMolecules[j+No] = AtomicFragments->ListOfMolecules[j];
-          AtomicFragments->ListOfMolecules[j] = NULL;
-          ReturnList->ListOfMolecules[j+No]->NoNonBonds = 0;
-          ReturnList->ListOfMolecules[j+No]->NoNonHydrogen = 1;
-          ReturnList->TEList[j+No] = AtomicFragments->TEList[j];
-        }
-        *out << Verbose(3) << "Freeing Atomic memory" << endl;
-        delete(AtomicFragments);
-        AtomicFragments = NULL;
-      } else {
-        *out << Verbose(2) << "AtomicFragments is NULL, filling list with whole molecule only." << endl; 
-        ReturnList = new MoleculeListClass(No, NumberOfTopAtoms);
-      }
-
-      // build the side pieces part of the list of molecules
-      *out << Verbose(1) << "Building side piece fragments and putting into list of molecules." << endl;
-      No = 0;
-      Binder = CurrentMolecule->first;
-      while(Binder->next != CurrentMolecule->last) {
-        Binder = Binder->next;
-#ifdef ADDHYDROGEN
-        if (Binder->HydrogenBond == 0)
-#endif        
-          if ((CutCyclic == SaturateBond) || (!Binder->Cyclic)) {
-            // split each bond into left and right side piece
-            if (Binder->Cyclic) {
-              molecule *dummy = NULL; // we lazily hand FragmentMoleculeByBond() a dummy as second fragment and ... 
-              *out << Verbose(1) << "Breaking up " << Binder->nr << "th and " << No << "th non-hydrogen, CYCLIC bond " << *Binder << " of molecule " << i << " [" << CurrentMolecule << "] into a single pieces." << endl; 
-              CurrentMolecule->FragmentMoleculeByBond(out, Binder, &(ReturnList->ListOfMolecules[No]), &(dummy), configuration->GetIsAngstroem(), CutCyclic);
-              delete(dummy);  // ... free the dummy which is due to the cyclic bond the exacy same fragment
-              *out << Verbose(1) << "Single fragment is " << ReturnList->ListOfMolecules[No] << "." << endl;
-      
-              // write down the necessary TE-summation in order to regain TE of whole molecule
-              *out << Verbose(2) << "Creating TEList for Fragment " << i << " of Bond " << No << "." << endl;
-#ifdef ADDHYDROGEN
-              ReturnList->TEList[No] = -TEList[i]/(1. -CurrentMolecule->NoNonBonds);
-#else              
-              ReturnList->TEList[No] = -TEList[i]/(1. -CurrentMolecule->BondCount);
-#endif
-            } else {
-              *out << Verbose(1) << "Breaking up " << Binder->nr << "th and " << No << "th non-hydrogen, non-cyclic bond " << *Binder << " of molecule " << i << " [" << CurrentMolecule << "] into left and right side pieces." << endl; 
-              CurrentMolecule->FragmentMoleculeByBond(out, Binder, &(ReturnList->ListOfMolecules[No]), &(ReturnList->ListOfMolecules[No+1]), configuration->GetIsAngstroem(), CutCyclic);
-              *out << Verbose(1) << "Left Fragment is " << ReturnList->ListOfMolecules[No] << ", right Fragment is " << ReturnList->ListOfMolecules[No+1] << "." << endl;
-      
-              // write down the necessary TE-summation in order to regain TE of whole molecule
-              *out << Verbose(2) << "Creating TEList for Fragment " << i << " of Bond " << No << "." << endl;
-              ReturnList->TEList[No] = -TEList[i]/(1. - NoBonds);
-              ReturnList->TEList[No+1] = -TEList[i]/(1. - NoBonds);
-            }
-            No += ((Binder->Cyclic) ? 1 : 2);
-          }
-      }
-
-      // Reduce this list
-      ReturnList->ReduceToUniqueList(out, NULL, bonddistance);
-
-      // recurse and receive List
-      *out << Verbose(1) << "Calling TopDown " << ReturnList<< " with filled list: ";
-      ReturnList->Output(out);
-      FragmentsList[i] = ReturnList->FragmentTopDown(out, BondDegree, bonddistance, configuration, CutCyclic);
-      *out << Verbose(1) << "Returning from TopDown " << ReturnList<< " with filled list: ";
-      ReturnList->Output(out);
-      
-      // remove list after we're done
-      *out << Verbose(2) << "Removing caller list." << endl;
-      delete(ReturnList);
-      ReturnList = NULL;
-      
-    } else {  // create a list with only a single molecule inside, transfer everything from "this" list to return list
-      *out << Verbose(1) << "Not breaking up molecule " << i << " in fragment list." << endl;
-      FragmentsList[i] = new MoleculeListClass(1, NumberOfTopAtoms);
-      FragmentsList[i]->ListOfMolecules[0] = ListOfMolecules[i];
-      ListOfMolecules[i] = NULL;
-      FragmentsList[i]->TEList[0] = TEList[i]; 
-    }
-    Num += FragmentsList[i]->NumberOfMolecules;
-  }
-  
-  // now, we have a list of MoleculeListClasses: combine all fragments list into one single list
-  *out << Verbose(2) << "Combining to a list of " << Num << " molecules and Memory cleanup of old list of lists." << endl; 
-  ReturnList = new MoleculeListClass(Num, NumberOfTopAtoms);
-  No = 0;
-  for(int i=0;i<NumberOfMolecules;i++) {
-    for(int j=0;j<FragmentsList[i]->NumberOfMolecules;j++) {
-      // transfer molecule
-      ReturnList->ListOfMolecules[No] = FragmentsList[i]->ListOfMolecules[j];
-      FragmentsList[i]->ListOfMolecules[j] = NULL;
-      // transfer TE factor
-      ReturnList->TEList[No] = FragmentsList[i]->TEList[j];
-      No++;
-    }
-    delete(FragmentsList[i]);
-    FragmentsList[i] = NULL;
-  }
-  Free((void **)&FragmentsList, "MoleculeListClass::FragmentTopDown: **FragmentsList");
-  
-  // Reduce the list to unique fragments
-  ReturnList->ReduceToUniqueList(out, NULL, bonddistance);
-  
-  // return pointer to list
-  *out << Verbose(0) << "Return with filled list: ";
-  ReturnList->Output(out);
-  
-  *out << Verbose(0) << "End of FragmentTopDown:" << this << "." << endl; 
-  return (ReturnList);
-}; 
 
 /******************************************* Class MoleculeLeafClass ************************************************/