/** \file parsing.cpp
 *
 * Declarations of various class functions for the parsing of value files.
 *
 */

// ======================================= INCLUDES ==========================================

#include "helpers.hpp"
#include "memoryallocator.hpp"
#include "parser.hpp"

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

// ======================================= FUNCTIONS ==========================================

/** Test if given filename can be opened.
 * \param filename name of file
 * \param test true - no error message, false - print error
 * \return given file exists
 */
bool FilePresent(const char *filename, bool test)
{
  ifstream input;

  input.open(filename, ios::in);
  if (input == NULL) {
    if (!test)
      Log() << Verbose(0) << endl << "Unable to open " << filename << ", is the directory correct?" << endl;
    return false;
  }
  input.close();
  return true;
};

/** Test the given parameters.
 * \param argc argument count
 * \param **argv arguments array
 * \return given inputdir is valid
 */
bool TestParams(int argc, char **argv)
{
  ifstream input;
  stringstream line;

  line << argv[1] << FRAGMENTPREFIX << KEYSETFILE;
  return FilePresent(line.str().c_str(), false);
};

// ======================================= CLASS MatrixContainer =============================

/** Constructor of MatrixContainer class.
 */
MatrixContainer::MatrixContainer() {
  Indices = NULL;
  Header = Malloc<char*>(1, "MatrixContainer::MatrixContainer: **Header");
  Matrix = Malloc<double**>(1, "MatrixContainer::MatrixContainer: ***Matrix"); // one more each for the total molecule
  RowCounter = Malloc<int>(1, "MatrixContainer::MatrixContainer: *RowCounter");
  ColumnCounter = Malloc<int>(1, "MatrixContainer::MatrixContainer: *ColumnCounter");
  Header[0] = NULL;
  Matrix[0] = NULL;
  RowCounter[0] = -1;
  MatrixCounter = 0;
  ColumnCounter[0] = -1;
};

/** Destructor of MatrixContainer class.
 */
