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

#include "molecules.hpp"

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

/** Constructor for MoleculeListClass.
 */
MoleculeListClass::MoleculeListClass()
{
  // empty lists
  ListOfMolecules.clear();
  MaxIndex = 1;
};

/** Destructor for MoleculeListClass.
 */
MoleculeListClass::~MoleculeListClass()
{
  cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
    cout << Verbose(4) << "ListOfMolecules: Freeing " << *ListRunner << "." << endl;
    delete (*ListRunner);
  }
  cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
  ListOfMolecules.clear(); // empty list
};

/** Insert a new molecule into the list and set its number.
 * \param *mol molecule to add to list.
 * \return true - add successful
 */
void MoleculeListClass::insert(molecule *mol)
{
  mol->IndexNr = MaxIndex++;
  ListOfMolecules.push_back(mol);
};

/** 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)
{
  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 = new int[Count];
      bList = new int[Count];

      // 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;
        }
      }
      delete[] (aList);
      delete[] (bList);
      return flag;
    }
  }
  return -1;
};

/** Output of a list of all molecules.
 * \param *out output stream
 */
void MoleculeListClass::Enumerate(ofstream *out)
{
  int i=1;
  element* Elemental = NULL;
  atom *Walker = NULL;
  int Counts[MAX_ELEMENTS];
  double size=0;
  Vector Origin;

  // header
  *out << "Index\tName\t\tAtoms\tFormula\tCenter\tSize" << endl;
  cout << Verbose(0) << "-----------------------------------------------" << endl;
  if (ListOfMolecules.size() == 0)
    *out << "\tNone" << endl;
  else {
    Origin.Zero();
    for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
      // reset element counts
      for (int j = 0; j<MAX_ELEMENTS;j++)
        Counts[j] = 0;
      // count atoms per element and determine size of bounding sphere
      size=0.;
      Walker = (*ListRunner)->start;
      while (Walker->next != (*ListRunner)->end) {
        Walker = Walker->next;
        Counts[Walker->type->Z]++;
        if (Walker->x.DistanceSquared(&Origin) > size)
          size = Walker->x.DistanceSquared(&Origin);
      }
      // output Index, Name, number of atoms, chemical formula
      *out << ((*ListRunner)->ActiveFlag ? "*" : " ") << (*ListRunner)->IndexNr << "\t" << (*ListRunner)->name << "\t\t" << (*ListRunner)->AtomCount << "\t";
      Elemental = (*ListRunner)->elemente->end;
      while(Elemental->previous != (*ListRunner)->elemente->start) {
        Elemental = Elemental->previous;
        if (Counts[Elemental->Z] != 0)
          *out << Elemental->symbol << Counts[Elemental->Z];
      }
      // Center and size
      *out << "\t" << (*ListRunner)->Center << "\t" << sqrt(size) << endl;
    }
  }
};

/** Returns the molecule with the given index \a index.
 * \param index index of the desired molecule
 * \return pointer to molecule structure, NULL if not found
 */
