/** \file MoleculeListClass.cpp
 * 
 * Function implementations for the class MoleculeListClass.
 * 
 */

#include "molecules.hpp"

/*********************************** Functions for class MoleculeListClass *************************/

/** Constructor for MoleculeListClass.
 */
MoleculeListClass::MoleculeListClass()
{
};

/** constructor for MoleculeListClass.
 * \param NumMolecules number of molecules to allocate for
 * \param NumAtoms number of atoms to allocate for
 */
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.;
    ListOfMolecules[i] = NULL;
  }
  NumberOfMolecules = NumMolecules;
  NumberOfTopAtoms = NumAtoms;
};


/** Destructor for MoleculeListClass.
 */
MoleculeListClass::~MoleculeListClass()
{
  cout << Verbose(3) << this << ": Freeing ListOfMolcules and TEList." << endl;
  for (int i=0;i<NumberOfMolecules;i++) {
    if (ListOfMolecules[i] != NULL) { // if NULL don't free
      cout << Verbose(4) << "ListOfMolecules: Freeing " << ListOfMolecules[i] << "." << endl;
      delete(ListOfMolecules[i]);
      ListOfMolecules[i] = NULL;
    }
  }
  cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
  Free((void **)&ListOfMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
  cout << Verbose(4) << "Freeing TEList." << endl;
  Free((void **)&TEList, "MoleculeListClass:MoleculeListClass: *TEList");
};

int MolCompare(const void *a, const void *b)
{
  int *aList = NULL, *bList = NULL;
  int Count, Counter, aCounter, bCounter;
  int flag;
  atom *aWalker = NULL;
  atom *bWalker = NULL;
  
  // sort each atom list and put the numbers into a list, then go through
  //cout << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl;
  if ( (**(molecule **)a).AtomCount < (**(molecule **)b).AtomCount ) {
    return -1;
  } else { if ((**(molecule **)a).AtomCount > (**(molecule **)b).AtomCount)
    return +1;
    else {
      Count = (**(molecule **)a).AtomCount;
      aList = (int *) Malloc(sizeof(int)*Count, "MolCompare: *aList");
      bList = (int *) Malloc(sizeof(int)*Count, "MolCompare: *bList");
  
      // fill the lists
      aWalker = (**(molecule **)a).start;
      bWalker = (**(molecule **)b).start;
      Counter = 0;
      aCounter = 0;
      bCounter = 0;
      while ((aWalker->next != (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
        aWalker = aWalker->next;
        bWalker = bWalker->next;
        if (aWalker->GetTrueFather() == NULL)
          aList[Counter] = Count + (aCounter++);
        else
          aList[Counter] = aWalker->GetTrueFather()->nr;
        if (bWalker->GetTrueFather() == NULL)
          bList[Counter] = Count + (bCounter++);
        else
          bList[Counter] = bWalker->GetTrueFather()->nr;
        Counter++;
      }
      // check if AtomCount was for real
      flag = 0;
      if ((aWalker->next == (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
        flag = -1;
      } else {
        if ((aWalker->next != (**(molecule **)a).end) && (bWalker->next == (**(molecule **)b).end))
          flag = 1;
      }
      if (flag == 0) {
        // sort the lists
        gsl_heapsort(aList,Count, sizeof(int), CompareDoubles); 
        gsl_heapsort(bList,Count, sizeof(int), CompareDoubles);
        // compare the lists 
        
        flag = 0;
        for(int i=0;i<Count;i++) {
          if (aList[i] < bList[i]) {
            flag = -1;
          } else {
            if (aList[i] > bList[i])
              flag = 1;
          }
          if (flag != 0)
            break;
        }
      }
      Free((void **)&aList, "MolCompare: *aList");
      Free((void **)&bList, "MolCompare: *bList");
      return flag;
    }
  }
  return  -1;
};

/** 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
 */
void MoleculeListClass::Output(ofstream *out)
{
  *out<< Verbose(1) << "MoleculeList: ";
  for (int i=0;i<NumberOfMolecules;i++)
    *out << ListOfMolecules[i] << "\t";
  *out << endl;
};

/** 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 *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)
{
  ofstream outputFragment;
  char FragmentName[255];
  char PathBackup[255];
  bool result = true;
  bool intermediateResult = true;
  atom *Walker = NULL;
  vector BoxDimension;
  char TEFilename[255];
  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->GetDefaultPath());

      // 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->GetDefaultPath());
      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->GetDefaultPath());
  // open file for the total energy factor
  sprintf(TEFilename, "%s/%s%s", PathBackup, prefix, "TE-Factors.dat");
  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;
  output.close();
  cout << " done." << endl;
  
  // open file for the force factors
  sprintf(TEFilename, "%s/%s%s", PathBackup, prefix, "Forces-Factors.dat");
  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
              cerr << "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;
  
  // printing final number
  cout << "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 ************************************************/

/** Constructor for MoleculeLeafClass root leaf.
 * \param *Up Leaf on upper level
 * \param *PreviousLeaf NULL - We are the first leaf on this level, otherwise points to previous in list
 */
//MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *Up = NULL, MoleculeLeafClass *Previous = NULL)
MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf = NULL)
{
//  if (Up != NULL)
//    if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?
//      Up->DownLeaf = this;
//  UpLeaf = Up;
//  DownLeaf = NULL;
  Leaf = NULL;
  previous = PreviousLeaf;
  if (previous != NULL) {
    MoleculeLeafClass *Walker = previous->next;
    previous->next = this;
    next = Walker;
  } else {
    next = NULL;
  }
};

/** Destructor for MoleculeLeafClass.
 */
MoleculeLeafClass::~MoleculeLeafClass()
{
//  if (DownLeaf != NULL) {// drop leaves further down
//    MoleculeLeafClass *Walker = DownLeaf;
//    MoleculeLeafClass *Next;
//    do {
//      Next = Walker->NextLeaf;
//      delete(Walker);
//      Walker = Next;
//    } while (Walker != NULL);
//    // Last Walker sets DownLeaf automatically to NULL
//  }
  // remove the leaf itself
  if (Leaf != NULL) {
    delete(Leaf);
    Leaf = NULL;
  }
  // remove this Leaf from level list
  if (previous != NULL)   
    previous->next = next;
//  } else { // we are first in list (connects to UpLeaf->DownLeaf)
//    if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))
//      NextLeaf->UpLeaf = UpLeaf;  // either null as we are top level or the upleaf of the first node
//    if (UpLeaf != NULL)
//      UpLeaf->DownLeaf = NextLeaf;  // either null as we are only leaf or NextLeaf if we are just the first
//  }
//  UpLeaf = NULL;
  if (next != NULL) // are we last in list
    next->previous = previous;
  next = NULL;
  previous = NULL;
};

/** Adds \a molecule leaf to the tree.
 * \param *ptr ptr to molecule to be added
 * \param *Previous previous MoleculeLeafClass referencing level and which on the level
 * \return true - success, false - something went wrong
 */
bool MoleculeLeafClass::AddLeaf(molecule *ptr, MoleculeLeafClass *Previous)
{
  return false;
};