MatrixContainer::~MatrixContainer() {
  if (Matrix != NULL) {
    for(int i=MatrixCounter;i--;) {
      if ((ColumnCounter != NULL) && (RowCounter != NULL)) {
          for(int j=RowCounter[i]+1;j--;)
            Free(&Matrix[i][j]);
        Free(&Matrix[i]);
      }
    }
    if ((ColumnCounter != NULL) && (RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
      for(int j=RowCounter[MatrixCounter]+1;j--;)
        Free(&Matrix[MatrixCounter][j]);
    if (MatrixCounter != 0)
      Free(&Matrix[MatrixCounter]);
    Free(&Matrix);
  }
  if (Indices != NULL)
    for(int i=MatrixCounter+1;i--;) {
      Free(&Indices[i]);
    }
  Free(&Indices);
  
  if (Header != NULL)
    for(int i=MatrixCounter+1;i--;)
      Free(&Header[i]);
  Free(&Header);
  Free(&RowCounter);
  Free(&ColumnCounter);
};

/** Either copies index matrix from another MatrixContainer or initialises with trivial mapping if NULL.
 * This either copies the index matrix or just maps 1 to 1, 2 to 2 and so on for all fragments.
 * \param *Matrix pointer to other MatrixContainer
 * \return true - copy/initialisation sucessful, false - dimension false for copying
 */
bool MatrixContainer::InitialiseIndices(class MatrixContainer *Matrix) 
{
  Log() << Verbose(0) << "Initialising indices";
  if (Matrix == NULL) {
    Log() << Verbose(0) << " with trivial mapping." << endl;
    Indices = Malloc<int*>(MatrixCounter + 1, "MatrixContainer::InitialiseIndices: **Indices");
    for(int i=MatrixCounter+1;i--;) {
      Indices[i] = Malloc<int>(RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
      for(int j=RowCounter[i];j--;)
        Indices[i][j] = j;
    }
  } else {
    Log() << Verbose(0) << " from other MatrixContainer." << endl;
    if (MatrixCounter != Matrix->MatrixCounter)
      return false;
    Indices = Malloc<int*>(MatrixCounter + 1, "MatrixContainer::InitialiseIndices: **Indices");
    for(int i=MatrixCounter+1;i--;) {
      if (RowCounter[i] != Matrix->RowCounter[i])
        return false;
      Indices[i] = Malloc<int>(Matrix->RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
      for(int j=Matrix->RowCounter[i];j--;) {
        Indices[i][j] = Matrix->Indices[i][j];
        //Log() << Verbose(0) << Indices[i][j] << "\t";
      }
      //Log() << Verbose(0) << endl;
    }
  }
  return true; 
};

/** Parsing a number of matrices.
 *    -# open the matrix file
 *    -# skip some lines (\a skiplines)
 *    -# scan header lines for number of columns
 *    -# scan lines for number of rows
 *    -# allocate matrix
 *    -# loop over found column and row counts and parse in each entry
 * \param *name directory with files
 * \param skiplines number of inital lines to skip
 * \param skiplines number of inital columns to skip
 * \param MatrixNr index number in Matrix array to parse into
 * \return parsing successful
 */
bool MatrixContainer::ParseMatrix(const char *name, int skiplines, int skipcolumns, int MatrixNr)
{
  ifstream input;
  stringstream line;
  string token;
  char filename[1023];

  input.open(name, ios::in);
  //Log() << Verbose(0) << "Opening " << name << " ... "  << input << endl;
  if (input == NULL) {
    eLog() << Verbose(0) << endl << "Unable to open " << name << ", is the directory correct?" << endl;
    performCriticalExit();
    return false;
  }

  // parse header
  Header[MatrixNr] = Malloc<char>(1024, "MatrixContainer::ParseMatrix: *Header[]");
  for (int m=skiplines+1;m--;)
    input.getline(Header[MatrixNr], 1023);
  
  // scan header for number of columns
  line.str(Header[MatrixNr]);
  for(int k=skipcolumns;k--;)
    line >> Header[MatrixNr];
  //Log() << Verbose(0) << line.str() << endl;
  ColumnCounter[MatrixNr]=0;
  while ( getline(line,token, '\t') ) {
    if (token.length() > 0)
      ColumnCounter[MatrixNr]++;
  }
  //Log() << Verbose(0) << line.str() << endl;
  //Log() << Verbose(0) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
  if (ColumnCounter[MatrixNr] == 0) {
    eLog() << Verbose(0) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
    performCriticalExit();
  }
  
  // scan rest for number of rows/lines
  RowCounter[MatrixNr]=-1;    // counts one line too much
  while (!input.eof()) {
    input.getline(filename, 1023);
    //Log() << Verbose(0) << "Comparing: " << strncmp(filename,"MeanForce",9) << endl;
    RowCounter[MatrixNr]++; // then line was not last MeanForce
    if (strncmp(filename,"MeanForce", 9) == 0) {
      break;
    }
  }
  //Log() << Verbose(0) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
  if (RowCounter[MatrixNr] == 0) {
    eLog() << Verbose(0) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
    performCriticalExit();
  }

  // allocate matrix if it's not zero dimension in one direction
  if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
    Matrix[MatrixNr] = Malloc<double*>(RowCounter[MatrixNr] + 1, "MatrixContainer::ParseMatrix: **Matrix[]");
  
    // parse in each entry for this matrix
    input.clear();
    input.seekg(ios::beg);
    for (int m=skiplines+1;m--;)
      input.getline(Header[MatrixNr], 1023);    // skip header
    line.str(Header[MatrixNr]);
    for(int k=skipcolumns;k--;)  // skip columns in header too
      line >> filename;
    strncpy(Header[MatrixNr], line.str().c_str(), 1023);  
    for(int j=0;j<RowCounter[MatrixNr];j++) {
      Matrix[MatrixNr][j] = Malloc<double>(ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[][]");
      input.getline(filename, 1023);
      stringstream lines(filename);
      //Log() << Verbose(0) << "Matrix at level " << j << ":";// << filename << endl;
      for(int k=skipcolumns;k--;)
        lines >> filename;
      for(int k=0;(k<ColumnCounter[MatrixNr]) && (!lines.eof());k++) {
        lines >> Matrix[MatrixNr][j][k];
        //Log() << Verbose(0) << " " << setprecision(2) << Matrix[MatrixNr][j][k];;
      }
      //Log() << Verbose(0) << endl;
      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = Malloc<double>(ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[RowCounter[MatrixNr]][]");
      for(int j=ColumnCounter[MatrixNr];j--;)
        Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
    }
  } else {
    eLog() << Verbose(1) << "Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl; 
  }
  input.close();
  return true;
};

/** Parsing a number of matrices.
 * -# First, count the number of matrices by counting lines in KEYSETFILE
 * -# Then,
 *    -# construct the fragment number
 *    -# open the matrix file
 *    -# skip some lines (\a skiplines)
 *    -# scan header lines for number of columns
 *    -# scan lines for number of rows
 *    -# allocate matrix
 *    -# loop over found column and row counts and parse in each entry
 * -# Finally, allocate one additional matrix (\a MatrixCounter) containing combined or temporary values
 * \param *name directory with files
 * \param *prefix prefix of each matrix file
 * \param *suffix suffix of each matrix file
 * \param skiplines number of inital lines to skip
 * \param skiplines number of inital columns to skip
 * \return parsing successful
 */
bool MatrixContainer::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
{
  char filename[1023];
  ifstream input;
  char *FragmentNumber = NULL;
  stringstream file;
  string token;

  // count the number of matrices
  MatrixCounter = -1; // we count one too much
  file << name << FRAGMENTPREFIX << KEYSETFILE;
  input.open(file.str().c_str(), ios::in);
  if (input == NULL) {
    Log() << Verbose(0) << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
    return false;
  }
  while (!input.eof()) {
    input.getline(filename, 1023);
    stringstream zeile(filename);
    MatrixCounter++;
  }
  input.close();
  Log() << Verbose(0) << "Determined " << MatrixCounter << " fragments." << endl;

  Log() << Verbose(0) << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl;
  Header = ReAlloc<char*>(Header, MatrixCounter + 1, "MatrixContainer::ParseFragmentMatrix: **Header"); // one more each for the total molecule
  Matrix = ReAlloc<double**>(Matrix, MatrixCounter + 1, "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
  RowCounter = ReAlloc<int>(RowCounter, MatrixCounter + 1, "MatrixContainer::ParseFragmentMatrix: *RowCounter");
  ColumnCounter = ReAlloc<int>(ColumnCounter, MatrixCounter + 1, "MatrixContainer::ParseFragmentMatrix: *ColumnCounter");
  for(int i=MatrixCounter+1;i--;) {
    Matrix[i] = NULL;
    Header[i] = NULL;
    RowCounter[i] = -1;
    ColumnCounter[i] = -1;
  }
  for(int i=0; i < MatrixCounter;i++) {
    // open matrix file
    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    file.str(" ");
    file << name << FRAGMENTPREFIX << FragmentNumber << prefix << suffix;
    if (!ParseMatrix(file.str().c_str(), skiplines, skipcolumns, i))
      return false;
    Free(&FragmentNumber);
  }
  return true;
};

/** Allocates and resets the memory for a number \a MCounter of matrices.
 * \param **GivenHeader Header line for each matrix
 * \param MCounter number of matrices
 * \param *RCounter number of rows for each matrix
 * \param *CCounter number of columns for each matrix
 * \return Allocation successful 
 */
bool MatrixContainer::AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter)
{
  MatrixCounter = MCounter;
  Header = Malloc<char*>(MatrixCounter + 1, "MatrixContainer::AllocateMatrix: *Header");
  Matrix = Malloc<double**>(MatrixCounter + 1, "MatrixContainer::AllocateMatrix: ***Matrix"); // one more each for the total molecule
  RowCounter = Malloc<int>(MatrixCounter + 1, "MatrixContainer::AllocateMatrix: *RowCounter");
  ColumnCounter = Malloc<int>(MatrixCounter + 1, "MatrixContainer::AllocateMatrix: *ColumnCounter");
  for(int i=MatrixCounter+1;i--;) {
    Header[i] = Malloc<char>(1024, "MatrixContainer::AllocateMatrix: *Header[i]");
    strncpy(Header[i], GivenHeader[i], 1023);
    RowCounter[i] = RCounter[i];
    ColumnCounter[i] = CCounter[i];
    Matrix[i] = Malloc<double*>(RowCounter[i] + 1, "MatrixContainer::AllocateMatrix: **Matrix[]");
    for(int j=RowCounter[i]+1;j--;) {
      Matrix[i][j] = Malloc<double>(ColumnCounter[i], "MatrixContainer::AllocateMatrix: *Matrix[][]");
      for(int k=ColumnCounter[i];k--;)
        Matrix[i][j][k] = 0.;
    }
  }
  return true;
};

/** Resets all values in MatrixContainer::Matrix.
 * \return true if successful
 */
bool MatrixContainer::ResetMatrix()
{
  for(int i=MatrixCounter+1;i--;)
    for(int j=RowCounter[i]+1;j--;)
      for(int k=ColumnCounter[i];k--;)
        Matrix[i][j][k] = 0.;
   return true;
};

/** Scans all elements of MatrixContainer::Matrix for greatest absolute value.
 * \return greatest value of MatrixContainer::Matrix
 */
double MatrixContainer::FindMaxValue()
{
  double max = Matrix[0][0][0];
  for(int i=MatrixCounter+1;i--;)
    for(int j=RowCounter[i]+1;j--;)
      for(int k=ColumnCounter[i];k--;)
        if (fabs(Matrix[i][j][k]) > max)
          max = fabs(Matrix[i][j][k]);
  if (fabs(max) < MYEPSILON)
    max += MYEPSILON;
 return max;
};

/** Scans all elements of MatrixContainer::Matrix for smallest absolute value.
 * \return smallest value of MatrixContainer::Matrix
 */
double MatrixContainer::FindMinValue()
{
  double min = Matrix[0][0][0];
  for(int i=MatrixCounter+1;i--;)
    for(int j=RowCounter[i]+1;j--;)
      for(int k=ColumnCounter[i];k--;)
        if (fabs(Matrix[i][j][k]) < min)
          min = fabs(Matrix[i][j][k]);
  if (fabs(min) < MYEPSILON)
    min += MYEPSILON;
  return min;
};

/** Sets all values in the last of MatrixContainer::Matrix to \a value.
 * \param value reset value
 * \param skipcolumns skip initial columns
 * \return true if successful
 */
bool MatrixContainer::SetLastMatrix(double value, int skipcolumns)
{
  for(int j=RowCounter[MatrixCounter]+1;j--;)
    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
      Matrix[MatrixCounter][j][k] = value;
   return true;
};

/** Sets all values in the last of MatrixContainer::Matrix to \a value.
 * \param **values matrix with each value (must have at least same dimensions!)
 * \param skipcolumns skip initial columns
 * \return true if successful
 */
bool MatrixContainer::SetLastMatrix(double **values, int skipcolumns)
{
  for(int j=RowCounter[MatrixCounter]+1;j--;)
    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
      Matrix[MatrixCounter][j][k] = values[j][k];
   return true;
};

/** Sums the entries with each factor and put into last element of \a ***Matrix.
 * Sums over "E"-terms to create the "F"-terms
 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
 * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
 * \param Order bond order
 * \return true if summing was successful
 */
bool MatrixContainer::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySets, int Order)
{
  // go through each order
  for (int CurrentFragment=0;CurrentFragment<KeySets.FragmentsPerOrder[Order];CurrentFragment++) {
    //Log() << Verbose(0) << "Current Fragment is " << CurrentFragment << "/" << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
    // then go per order through each suborder and pick together all the terms that contain this fragment
    for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
      for (int j=0;j<KeySets.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
        if (KeySets.Contains(KeySets.OrderSet[Order][CurrentFragment], KeySets.OrderSet[SubOrder][j])) {
          //Log() << Verbose(0) << "Current other fragment is " << j << "/" << KeySets.OrderSet[SubOrder][j] << "." << endl;
          // if the fragment's indices are all in the current fragment
          for(int k=0;k<RowCounter[ KeySets.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
            int m = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][k];
            //Log() << Verbose(0) << "Current index is " << k << "/" << m << "." << endl;
            if (m != -1) { // if it's not an added hydrogen
              for (int l=0;l<RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
                //Log() << Verbose(0) << "Comparing " << m << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
                if (m == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l]) {
                  m = l;
                  break;
                }
              }
              //Log() << Verbose(0) << "Corresponding index in CurrentFragment is " << m << "." << endl;
              if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
                eLog() << Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current force index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl;
                performCriticalExit();
                return false;
              }
              if (Order == SubOrder) { // equal order is always copy from Energies
                for(int l=ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l--;) // then adds/subtract each column
                  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
              } else {
                for(int l=ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l--;)
                  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
              }
            }
            //if ((ColumnCounter[ KeySets.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
             //Log() << Verbose(0) << "Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
          }
        } else {
          //Log() << Verbose(0) << "Fragment " << KeySets.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
        }
      }
    }
   //Log() << Verbose(0) << "Final Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << KeySets.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][KeySets.AtomCounter[0]-1][1] << endl;
  }

  return true;
};

/** Writes the summed total fragment terms \f$F_{ij}\f$ to file.
 * \param *name inputdir
 * \param *prefix prefix before \a EnergySuffix
 * \return file was written
 */
bool MatrixContainer::WriteTotalFragments(const char *name, const char *prefix)
{
  ofstream output;
  char *FragmentNumber = NULL;

  Log() << Verbose(0) << "Writing fragment files." << endl;
  for(int i=0;i<MatrixCounter;i++) {
    stringstream line;
    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    line << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix;
    Free(&FragmentNumber);
    output.open(line.str().c_str(), ios::out);
    if (output == NULL) {
      eLog() << Verbose(0) << "Unable to open output energy file " << line.str() << "!" << endl;
      performCriticalExit();
      return false;
    }
    output << Header[i] << endl;
    for(int j=0;j<RowCounter[i];j++) {
      for(int k=0;k<ColumnCounter[i];k++)
        output << scientific << Matrix[i][j][k] << "\t";
      output << endl;
    }
    output.close();
  }
  return true;
};

/** Writes the summed total values in the last matrix to file.
 * \param *name inputdir
 * \param *prefix prefix
 * \param *suffix suffix
 * \return file was written
 */
bool MatrixContainer::WriteLastMatrix(const char *name, const char *prefix, const char *suffix)
{
  ofstream output;
  stringstream line;

  Log() << Verbose(0) << "Writing matrix values of " << suffix << "." << endl;
  line << name << prefix << suffix;
  output.open(line.str().c_str(), ios::out);
  if (output == NULL) {
    eLog() << Verbose(0) << "Unable to open output matrix file " << line.str() << "!" << endl;
    performCriticalExit();
    return false;
  }
  output << Header[MatrixCounter] << endl;
  for(int j=0;j<RowCounter[MatrixCounter];j++) {
    for(int k=0;k<ColumnCounter[MatrixCounter];k++)
      output << scientific << Matrix[MatrixCounter][j][k] << "\t";
    output << endl;
  }
  output.close();
  return true;
};

// ======================================= CLASS EnergyMatrix =============================

/** Create a trivial energy index mapping.
 * This just maps 1 to 1, 2 to 2 and so on for all fragments.
 * \return creation sucessful
 */
bool EnergyMatrix::ParseIndices()
{
  Log() << Verbose(0) << "Parsing energy indices." << endl;
  Indices = Malloc<int*>(MatrixCounter + 1, "EnergyMatrix::ParseIndices: **Indices");
  for(int i=MatrixCounter+1;i--;) {
    Indices[i] = Malloc<int>(RowCounter[i], "EnergyMatrix::ParseIndices: *Indices[]");
    for(int j=RowCounter[i];j--;)
      Indices[i][j] = j;
  }
  return true;
};

/** Sums the energy with each factor and put into last element of \a EnergyMatrix::Matrix.
 * Sums over the "F"-terms in ANOVA decomposition.
 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
 * \param CorrectionFragments MatrixContainer with hydrogen saturation correction per fragments
 * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
 * \param Order bond order
 * \parsm sign +1 or -1
 * \return true if summing was successful
 */
bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySets, int Order, double sign)
{
  // sum energy
  if (CorrectionFragments == NULL)
    for(int i=KeySets.FragmentsPerOrder[Order];i--;)
      for(int j=RowCounter[ KeySets.OrderSet[Order][i] ];j--;)
        for(int k=ColumnCounter[ KeySets.OrderSet[Order][i] ];k--;)
          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySets.OrderSet[Order][i] ][j][k];
  else
    for(int i=KeySets.FragmentsPerOrder[Order];i--;)
      for(int j=RowCounter[ KeySets.OrderSet[Order][i] ];j--;)
        for(int k=ColumnCounter[ KeySets.OrderSet[Order][i] ];k--;)
          Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySets.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySets.OrderSet[Order][i] ][j][k]);
  return true;
};