molecule * MoleculeListClass::ReturnIndex(int index)
{
  for(MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
    if ((*ListRunner)->IndexNr == index)
      return (*ListRunner);
  return NULL;
};

/** Simple merge of two molecules into one.
 * \param *mol destination molecule
 * \param *srcmol source molecule
 * \return true - merge successful, false - merge failed (probably due to non-existant indices
 */
bool MoleculeListClass::SimpleMerge(molecule *mol, molecule *srcmol)
{
  if (srcmol == NULL)
    return false;

  // put all molecules of src into mol
  atom *Walker = srcmol->start;
  atom *NextAtom = Walker->next;
  while (NextAtom != srcmol->end) {
    Walker = NextAtom;
    NextAtom = Walker->next;
    srcmol->UnlinkAtom(Walker);
    mol->AddAtom(Walker);
  }

  // remove src
  ListOfMolecules.remove(srcmol);
  delete(srcmol);
  return true;
};

/** Simple add of one molecules into another.
 * \param *mol destination molecule
 * \param *srcmol source molecule
 * \return true - merge successful, false - merge failed (probably due to non-existant indices
 */
bool MoleculeListClass::SimpleAdd(molecule *mol, molecule *srcmol)
{
  if (srcmol == NULL)
    return false;

  // put all molecules of src into mol
  atom *Walker = srcmol->start;
  atom *NextAtom = Walker->next;
  while (NextAtom != srcmol->end) {
    Walker = NextAtom;
    NextAtom = Walker->next;
    Walker = mol->AddCopyAtom(Walker);
    Walker->father = Walker;
  }

  return true;
};

/** Simple merge of a given set of molecules into one.
 * \param *mol destination molecule
 * \param *src index of set of source molecule
 * \param N number of source molecules
 * \return true - merge successful, false - some merges failed (probably due to non-existant indices)
 */
bool MoleculeListClass::SimpleMultiMerge(molecule *mol, int *src, int N)
{
  bool status = true;
  // check presence of all source molecules
  for (int i=0;i<N;i++) {
    molecule *srcmol = ReturnIndex(src[i]);
    status = status && SimpleMerge(mol, srcmol);
  }
  return status;
};

/** Simple add of a given set of molecules into one.
 * \param *mol destination molecule
 * \param *src index of set of source molecule
 * \param N number of source molecules
 * \return true - merge successful, false - some merges failed (probably due to non-existant indices)
 */
bool MoleculeListClass::SimpleMultiAdd(molecule *mol, int *src, int N)
{
  bool status = true;
  // check presence of all source molecules
  for (int i=0;i<N;i++) {
    molecule *srcmol = ReturnIndex(src[i]);
    status = status && SimpleAdd(mol, srcmol);
  }
  return status;
};

/** Scatter merge of a given set of molecules into one.
 * Scatter merge distributes the molecules in such a manner that they don't overlap.
 * \param *mol destination molecule
 * \param *src index of set of source molecule
 * \param N number of source molecules
 * \return true - merge successful, false - merge failed (probably due to non-existant indices
 * \TODO find scatter center for each src molecule
 */
bool MoleculeListClass::ScatterMerge(molecule *mol, int *src, int N)
{
  // check presence of all source molecules
  for (int i=0;i<N;i++) {
    // get pointer to src molecule
    molecule *srcmol = ReturnIndex(src[i]);
    if (srcmol == NULL)
      return false;
  }
  // adapt each Center
  for (int i=0;i<N;i++) {
    // get pointer to src molecule
    molecule *srcmol = ReturnIndex(src[i]);
    //srcmol->Center.Zero();
    srcmol->Translate(&srcmol->Center);
  }
  // perform a simple multi merge
  SimpleMultiMerge(mol, src, N);
  return true;
};

/** Embedding merge of a given set of molecules into one.
 * Embedding merge inserts one molecule into the other.
 * \param *mol destination molecule
 * \param *srcmol source molecule
 * \return true - merge successful, false - merge failed (probably due to non-existant indices
 * \TODO find embedding center
 */
bool MoleculeListClass::EmbedMerge(molecule *mol, molecule *srcmol)
{
  if (srcmol == NULL)
    return false;

  // calculate center for merge
  srcmol->Center.CopyVector(mol->FindEmbeddingHole((ofstream *)&cout, srcmol));
  srcmol->Center.Zero();

  // perform simple merge
  SimpleMerge(mol, srcmol);
  return true;
};

/** Simple output of the pointers in ListOfMolecules.
 * \param *out output stream
 */
void MoleculeListClass::Output(ofstream *out)
{
  *out << Verbose(1) << "MoleculeList: ";
  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
    *out << *ListRunner << "\t";
  *out << endl;
};

/** Calculates necessary hydrogen correction due to unwanted interaction between saturated ones.
 * If for a pair of two hydrogen atoms a and b, at least is a saturated one, and a and b are not
 * bonded to the same atom, then we add for this pair a correction term constructed from a Morse
 * potential function fit to QM calculations with respecting to the interatomic hydrogen distance.
 * \param *out output stream for debugging
 * \param *path path to file
 */
bool MoleculeListClass::AddHydrogenCorrection(ofstream *out, char *path)
{
  atom *Walker = NULL;
  atom *Runner = NULL;
  double ***FitConstant = NULL, **correction = NULL;
  int a, b;
  ofstream output;
  ifstream input;
  string line;
  stringstream zeile;
  double distance;
  char ParsedLine[1023];
  double tmp;
  char *FragmentNumber = NULL;

  cout << Verbose(1) << "Saving hydrogen saturation correction ... ";
  // 0. parse in fit constant files that should have the same dimension as the final energy files
  // 0a. find dimension of matrices with constants
  line = path;
  line.append("/");
  line += FRAGMENTPREFIX;
  line += "1";
  line += FITCONSTANTSUFFIX;
  input.open(line.c_str());
  if (input == NULL) {
    cerr << endl << "Unable to open " << line << ", is the directory correct?"
        << endl;
    return false;
  }
  a = 0;
  b = -1; // we overcount by one
  while (!input.eof()) {
    input.getline(ParsedLine, 1023);
    zeile.str(ParsedLine);
    int i = 0;
    while (!zeile.eof()) {
      zeile >> distance;
      i++;
    }
    if (i > a)
      a = i;
    b++;
  }
  cout << "I recognized " << a << " columns and " << b << " rows, ";
  input.close();

  // 0b. allocate memory for constants
  FitConstant = (double ***) Malloc(sizeof(double **) * 3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
  for (int k = 0; k < 3; k++) {
    FitConstant[k] = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
    for (int i = a; i--;) {
      FitConstant[k][i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
    }
  }
  // 0c. parse in constants
  for (int i = 0; i < 3; i++) {
    line = path;
    line.append("/");
    line += FRAGMENTPREFIX;
    sprintf(ParsedLine, "%d", i + 1);
    line += ParsedLine;
    line += FITCONSTANTSUFFIX;
    input.open(line.c_str());
    if (input == NULL) {
      cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
      return false;
    }
    int k = 0, l;
    while ((!input.eof()) && (k < b)) {
      input.getline(ParsedLine, 1023);
      //cout << "Current Line: " << ParsedLine << endl;
      zeile.str(ParsedLine);
      zeile.clear();
      l = 0;
      while ((!zeile.eof()) && (l < a)) {
        zeile >> FitConstant[i][l][k];
        //cout << FitConstant[i][l][k] << "\t";
        l++;
      }
      //cout << endl;
      k++;
    }
    input.close();
  }
  for (int k = 0; k < 3; k++) {
    cout << "Constants " << k << ":" << endl;
    for (int j = 0; j < b; j++) {
      for (int i = 0; i < a; i++) {
        cout << FitConstant[k][i][j] << "\t";
      }
      cout << endl;
    }
    cout << endl;
  }

  // 0d. allocate final correction matrix
  correction = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **correction");
  for (int i = a; i--;)
    correction[i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");

  // 1a. go through every molecule in the list
  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
    // 1b. zero final correction matrix
    for (int k = a; k--;)
      for (int j = b; j--;)
        correction[k][j] = 0.;
    // 2. take every hydrogen that is a saturated one
    Walker = (*ListRunner)->start;
    while (Walker->next != (*ListRunner)->end) {
      Walker = Walker->next;
      //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
      if ((Walker->type->Z == 1) && ((Walker->father == NULL)
          || (Walker->father->type->Z != 1))) { // if it's a hydrogen
        Runner = (*ListRunner)->start;
        while (Runner->next != (*ListRunner)->end) {
          Runner = Runner->next;
          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
          // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
          if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && ((*ListRunner)->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != (*ListRunner)->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) { // (hydrogens have only one bonding partner!)
            // 4. evaluate the morse potential for each matrix component and add up
            distance = Runner->x.Distance(&Walker->x);
            //cout << "Fragment " << (*ListRunner)->name << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
            for (int k = 0; k < a; k++) {
              for (int j = 0; j < b; j++) {
                switch (k) {
                  case 1:
                  case 7:
                  case 11:
                    tmp = pow(FitConstant[0][k][j] * (1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]))), 2);
                    break;
                  default:
                    tmp = FitConstant[0][k][j] * pow(distance, FitConstant[1][k][j]) + FitConstant[2][k][j];
                };
                correction[k][j] -= tmp; // ground state is actually lower (disturbed by additional interaction)
                //cout << tmp << "\t";
              }
              //cout << endl;
            }
            //cout << endl;
          }
        }
      }
    }
    // 5. write final matrix to file
    line = path;
    line.append("/");
    line += FRAGMENTPREFIX;
    FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), (*ListRunner)->IndexNr);
    line += FragmentNumber;
    delete (FragmentNumber);
    line += HCORRECTIONSUFFIX;
    output.open(line.c_str());
    output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
    for (int j = 0; j < b; j++) {
      for (int i = 0; i < a; i++)
        output << correction[i][j] << "\t";
      output << endl;
    }
    output.close();
  }
  line = path;
  line.append("/");
  line += HCORRECTIONSUFFIX;
  output.open(line.c_str());
  output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
  for (int j = 0; j < b; j++) {
    for (int i = 0; i < a; i++)
      output << 0 << "\t";
    output << endl;
  }
  output.close();
  // 6. free memory of parsed matrices
  FitConstant = (double ***) Malloc(sizeof(double **) * a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
  for (int k = 0; k < 3; k++) {
    FitConstant[k] = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
    for (int i = a; i--;) {
      FitConstant[k][i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
    }
  }
  cout << "done." << endl;
  return true;
};

/** 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 (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
      runner = (*ListRunner)->elemente->start;
      while (runner->next != (*ListRunner)->elemente->end) { // go through every element
        runner = runner->next;
        if ((*ListRunner)->ElementsInMolecule[runner->Z]) { // if this element got atoms
          Walker = (*ListRunner)->start;
          while (Walker->next != (*ListRunner)->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
 * \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(ofstream *out,
    config *configuration, int *SortIndex)
{
  ofstream outputFragment;
  char FragmentName[MAXSTRINGSIZE];
  char PathBackup[MAXSTRINGSIZE];
  bool result = true;
  bool intermediateResult = true;
  atom *Walker = NULL;
  Vector BoxDimension;
  char *FragmentNumber = NULL;
  char *path = NULL;
  int FragmentCounter = 0;
  ofstream output;

  // store the fragments as config and as xyz
  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
    // save default path as it is changed for each fragment
    path = configuration->GetDefaultPath();
    if (path != NULL)
      strcpy(PathBackup, path);
    else
      cerr << "OutputConfigForListOfFragments: NULL default path obtained from config!" << endl;

    // correct periodic
    (*ListRunner)->ScanForPeriodicCorrection(out);

    // output xyz file
    FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), 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 << "/" << FragmentCounter - 1 << " as XYZ ...";
    if ((intermediateResult = (*ListRunner)->OutputXYZ(&outputFragment)))
      *out << " done." << endl;
    else
      *out << " failed." << endl;
    result = result && intermediateResult;
    outputFragment.close();
    outputFragment.clear();

    // list atoms in fragment for debugging
    *out << Verbose(2) << "Contained atoms: ";
    Walker = (*ListRunner)->start;
    while (Walker->next != (*ListRunner)->end) {
      Walker = Walker->next;
      *out << Walker->Name << " ";
    }
    *out << endl;

    // center on edge
    (*ListRunner)->CenterEdge(out, &BoxDimension);
    (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary
    int j = -1;
    for (int k = 0; k < NDIM; k++) {
      j += k + 1;
      BoxDimension.x[k] = 2.5 * (configuration->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
      (*ListRunner)->cell_size[j] += BoxDimension.x[k] * 2.;
    }
    (*ListRunner)->Translate(&BoxDimension);

    // also calculate necessary orbitals
    (*ListRunner)->CountElements(); // this is a bugfix, atoms should shoulds actually be added correctly to this fragment
    (*ListRunner)->CalculateOrbitals(*configuration);

    // change path in config
    //strcpy(PathBackup, configuration->configpath);
    sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
    configuration->SetDefaultPath(FragmentName);

    // and save as config
    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as config ...";
    if ((intermediateResult = configuration->Save(FragmentName, (*ListRunner)->elemente, (*ListRunner))))
      *out << " done." << endl;
    else
      *out << " failed." << endl;
    result = result && intermediateResult;

    // restore old config
    configuration->SetDefaultPath(PathBackup);

    // and save as mpqc input file
    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as mpqc input ...";
    if ((intermediateResult = configuration->SaveMPQC(FragmentName, (*ListRunner))))
      *out << " done." << endl;
    else
      *out << " failed." << endl;

    result = result && intermediateResult;
    //outputFragment.close();
    //outputFragment.clear();
    delete (FragmentNumber);
    //Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
  }
  cout << " done." << endl;

  // printing final number
  *out << "Final number of fragments: " << FragmentCounter << "." << endl;

  return result;
};

/** Counts the number of molecules with the molecule::ActiveFlag set.
 * \return number of molecules with ActiveFlag set to true.
 */
int MoleculeListClass::NumberOfActiveMolecules()
{
  int count = 0;
  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
    count += ((*ListRunner)->ActiveFlag ? 1 : 0);
  return count;
};


/******************************************* 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;
};

/** Fills the bond structure of this chain list subgraphs that are derived from a complete \a *reference molecule.
 * Calls this routine in each MoleculeLeafClass::next subgraph if it's not NULL.
 * \param *out output stream for debugging
 * \param *reference reference molecule with the bond structure to be copied
 * \param &FragmentCounter Counter needed to address \a **ListOfLocalAtoms
 * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in \a *reference, may be NULL on start, then it is filled
 * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
 * \return true - success, false - faoilure
 */
bool MoleculeLeafClass::FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList)
{
  atom *Walker = NULL, *OtherWalker = NULL;
  bond *Binder = NULL;
  bool status = true;
  int AtomNo;

  *out << Verbose(1) << "Begin of FillBondStructureFromReference." << endl;
  // fill ListOfLocalAtoms if NULL was given
  if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
    *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
    return false;
  }

  if (status) {
    *out << Verbose(1) << "Creating adjacency list for subgraph " << this
        << "." << endl;
    Walker = Leaf->start;
    while (Walker->next != Leaf->end) {
      Walker = Walker->next;
      AtomNo = Walker->GetTrueFather()->nr; // global id of the current walker
      for (int i = 0; i < reference->NumberOfBondsPerAtom[AtomNo]; i++) { // go through father's bonds and copy them all
        Binder = reference->ListOfBondsPerAtom[AtomNo][i];
        OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->GetTrueFather())->nr]; // local copy of current bond partner of walker
        if (OtherWalker != NULL) {
          if (OtherWalker->nr > Walker->nr)
            Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);
        } else {
          *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
          status = false;
        }
      }
    }
    Leaf->CreateListOfBondsPerAtom(out);
    FragmentCounter++;
    if (next != NULL)
      status = next->FillBondStructureFromReference(out, reference, FragmentCounter, ListOfLocalAtoms);
    FragmentCounter--;
  }

  if ((FreeList) && (ListOfLocalAtoms != NULL)) {
    // free the index lookup list
    Free((void **) &ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::FillBondStructureFromReference - **ListOfLocalAtoms[]");
    if (FragmentCounter == 0) // first fragments frees the initial pointer to list
      Free((void **) &ListOfLocalAtoms, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
  }
  FragmentCounter--;
  *out << Verbose(1) << "End of FillBondStructureFromReference." << endl;
  return status;
};

/** Fills the root stack for sites to be used as root in fragmentation depending on order or adaptivity criteria
 * Again, as in \sa FillBondStructureFromReference steps recursively through each Leaf in this chain list of molecule's.
 * \param *out output stream for debugging
 * \param *&RootStack stack to be filled
 * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site
 * \param &FragmentCounter counts through the fragments in this MoleculeLeafClass
 * \return true - stack is non-empty, fragmentation necessary, false - stack is empty, no more sites to update
 */
bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out,
    KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter)
{
  atom *Walker = NULL, *Father = NULL;

  if (RootStack != NULL) {
    // find first root candidates
    if (&(RootStack[FragmentCounter]) != NULL) {
      RootStack[FragmentCounter].clear();
      Walker = Leaf->start;
      while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
        Walker = Walker->next;
        Father = Walker->GetTrueFather();
        if (AtomMask[Father->nr]) // apply mask
#ifdef ADDHYDROGEN
          if (Walker->type->Z != 1) // skip hydrogen
#endif
          RootStack[FragmentCounter].push_front(Walker->nr);
      }
      if (next != NULL)
        next->FillRootStackForSubgraphs(out, RootStack, AtomMask, ++FragmentCounter);
    } else {
      *out << Verbose(1) << "Rootstack[" << FragmentCounter << "] is NULL." << endl;
      return false;
    }
    FragmentCounter--;
    return true;
  } else {
    *out << Verbose(1) << "Rootstack is NULL." << endl;
    return false;
  }
};

