Index: src/Fragmentation/parser.cpp
===================================================================
--- src/Fragmentation/parser.cpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
+++ src/Fragmentation/parser.cpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
@@ -0,0 +1,1214 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
+/** \file parsing.cpp
+ *
+ * Declarations of various class functions for the parsing of value files.
+ *
+ */
+
+// ======================================= INCLUDES ==========================================
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include <cstring>
+#include <fstream>
+#include <iomanip>
+
+#include "parser.hpp"
+#include "CodePatterns/Log.hpp"
+#include "CodePatterns/Verbose.hpp"
+#include "Helpers/defs.hpp"
+#include "Helpers/helpers.hpp"
+
+using namespace std;
+
+// ======================================= 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.fail()) {
+    if (!test)
+      DoLog(0) && (Log() << Verbose(0) << endl << "FilePresent: 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);
+};
+
+/** Returns a string with \a i prefixed with 0s to match order of total number of molecules in digits.
+ * \param FragmentNumber total number of fragments to determine necessary number of digits
+ * \param digits number to create with 0 prefixed
+ * \return allocated(!) char array with number in digits, ten base.
+ */
+inline char *FixedDigitNumber(const int FragmentNumber, const int digits)
+{
+  char *returnstring;
+  int number = FragmentNumber;
+  int order = 0;
+  while (number != 0) { // determine number of digits needed
+    number = (int)floor(((double)number / 10.));
+    order++;
+    //Log() << Verbose(0) << "Number is " << number << ", order is " << order << "." << endl;
+  }
+  // allocate string
+  returnstring = new char[order + 2];
+  // terminate  and fill string array from end backward
+  returnstring[order] = '\0';
+  number = digits;
+  for (int i=order;i--;) {
+    returnstring[i] = '0' + (char)(number % 10);
+    number = (int)floor(((double)number / 10.));
+  }
+  //Log() << Verbose(0) << returnstring << endl;
+  return returnstring;
+};
+
+// ======================================= CLASS MatrixContainer =============================
+
+/** Constructor of MatrixContainer class.
+ */
+MatrixContainer::MatrixContainer() :
+  Indices(NULL)
+{
+  Header = new char*[1];
+  Matrix = new double**[1]; // one more each for the total molecule
+  RowCounter = new int[1];
+  ColumnCounter = new int[1];
+  Header[0] = NULL;
+  Matrix[0] = NULL;
+  RowCounter[0] = -1;
+  MatrixCounter = 0;
+  ColumnCounter[0] = -1;
+};
+
+/** Destructor of MatrixContainer class.
+ */
+MatrixContainer::~MatrixContainer() {
+  if (Matrix != NULL) {
+    // free Matrix[MatrixCounter]
+    if ((ColumnCounter != NULL) && (RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
+      for(int j=RowCounter[MatrixCounter]+1;j--;)
+        delete[](Matrix[MatrixCounter][j]);
+    //if (MatrixCounter != 0)
+    delete[](Matrix[MatrixCounter]);
+    // free all matrices but ultimate one
+    for(int i=MatrixCounter;i--;) {
+      if ((ColumnCounter != NULL) && (RowCounter != NULL)) {
+          for(int j=RowCounter[i]+1;j--;)
+            delete[](Matrix[i][j]);
+          delete[](Matrix[i]);
+      }
+    }
+    delete[](Matrix);
+  }
+  // free indices
+  if (Indices != NULL)
+    for(int i=MatrixCounter+1;i--;) {
+      delete[](Indices[i]);
+    }
+  delete[](Indices);
+  
+  // free header and counters
+  if (Header != NULL)
+    for(int i=MatrixCounter+1;i--;)
+      delete[](Header[i]);
+  delete[](Header);
+  delete[](RowCounter);
+  delete[](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) 
+{
+  DoLog(0) && (Log() << Verbose(0) << "Initialising indices");
+  if (Matrix == NULL) {
+    DoLog(0) && (Log() << Verbose(0) << " with trivial mapping." << endl);
+    Indices = new int*[MatrixCounter + 1];
+    for(int i=MatrixCounter+1;i--;) {
+      Indices[i] = new int[RowCounter[i]];
+      for(int j=RowCounter[i];j--;)
+        Indices[i][j] = j;
+    }
+  } else {
+    DoLog(0) && (Log() << Verbose(0) << " from other MatrixContainer." << endl);
+    if (MatrixCounter != Matrix->MatrixCounter)
+      return false;
+    Indices = new int*[MatrixCounter + 1];
+    for(int i=MatrixCounter+1;i--;) {
+      if (RowCounter[i] != Matrix->RowCounter[i])
+        return false;
+      Indices[i] = new int[Matrix->RowCounter[i]];
+      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 &input input stream
+ * \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(std::istream &input, int skiplines, int skipcolumns, int MatrixNr)
+{
+  stringstream line;
+  string token;
+  char filename[1023];
+
+  if (input.fail()) {
+    DoeLog(1) && (eLog()<< Verbose(1) << endl << "MatrixContainer::ParseMatrix: Unable to parse istream." << endl);
+    //performCriticalExit();
+    return false;
+  }
+
+  // parse header
+  if (Header[MatrixNr] != NULL)
+    delete[] Header[MatrixNr];
+  Header[MatrixNr] = new char[1024];
+  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(1) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
+  if (ColumnCounter[MatrixNr] == 0) {
+    DoeLog(0) && (eLog()<< Verbose(0) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from ostream." << 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(1) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from input stream." << endl;
+  if (RowCounter[MatrixNr] == 0) {
+    DoeLog(0) && (eLog()<< Verbose(0) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from input stream." << endl);
+    performCriticalExit();
+  }
+
+  // allocate matrix if it's not zero dimension in one direction
+  if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
+    if (Matrix[MatrixNr] != NULL)
+      delete[] Matrix[MatrixNr];
+    Matrix[MatrixNr] = new double*[RowCounter[MatrixNr] + 1];
+    for(int j=0;j<RowCounter[MatrixNr]+1;j++)
+      Matrix[MatrixNr][j] = 0;
+
+    // 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]);
+    Log() << Verbose(0) << "Header: " << line.str() << endl;
+    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++) {
+      if (Matrix[MatrixNr][j] != NULL)
+        delete[] Matrix[MatrixNr][j];
+      Matrix[MatrixNr][j] = new double[ColumnCounter[MatrixNr]];
+      for(int k=0;k<ColumnCounter[MatrixNr];k++)
+        Matrix[MatrixNr][j][k] = 0;
+
+      input.getline(filename, 1023);
+      stringstream lines(filename);
+      Log() << Verbose(2) << "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(1) << " " << setprecision(2) << Matrix[MatrixNr][j][k] << endl;
+      }
+      if (Matrix[MatrixNr][ RowCounter[MatrixNr] ] != NULL)
+        delete[] Matrix[MatrixNr][ RowCounter[MatrixNr] ];
+      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = new double[ColumnCounter[MatrixNr]];
+      for(int j=ColumnCounter[MatrixNr];j--;)
+        Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
+    }
+  } else {
+    DoeLog(1) && (eLog()<< Verbose(1) << "Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl);
+  }
+  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.bad()) {
+    DoLog(0) && (Log() << Verbose(0) << endl << "MatrixContainer::ParseFragmentMatrix: 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();
+  DoLog(0) && (Log() << Verbose(0) << "Determined " << MatrixCounter << " fragments." << endl);
+
+  DoLog(0) && (Log() << Verbose(0) << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl);
+  delete[](Header);
+  delete[](Matrix);
+  delete[](RowCounter);
+  delete[](ColumnCounter);
+  Header = new char*[MatrixCounter + 1]; // one more each for the total molecule
+  Matrix = new double**[MatrixCounter + 1]; // one more each for the total molecule
+  RowCounter = new int[MatrixCounter + 1];
+  ColumnCounter = new int[MatrixCounter + 1];
+  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;
+    std::ifstream input(file.str().c_str());
+    DoLog(0) &&( Log() << Verbose(0) << "Opening " << file.str() << " ... " << endl);
+    if (!ParseMatrix(input, skiplines, skipcolumns, i)) {
+      input.close();
+      return false;
+    }
+    input.close();
+    delete[](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 = new char*[MatrixCounter + 1];
+  Matrix = new double**[MatrixCounter + 1]; // one more each for the total molecule
+  RowCounter = new int[MatrixCounter + 1];
+  ColumnCounter = new int[MatrixCounter + 1];
+  for(int i=MatrixCounter+1;i--;) {
+    Header[i] = new char[1024];
+    strncpy(Header[i], GivenHeader[i], 1023);
+    RowCounter[i] = RCounter[i];
+    ColumnCounter[i] = CCounter[i];
+    Matrix[i] = new double*[RowCounter[i] + 1];
+    for(int j=RowCounter[i]+1;j--;) {
+      Matrix[i][j] = new double[ColumnCounter[i]];
+      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] ]) {
+                DoeLog(0) && (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;
+
+  DoLog(0) && (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;
+    delete[](FragmentNumber);
+    output.open(line.str().c_str(), ios::out);
+    if (output == NULL) {
+      DoeLog(0) && (eLog()<< Verbose(0) << "MatrixContainer::WriteTotalFragments: 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;
+
+  DoLog(0) && (Log() << Verbose(0) << "Writing matrix values of " << suffix << "." << endl);
+  line << name << prefix << suffix;
+  output.open(line.str().c_str(), ios::out);
+  if (output == NULL) {
+    DoeLog(0) && (eLog()<< Verbose(0) << "MatrixContainer::WriteLastMatrix: 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;
+};
+
+std::ostream & operator << (std::ostream &ost, const MatrixContainer &m)
+{
+  for (int i=0;i<=m.MatrixCounter;++i) {
+    ost << "Printing matrix " << i << ":" << std::endl;
+    for (int j=0;j<=m.RowCounter[i];++j) {
+      for (int k=0;k<m.ColumnCounter[i];++k) {
+        ost << m.Matrix[i][j][k] << " ";
+      }
+      ost << std::endl;
+    }
+  }
+  ost << std::endl;
+  return ost;
+}
+
+// ======================================= 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()
+{
+  DoLog(0) && (Log() << Verbose(0) << "Parsing energy indices." << endl);
+  Indices = new int*[MatrixCounter + 1];
+  for(int i=MatrixCounter+1;i--;) {
+    Indices[i] = new int[RowCounter[i]];
+    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
+    DoLog(0) && (Log() << Verbose(0) << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl);
+    Matrix[MatrixCounter] = new double*[RowCounter[MatrixCounter] + 1];
+    for(int j=0;j<=RowCounter[MatrixCounter];j++)
+      Matrix[MatrixCounter][j] = new double[ColumnCounter[MatrixCounter]];
+    
+    // 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));
+    std::ifstream input(filename);
+    ParseMatrix(input, skiplines, skipcolumns, MatrixCounter);
+    input.close();
+  }
+  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;
+
+  DoLog(0) && (Log() << Verbose(0) << "Parsing force indices for " << MatrixCounter << " matrices." << endl);
+  Indices = new int*[MatrixCounter + 1];
+  line << name << FRAGMENTPREFIX << FORCESFILE;
+  input.open(line.str().c_str(), ios::in);
+  //Log() << Verbose(0) << "Opening " << line.str() << " ... "  << input << endl;
+  if (input.fail()) {
+    DoLog(0) && (Log() << Verbose(0) << endl << "ForceMatrix::ParseIndices: 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] = new int[RowCounter[i]];
+    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
+    //Log() << Verbose(0) << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
+    delete[](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] = new int[RowCounter[MatrixCounter]];
+  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]) {
+        DoeLog(0) && (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.fail()) {
+      DoLog(0) && (Log() << Verbose(0) << endl << "ForceMatrix::ParseFragmentMatrix: 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: " << getNr() << "." << 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
+    DoLog(0) && (Log() << Verbose(0) << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl);
+    Matrix[MatrixCounter] = new double*[RowCounter[MatrixCounter] + 1];
+    for(int j=0;j<=RowCounter[MatrixCounter];j++)
+      Matrix[MatrixCounter][j] = new double[ColumnCounter[MatrixCounter]];
+
+    // 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));
+    std::ifstream input(filename);
+    ParseMatrix(input, skiplines, skipcolumns, MatrixCounter);
+    input.close();
+  }
+ 
+
+  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;
+  
+  DoLog(0) && (Log() << Verbose(0) << "Parsing hessian indices for " << MatrixCounter << " matrices." << endl);
+  Indices = new int*[MatrixCounter + 1];
+  line << name << FRAGMENTPREFIX << FORCESFILE;
+  input.open(line.str().c_str(), ios::in);
+  //Log() << Verbose(0) << "Opening " << line.str() << " ... "  << input << endl;
+  if (input.fail()) {
+    DoLog(0) && (Log() << Verbose(0) << endl << "HessianMatrix::ParseIndices: 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] = new int[RowCounter[i]];
+    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
+    //Log() << Verbose(0) << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
+    delete[](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] = new int[RowCounter[MatrixCounter]];
+  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]) {
+        DoeLog(0) && (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]) {
+            DoeLog(0) && (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] ]) {
+                DoeLog(0) && (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] ]) {
+                    DoeLog(0) && (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.fail()) {
+      DoLog(0) && (Log() << Verbose(0) << endl << "HessianMatrix::ParseFragmentMatrix: 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: " << getNr() << "." << 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
+    DoLog(0) && (Log() << Verbose(0) << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl);
+    Matrix[MatrixCounter] = new double*[RowCounter[MatrixCounter] + 1];
+    for(int j=0;j<=RowCounter[MatrixCounter];j++)
+      Matrix[MatrixCounter][j] = new double[ColumnCounter[MatrixCounter]];
+
+    // 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));
+    std::ifstream input(filename);
+    ParseMatrix(input, skiplines, skipcolumns, MatrixCounter);
+    input.close();
+  }
+
+
+  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--;)
+    delete[](KeySets[i]);
+  for(int i=Order;i--;)
+    delete[](OrderSet[i]);
+  delete[](KeySets);
+  delete[](OrderSet);
+  delete[](AtomCounter);
+  delete[](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;
+  DoLog(0) && (Log() << Verbose(0) << "Parsing key sets." << endl);
+  KeySets = new int*[FragmentCounter];
+  for(int i=FragmentCounter;i--;)
+    KeySets[i] = NULL;
+  file << name << FRAGMENTPREFIX << KEYSETFILE;
+  input.open(file.str().c_str(), ios::in);
+  if (input.fail()) {
+    DoLog(0) && (Log() << Verbose(0) << endl << "KeySetsContainer::ParseKeySets: Unable to open " << file.str() << ", is the directory correct?" << endl);
+    return false;
+  }
+
+  AtomCounter = new int[FragmentCounter];
+  for(int i=0;(i<FragmentCounter) && (!input.eof());i++) {
+    stringstream line;
+    AtomCounter[i] = ACounter[i];
+    // parse the values
+    KeySets[i] = new int[AtomCounter[i]];
+    for(int j=AtomCounter[i];j--;)
+      KeySets[i][j] = -1;
+    FragmentNumber = FixedDigitNumber(FragmentCounter, i);
+    //Log() << Verbose(0) << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
+    delete[](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;
+
+  DoLog(0) && (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;
+  }
+  DoLog(0) && (Log() << Verbose(0) << "Found Order is " << Order << "." << endl);
+
+  // scan through all to determine fragments per order
+  FragmentsPerOrder = new int[Order];
+  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++)
+    DoLog(0) && (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 = new int*[Order];
+  for(int i=Order;i--;)
+    OrderSet[i] = new int[FragmentsPerOrder[i]];
+  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]++;
+  }
+  DoLog(0) && (Log() << Verbose(0) << "Printing OrderSet." << endl);
+  for(int i=0;i<Order;i++) {
+    for (int j=0;j<FragmentsPerOrder[i];j++) {
+      DoLog(0) && (Log() << Verbose(0) << " " << OrderSet[i][j]);
+    }
+    DoLog(0) && (Log() << Verbose(0) << endl);
+  }
+  DoLog(0) && (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 =============================================