/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
 * \param *name directory with files
 * \param *prefix prefix of each matrix file
 * \param *suffix suffix of each matrix file
 * \param skiplines number of inital lines to skip
 * \param skiplines number of inital columns to skip
 * \return parsing successful
 */
bool EnergyMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
{
  char filename[1024];
  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);

  if (status) {
    // count maximum of columns
    RowCounter[MatrixCounter] = 0;
    ColumnCounter[MatrixCounter] = 0;
    for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
      if (RowCounter[j] > RowCounter[MatrixCounter])
        RowCounter[MatrixCounter] = RowCounter[j];
      if (ColumnCounter[j] > ColumnCounter[MatrixCounter])  // take maximum of all for last matrix
        ColumnCounter[MatrixCounter] = ColumnCounter[j];
    }
    // allocate last plus one matrix
    Log() << Verbose(0) << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
    Matrix[MatrixCounter] = Malloc<double*>(RowCounter[MatrixCounter] + 1, "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    for(int j=0;j<=RowCounter[MatrixCounter];j++)
      Matrix[MatrixCounter][j] = Malloc<double>(ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    
    // try independently to parse global energysuffix file if present
    strncpy(filename, name, 1023);
    strncat(filename, prefix, 1023-strlen(filename));
    strncat(filename, suffix.c_str(), 1023-strlen(filename));
    ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
  }
  return status;
};