/** Fills a lookup list of father's Atom::nr -> atom for each subgraph.
 * \param *out output stream fro debugging
 * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
 * \param FragmentCounter counts the fragments as we move along the list
 * \param GlobalAtomCount number of atoms in the complete molecule
 * \param &FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
 * \return true - succes, false - failure
 */
bool MoleculeLeafClass::FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, const int FragmentCounter, const int GlobalAtomCount, bool &FreeList)
{
  bool status = true;

  int Counter = Count();
  if (ListOfLocalAtoms == NULL) { // allocated initial pointer
    // allocate and set each field to NULL
    ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **) * Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
    if (ListOfLocalAtoms != NULL) {
      for (int i = Counter; i--;)
        ListOfLocalAtoms[i] = NULL;
      FreeList = FreeList && true;
    } else
      status = false;
  }

  if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph
    status = status && CreateFatherLookupTable(out, Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
    FreeList = FreeList && true;
  }

  return status;
};

/** The indices per keyset are compared to the respective father's Atom::nr in each subgraph and thus put into \a **&FragmentList.
 * \param *out output stream fro debugging
 * \param *reference reference molecule with the bond structure to be copied
 * \param *KeySetList list with all keysets
 * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
 * \param **&FragmentList list to be allocated and returned
 * \param &FragmentCounter counts the fragments as we move along the list
 * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
 * \retuen true - success, false - failure
 */