// ======================================= CLASS ForceMatrix =============================

/** Parsing force Indices of each fragment
 * \param *name directory with \a ForcesFile
 * \return parsing successful
 */
bool ForceMatrix::ParseIndices(const char *name)
{
  ifstream input;
  char *FragmentNumber = NULL;
  char filename[1023];
  stringstream line;

  Log() << Verbose(0) << "Parsing force indices for " << MatrixCounter << " matrices." << endl;
  Indices = Malloc<int*>(MatrixCounter + 1, "ForceMatrix::ParseIndices: **Indices");
  line << name << FRAGMENTPREFIX << FORCESFILE;
  input.open(line.str().c_str(), ios::in);
  //Log() << Verbose(0) << "Opening " << line.str() << " ... "  << input << endl;
  if (input == NULL) {
    Log() << Verbose(0) << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
    return false;
  }
  for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
    // get the number of atoms for this fragment
    input.getline(filename, 1023);
    line.str(filename);
    // parse the values
    Indices[i] = Malloc<int>(RowCounter[i], "ForceMatrix::ParseIndices: *Indices[]");
    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    //Log() << Verbose(0) << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
    Free(&FragmentNumber);
    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
      line >> Indices[i][j];
      //Log() << Verbose(0) << " " << Indices[i][j];
    }
    //Log() << Verbose(0) << endl;
  }
  Indices[MatrixCounter] = Malloc<int>(RowCounter[MatrixCounter], "ForceMatrix::ParseIndices: *Indices[]");
  for(int j=RowCounter[MatrixCounter];j--;) {
    Indices[MatrixCounter][j] = j;
  }
  input.close();
  return true;
};


/** Sums the forces and puts into last element of \a ForceMatrix::Matrix.
 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
 * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
 * \param Order bond order
 *  \param sign +1 or -1
 * \return true if summing was successful
 */
bool ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySets, int Order, double sign)
{
  int FragmentNr;
  // sum forces
  for(int i=0;i<KeySets.FragmentsPerOrder[Order];i++) {
    FragmentNr = KeySets.OrderSet[Order][i];
    for(int l=0;l<RowCounter[ FragmentNr ];l++) {
      int j = Indices[ FragmentNr ][l];
      if (j > RowCounter[MatrixCounter]) {
        eLog() << Verbose(0) << "Current force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
        performCriticalExit();
        return false;
      }
      if (j != -1) {
        //if (j == 0) Log() << Verbose(0) << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
        for(int k=2;k<ColumnCounter[MatrixCounter];k++)
          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
      }
    }
  }
  return true;
};


/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
 * \param *name directory with files
 * \param *prefix prefix of each matrix file
 * \param *suffix suffix of each matrix file
 * \param skiplines number of inital lines to skip
 * \param skiplines number of inital columns to skip
 * \return parsing successful
 */
bool ForceMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
{
  char filename[1023];
  ifstream input;
  stringstream file;
  int nr;
  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);

  if (status) {
    // count number of atoms for last plus one matrix
    file << name << FRAGMENTPREFIX << KEYSETFILE;
    input.open(file.str().c_str(), ios::in);
    if (input == NULL) {
      Log() << Verbose(0) << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
      return false;
    }
    RowCounter[MatrixCounter] = 0;
    while (!input.eof()) {
      input.getline(filename, 1023);
      stringstream zeile(filename);
      while (!zeile.eof()) {
        zeile >> nr;
        //Log() << Verbose(0) << "Current index: " << nr << "." << endl;
        if (nr > RowCounter[MatrixCounter])
          RowCounter[MatrixCounter] = nr;
      }
    }
    RowCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
    input.close();

    ColumnCounter[MatrixCounter] = 0;
    for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
      if (ColumnCounter[j] > ColumnCounter[MatrixCounter])  // take maximum of all for last matrix
        ColumnCounter[MatrixCounter] = ColumnCounter[j];
    }
  
    // allocate last plus one matrix
    Log() << Verbose(0) << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
    Matrix[MatrixCounter] = Malloc<double*>(RowCounter[MatrixCounter] + 1, "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    for(int j=0;j<=RowCounter[MatrixCounter];j++)
      Matrix[MatrixCounter][j] = Malloc<double>(ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");

    // try independently to parse global forcesuffix file if present
    strncpy(filename, name, 1023);
    strncat(filename, prefix, 1023-strlen(filename));
    strncat(filename, suffix.c_str(), 1023-strlen(filename));
    ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
  }
 

  return status;
};