bool MoleculeLeafClass::AssignKeySetsToFragment(ofstream *out,
    molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms,
    Graph **&FragmentList, int &FragmentCounter, bool FreeList)
{
  bool status = true;
  int KeySetCounter = 0;

  *out << Verbose(1) << "Begin of AssignKeySetsToFragment." << endl;
  // fill ListOfLocalAtoms if NULL was given
  if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
    *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
    return false;
  }

  // allocate fragment list
  if (FragmentList == NULL) {
    KeySetCounter = Count();
    FragmentList = (Graph **) Malloc(sizeof(Graph *) * KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");
    for (int i = KeySetCounter; i--;)
      FragmentList[i] = NULL;
    KeySetCounter = 0;
  }

  if ((KeySetList != NULL) && (KeySetList->size() != 0)) { // if there are some scanned keysets at all
    // assign scanned keysets
    if (FragmentList[FragmentCounter] == NULL)
      FragmentList[FragmentCounter] = new Graph;
    KeySet *TempSet = new KeySet;
    for (Graph::iterator runner = KeySetList->begin(); runner != KeySetList->end(); runner++) { // key sets contain global numbers!
      if (ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*((*runner).first.begin()))->nr] != NULL) {// 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][reference->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) {// if there are no keysets, delete the list
      *out << Verbose(1) << "KeySetCounter is zero, deleting FragmentList." << endl;
      delete (FragmentList[FragmentCounter]);
    } else
      *out << Verbose(1) << KeySetCounter << " keysets were assigned to subgraph " << FragmentCounter << "." << endl;
    FragmentCounter++;
    if (next != NULL)
      next->AssignKeySetsToFragment(out, reference, KeySetList, ListOfLocalAtoms, FragmentList, FragmentCounter, FreeList);
    FragmentCounter--;
  } else
    *out << Verbose(1) << "KeySetList is NULL or empty." << endl;

  if ((FreeList) && (ListOfLocalAtoms != NULL)) {
    // free the index lookup list
    Free((void **) &ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::AssignKeySetsToFragment - **ListOfLocalAtoms[]");
    if (FragmentCounter == 0) // first fragments frees the initial pointer to list
      Free((void **) &ListOfLocalAtoms, "MoleculeLeafClass::AssignKeySetsToFragment - ***ListOfLocalAtoms");
  }
  *out << Verbose(1) << "End of AssignKeySetsToFragment." << endl;
  return status;
};

/** Translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
 * \param *out output stream for debugging
 * \param **FragmentList Graph with local numbers per fragment
 * \param &FragmentCounter counts the fragments as we move along the list
 * \param &TotalNumberOfKeySets global key set counter
 * \param &TotalGraph Graph to be filled with global numbers
 */
void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out,
    Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets,
    Graph &TotalGraph)
{
  *out << Verbose(1) << "Begin of TranslateIndicesToGlobalIDs." << endl;
  KeySet *TempSet = new KeySet;
  if (FragmentList[FragmentCounter] != NULL) {
    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((Leaf->FindAtom(*sprinter))->GetTrueFather()->nr);
      TotalGraph.insert(GraphPair(*TempSet, pair<int, double> (TotalNumberOfKeySets++, (*runner).second.second)));
      TempSet->clear();
    }
    delete (TempSet);
  } else {
    *out << Verbose(1) << "FragmentList is NULL." << endl;
  }
  if (next != NULL)
    next->TranslateIndicesToGlobalIDs(out, FragmentList, ++FragmentCounter, TotalNumberOfKeySets, TotalGraph);
  FragmentCounter--;
  *out << Verbose(1) << "End of TranslateIndicesToGlobalIDs." << endl;
};

/** Simply counts the number of items in the list, from given MoleculeLeafClass.
 * \return number of items
 */
int MoleculeLeafClass::Count() const
{
  if (next != NULL)
    return next->Count() + 1;
  else
    return 1;
};