// ======================================= CLASS HessianMatrix =============================

/** Parsing force Indices of each fragment
 * \param *name directory with \a ForcesFile
 * \return parsing successful 
 */
bool HessianMatrix::ParseIndices(char *name) 
{
  ifstream input;
  char *FragmentNumber = NULL;
  char filename[1023];
  stringstream line;
  
  Log() << Verbose(0) << "Parsing hessian indices for " << MatrixCounter << " matrices." << endl;
  Indices = Malloc<int*>(MatrixCounter + 1, "HessianMatrix::ParseIndices: **Indices");
  line << name << FRAGMENTPREFIX << FORCESFILE;
  input.open(line.str().c_str(), ios::in);
  //Log() << Verbose(0) << "Opening " << line.str() << " ... "  << input << endl;
  if (input == NULL) {
    Log() << Verbose(0) << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
    return false;
  }
  for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
    // get the number of atoms for this fragment
    input.getline(filename, 1023);
    line.str(filename);
    // parse the values
    Indices[i] = Malloc<int>(RowCounter[i], "HessianMatrix::ParseIndices: *Indices[]");
    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    //Log() << Verbose(0) << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
    Free(&FragmentNumber);
    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
      line >> Indices[i][j];
      //Log() << Verbose(0) << " " << Indices[i][j];
    }
    //Log() << Verbose(0) << endl;
  }
  Indices[MatrixCounter] = Malloc<int>(RowCounter[MatrixCounter], "HessianMatrix::ParseIndices: *Indices[]");
  for(int j=RowCounter[MatrixCounter];j--;) {
    Indices[MatrixCounter][j] = j;
  }
  input.close();
  return true;
};


/** Sums the hessian entries and puts into last element of \a HessianMatrix::Matrix.
 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
 * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
 * \param Order bond order
 *  \param sign +1 or -1
 * \return true if summing was successful
 */
bool HessianMatrix::SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySets, int Order, double sign)
{
  int FragmentNr;
  // sum forces
  for(int i=0;i<KeySets.FragmentsPerOrder[Order];i++) {
    FragmentNr = KeySets.OrderSet[Order][i];
    for(int l=0;l<RowCounter[ FragmentNr ];l++) {
      int j = Indices[ FragmentNr ][l];
      if (j > RowCounter[MatrixCounter]) {
        eLog() << Verbose(0) << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
        performCriticalExit();
        return false;
      }
      if (j != -1) {
        for(int m=0;m<ColumnCounter[ FragmentNr ];m++) {
          int k = Indices[ FragmentNr ][m];
          if (k > ColumnCounter[MatrixCounter]) {
            eLog() << Verbose(0) << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
            performCriticalExit();
            return false;
          }
          if (k != -1) {
            //Log() << Verbose(0) << "Adding " << sign*Fragments.Matrix[ FragmentNr ][l][m] << " from [" << l << "][" << m << "] onto [" << j << "][" << k << "]." << endl;
            Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][m];
          }
        }
      }
    }
  }
  return true;
};

/** Constructor for class HessianMatrix.
 */
HessianMatrix::HessianMatrix() : MatrixContainer() 
{
   IsSymmetric = true;
}

/** Sums the hessian entries with each factor and put into last element of \a ***Matrix.
 * Sums over "E"-terms to create the "F"-terms
 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
 * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
 * \param Order bond order
 * \return true if summing was successful
 */
bool HessianMatrix::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySets, int Order)
{
  // go through each order
  for (int CurrentFragment=0;CurrentFragment<KeySets.FragmentsPerOrder[Order];CurrentFragment++) {
    //Log() << Verbose(0) << "Current Fragment is " << CurrentFragment << "/" << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
    // then go per order through each suborder and pick together all the terms that contain this fragment
    for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
      for (int j=0;j<KeySets.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
        if (KeySets.Contains(KeySets.OrderSet[Order][CurrentFragment], KeySets.OrderSet[SubOrder][j])) {
          //Log() << Verbose(0) << "Current other fragment is " << j << "/" << KeySets.OrderSet[SubOrder][j] << "." << endl;
          // if the fragment's indices are all in the current fragment
          for(int k=0;k<RowCounter[ KeySets.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
            int m = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][k];
            //Log() << Verbose(0) << "Current row index is " << k << "/" << m << "." << endl;
            if (m != -1) { // if it's not an added hydrogen
              for (int l=0;l<RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
                //Log() << Verbose(0) << "Comparing " << m << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
                if (m == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l]) {
                  m = l;
                  break;  
                }
              }
              //Log() << Verbose(0) << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl;
              if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
                eLog() << Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl;
                performCriticalExit();
                return false;
              }
              
              for(int l=0;l<ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l++) {
                int n = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][l];
                //Log() << Verbose(0) << "Current column index is " << l << "/" << n << "." << endl;
                if (n != -1) { // if it's not an added hydrogen
                  for (int p=0;p<ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment
                    //Log() << Verbose(0) << "Comparing " << n << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][p] << "." << endl;
                    if (n == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][p]) {
                      n = p;
                      break;  
                    }
                  }
                  //Log() << Verbose(0) << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl;
                  if (n > ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
                    eLog() << Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl;
                    performCriticalExit();
                    return false;
                  }
                  if (Order == SubOrder) { // equal order is always copy from Energies
                    //Log() << Verbose(0) << "Adding " << MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
                    Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
                  } else {
                    //Log() << Verbose(0) << "Subtracting " << Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
                    Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
                  }
                }
              }
            }
            //if ((ColumnCounter[ KeySets.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
             //Log() << Verbose(0) << "Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
          }
        } else {
          //Log() << Verbose(0) << "Fragment " << KeySets.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
        }
      }
    }
   //Log() << Verbose(0) << "Final Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << KeySets.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][KeySets.AtomCounter[0]-1][1] << endl;
  }
  
  return true;
};

/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
 * \param *name directory with files
 * \param *prefix prefix of each matrix file
 * \param *suffix suffix of each matrix file
 * \param skiplines number of inital lines to skip
 * \param skiplines number of inital columns to skip
 * \return parsing successful
 */ 
bool HessianMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
{
  char filename[1023];
  ifstream input;
  stringstream file;
  int nr;
  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);

  if (status) {
    // count number of atoms for last plus one matrix
    file << name << FRAGMENTPREFIX << KEYSETFILE;
    input.open(file.str().c_str(), ios::in);
    if (input == NULL) {
      Log() << Verbose(0) << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
      return false;
    }
    RowCounter[MatrixCounter] = 0;
    ColumnCounter[MatrixCounter] = 0;
    while (!input.eof()) {
      input.getline(filename, 1023);
      stringstream zeile(filename);
      while (!zeile.eof()) {
        zeile >> nr;
        //Log() << Verbose(0) << "Current index: " << nr << "." << endl;
        if (nr > RowCounter[MatrixCounter]) {
          RowCounter[MatrixCounter] = nr;
          ColumnCounter[MatrixCounter] = nr;
        }
      }
    }
    RowCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
    ColumnCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
    input.close();
  
    // allocate last plus one matrix
    Log() << Verbose(0) << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
    Matrix[MatrixCounter] = Malloc<double*>(RowCounter[MatrixCounter] + 1, "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    for(int j=0;j<=RowCounter[MatrixCounter];j++)
      Matrix[MatrixCounter][j] = Malloc<double>(ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");

    // try independently to parse global forcesuffix file if present
    strncpy(filename, name, 1023);
    strncat(filename, prefix, 1023-strlen(filename));
    strncat(filename, suffix.c_str(), 1023-strlen(filename));
    ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
  }


  return status;
};

// ======================================= CLASS KeySetsContainer =============================

/** Constructor of KeySetsContainer class.
 */
KeySetsContainer::KeySetsContainer() {
  KeySets = NULL;
  AtomCounter = NULL;
  FragmentCounter = 0;
  Order = 0;
  FragmentsPerOrder = 0;
  OrderSet = NULL;
};

/** Destructor of KeySetsContainer class.
 */
KeySetsContainer::~KeySetsContainer() {
  for(int i=FragmentCounter;i--;)
    Free(&KeySets[i]);
  for(int i=Order;i--;)
    Free(&OrderSet[i]);
  Free(&KeySets);
  Free(&OrderSet);
  Free(&AtomCounter);
  Free(&FragmentsPerOrder);
};

/** Parsing KeySets into array.
 * \param *name directory with keyset file
 * \param *ACounter number of atoms per fragment
 * \param FCounter number of fragments
 * \return parsing succesful
 */
bool KeySetsContainer::ParseKeySets(const char *name, const int *ACounter, const int FCounter) {
  ifstream input;
  char *FragmentNumber = NULL;
  stringstream file;
  char filename[1023];

  FragmentCounter = FCounter;
  Log() << Verbose(0) << "Parsing key sets." << endl;
  KeySets = Malloc<int*>(FragmentCounter, "KeySetsContainer::ParseKeySets: **KeySets");
  for(int i=FragmentCounter;i--;)
    KeySets[i] = NULL;
  file << name << FRAGMENTPREFIX << KEYSETFILE;
  input.open(file.str().c_str(), ios::in);
  if (input == NULL) {
    Log() << Verbose(0) << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
    return false;
  }

  AtomCounter = Malloc<int>(FragmentCounter, "KeySetsContainer::ParseKeySets: *RowCounter");
  for(int i=0;(i<FragmentCounter) && (!input.eof());i++) {
    stringstream line;
    AtomCounter[i] = ACounter[i];
    // parse the values
    KeySets[i] = Malloc<int>(AtomCounter[i], "KeySetsContainer::ParseKeySets: *KeySets[]");
    for(int j=AtomCounter[i];j--;)
      KeySets[i][j] = -1;
    FragmentNumber = FixedDigitNumber(FragmentCounter, i);
    //Log() << Verbose(0) << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
    Free(&FragmentNumber);
    input.getline(filename, 1023);
    line.str(filename);
    for(int j=0;(j<AtomCounter[i]) && (!line.eof());j++) {
      line >> KeySets[i][j];
      //Log() << Verbose(0) << " " << KeySets[i][j];
    }
    //Log() << Verbose(0) << endl;
  }
  input.close();
  return true;
};

/** Parse many body terms, associating each fragment to a certain bond order.
 * \return parsing succesful
 */
bool KeySetsContainer::ParseManyBodyTerms()
{
  int Counter;

  Log() << Verbose(0) << "Creating Fragment terms." << endl;
  // scan through all to determine maximum order
  Order=0;
  for(int i=FragmentCounter;i--;) {
    Counter=0;
    for(int j=AtomCounter[i];j--;)
      if (KeySets[i][j] != -1)
        Counter++;
    if (Counter > Order)
      Order = Counter;
  }
  Log() << Verbose(0) << "Found Order is " << Order << "." << endl;

  // scan through all to determine fragments per order
  FragmentsPerOrder = Malloc<int>(Order, "KeySetsContainer::ParseManyBodyTerms: *FragmentsPerOrder");
  for(int i=Order;i--;)
    FragmentsPerOrder[i] = 0;
  for(int i=FragmentCounter;i--;) {
    Counter=0;
    for(int j=AtomCounter[i];j--;)
      if (KeySets[i][j] != -1)
        Counter++;
    FragmentsPerOrder[Counter-1]++;
  }
  for(int i=0;i<Order;i++)
    Log() << Verbose(0) << "Found No. of Fragments of Order " << i+1 << " is " << FragmentsPerOrder[i] << "." << endl;

  // scan through all to gather indices to each order set
  OrderSet = Malloc<int*>(Order, "KeySetsContainer::ParseManyBodyTerms: **OrderSet");
  for(int i=Order;i--;)
    OrderSet[i] = Malloc<int>(FragmentsPerOrder[i], "KeySetsContainer::ParseManyBodyTermsKeySetsContainer::ParseManyBodyTerms: *OrderSet[]");
  for(int i=Order;i--;)
    FragmentsPerOrder[i] = 0;
  for(int i=FragmentCounter;i--;) {
    Counter=0;
    for(int j=AtomCounter[i];j--;)
      if (KeySets[i][j] != -1)
        Counter++;
    OrderSet[Counter-1][FragmentsPerOrder[Counter-1]] = i;
    FragmentsPerOrder[Counter-1]++;
  }
  Log() << Verbose(0) << "Printing OrderSet." << endl;
  for(int i=0;i<Order;i++) {
    for (int j=0;j<FragmentsPerOrder[i];j++) {
      Log() << Verbose(0) << " " << OrderSet[i][j];
    }
    Log() << Verbose(0) << endl;
  }
  Log() << Verbose(0) << endl;


  return true;
};

/** Compares each entry in \a *SmallerSet if it is containted in \a *GreaterSet.
 * \param GreaterSet index to greater set
 * \param SmallerSet index to smaller set
 * \return true if all keys of SmallerSet contained in GreaterSet
 */
bool KeySetsContainer::Contains(const int GreaterSet, const int SmallerSet)
{
  bool result = true;
  bool intermediate;
  if ((GreaterSet < 0) || (SmallerSet < 0) || (GreaterSet > FragmentCounter) || (SmallerSet > FragmentCounter)) // index out of bounds
    return false;
  for(int i=AtomCounter[SmallerSet];i--;) {
    intermediate = false;
    for (int j=AtomCounter[GreaterSet];j--;)
      intermediate = (intermediate || ((KeySets[SmallerSet][i] == KeySets[GreaterSet][j]) || (KeySets[SmallerSet][i] == -1)));
    result = result && intermediate;
  }

  return result;
};


// ======================================= END =============================================