Index: src/joiner.cpp
===================================================================
--- src/joiner.cpp	(revision 450d63e3aecf24a3631d6e9db7a79102c715e4e8)
+++ src/joiner.cpp	(revision 85d278634d82c0d87206785b8cb85c4c6f794e06)
@@ -19,9 +19,15 @@
   periodentafel *periode = NULL; // and a period table of all elements
   EnergyMatrix Energy;
+  EnergyMatrix EnergyFragments;
+  
   EnergyMatrix Hcorrection;
+  EnergyMatrix HcorrectionFragments;
+  
   ForceMatrix Force;
-  EnergyMatrix EnergyFragments;
-  EnergyMatrix HcorrectionFragments;
   ForceMatrix ForceFragments;
+
+  HessianMatrix Hessian;
+  HessianMatrix HessianFragments;
+  
   ForceMatrix Shielding;
   ForceMatrix ShieldingPAS;
@@ -32,4 +38,5 @@
   char *dir = NULL;
   bool Hcorrected = true;
+  bool NoHessian = false;
 
   cout << "Joiner" << endl;
@@ -63,4 +70,8 @@
   Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX, 0,0);
   if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix, 0,0)) return 1;
+  if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix, 0,0)) {
+    NoHessian = true;
+    cout << "No hessian matrices found, skipping those.";
+  }
   if (periode != NULL) { // also look for PAS values
     if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
@@ -74,4 +85,7 @@
   // ---------- Parse the Force indices into an array ---------------
   if (!Force.ParseIndices(argv[1])) return 1;
+
+  // ---------- Parse the Hessian (=force) indices into an array ---------------
+  if (!Hessian.ParseIndices(argv[1])) return 1;
 
   // ---------- Parse the shielding indices into an array ---------------
@@ -88,4 +102,6 @@
   if (Hcorrected)  HcorrectionFragments.AllocateMatrix(Hcorrection.Header, Hcorrection.MatrixCounter, Hcorrection.RowCounter, Hcorrection.ColumnCounter);
   if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
+  if (!NoHessian)
+    if (!HessianFragments.AllocateMatrix(Hessian.Header, Hessian.MatrixCounter, Hessian.RowCounter, Hessian.ColumnCounter)) return 1;
   if (periode != NULL) { // also look for PAS values
     if (!ShieldingFragments.AllocateMatrix(Shielding.Header, Shielding.MatrixCounter, Shielding.RowCounter, Shielding.ColumnCounter)) return 1;
@@ -96,4 +112,6 @@
   if(!Energy.SetLastMatrix(0., 0)) return 1;
   if(!Force.SetLastMatrix(0., 2)) return 1;
+  if (!NoHessian)
+    if (!Hessian.SetLastMatrix(0., 0)) return 1;
   if (periode != NULL) { // also look for PAS values
     if(!Shielding.SetLastMatrix(0., 2)) return 1;
@@ -118,4 +136,10 @@
     if (!ForceFragments.SumSubManyBodyTerms(Force, KeySet, BondOrder)) return 1;
     if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return 1;
+    // --------- sum up Hessian --------------------
+    if (!NoHessian) {
+      cout << "Summing Hessian of order " << BondOrder+1 << " ..." << endl;
+      if (!HessianFragments.SumSubManyBodyTerms(Hessian, KeySet, BondOrder)) return 1;
+      if (!Hessian.SumSubHessians(HessianFragments, KeySet, BondOrder, 1.)) return 1;
+    }
     if (periode != NULL) { // also look for PAS values
       cout << "Summing shieldings of order " << BondOrder+1 << " ..." << endl;
@@ -134,4 +158,6 @@
     // forces
     if (!Force.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ForcesSuffix)) return 1;
+    // hessian
+    if (!Hessian.WriteLastMatrix(argv[1], (prefix.str()).c_str(), HessianSuffix)) return 1;
     // shieldings
     if (periode != NULL) { // also look for PAS values
@@ -153,4 +179,9 @@
   if (!ForceFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
   if (!CreateDataFragment(EnergyFragments, KeySet, argv[1], FRAGMENTPREFIX ENERGYPERFRAGMENT, "fragment energy versus the Fragment No", "today", CreateEnergy)) return 1;
+  if (!NoHessian) {
+    prefix.str(" ");
+    prefix << dir << HessianFragmentSuffix;
+    if (!HessianFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+  }
   if (periode != NULL) { // also look for PAS values
     prefix.str(" ");
@@ -166,4 +197,6 @@
   if (Hcorrected) Hcorrection.WriteLastMatrix(argv[1], dir, HcorrectionFragmentSuffix);
   if (!Force.WriteLastMatrix(argv[1], dir, ForceFragmentSuffix)) return 1;
+  if (!NoHessian)
+    if (!Hessian.WriteLastMatrix(argv[1], dir, HessianFragmentSuffix)) return 1;
   if (periode != NULL) { // also look for PAS values
     if (!Shielding.WriteLastMatrix(argv[1], dir, ShieldingFragmentSuffix)) return 1;
Index: src/parser.cpp
===================================================================
--- src/parser.cpp	(revision 450d63e3aecf24a3631d6e9db7a79102c715e4e8)
+++ src/parser.cpp	(revision 85d278634d82c0d87206785b8cb85c4c6f794e06)
@@ -673,4 +673,145 @@
 };
 
+// ======================================= 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;
+  
+  cout << "Parsing force indices for " << MatrixCounter << " matrices." << endl;
+  Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "HessianMatrix::ParseIndices: **Indices");
+  line << name << FRAGMENTPREFIX << FORCESFILE;
+  input.open(line.str().c_str(), ios::in);
+  //cout << "Opening " << line.str() << " ... "  << input << endl;
+  if (input == NULL) {
+    cout << 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] = (int *) Malloc(sizeof(int)*RowCounter[i], "HessianMatrix::ParseIndices: *Indices[]");
+    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
+    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
+    Free((void **)&FragmentNumber, "HessianMatrix::ParseIndices: *FragmentNumber");
+    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
+      line >> Indices[i][j];
+      //cout << " " << Indices[i][j];
+    }
+    //cout << endl;
+  }
+  Indices[MatrixCounter] = (int *) Malloc(sizeof(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 KeySet 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 &KeySet, int Order, double sign)
+{
+  int FragmentNr;
+  // sum forces
+  for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++) {
+    FragmentNr = KeySet.OrderSet[Order][i];
+    for(int l=0;l<RowCounter[ FragmentNr ];l++) {
+      int j = Indices[ FragmentNr ][l];
+      if (j > RowCounter[MatrixCounter]) {
+        cerr << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
+        return false;
+      }
+      if (j != -1) {
+        //for(int m=0;m<ColumnCounter[ FragmentNr ];m++) {
+        for(int m=0;m<ColumnCounter;m++) {
+          int k = Indices[ FragmentNr ][m];
+//          if (k > ColumnCounter[MatrixCounter]) {
+//            cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << "!" << endl;
+          if (k > ColumnCounter) {
+            cerr << "Current hessian index " << k << " is greater than " << ColumnCounter << "!" << endl;
+            return false;
+          }
+          if (k != -1) {
+            Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][m];
+          }
+        }
+      }
+    }
+  }
+  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(char *name, 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) {
+      cout << 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;
+        //cout << "Current index: " << nr << "." << endl;
+        if (nr > RowCounter[MatrixCounter])
+          RowCounter[MatrixCounter] = nr;
+      }
+    }
+    RowCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
+    input.close();
+  
+    // allocate last plus one matrix
+    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter << " columns." << endl;
+    Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
+    for(int j=0;j<=RowCounter[MatrixCounter];j++)
+      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "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 =============================
 
Index: src/parser.hpp
===================================================================
--- src/parser.hpp	(revision 450d63e3aecf24a3631d6e9db7a79102c715e4e8)
+++ src/parser.hpp	(revision 85d278634d82c0d87206785b8cb85c4c6f794e06)
@@ -19,6 +19,12 @@
 
 #define EnergySuffix ".energy.all"
+#define EnergyFragmentSuffix ".energyfragment.all"
+#define ForcesSuffix ".forces.all"
+#define ForceFragmentSuffix ".forcefragment.all"
 #define HcorrectionSuffix ".Hcorrection.all"
-#define ForcesSuffix ".forces.all"
+#define HcorrectionFragmentSuffix ".Hcorrectionfragment.all"
+#define HessianSuffix ".hessian_xx.all"
+#define HessianFragmentSuffix ".hessianfragment_xx.all"
+#define OrderSuffix ".Order"
 #define ShieldingSuffix ".sigma_all.csv"
 #define ShieldingPASSuffix ".sigma_all_PAS.csv"
@@ -26,8 +32,4 @@
 #define ShieldingPASFragmentSuffix ".sigma_all_PAS_fragment.all"
 #define TimeSuffix ".speed"
-#define EnergyFragmentSuffix ".energyfragment.all"
-#define HcorrectionFragmentSuffix ".Hcorrectionfragment.all"
-#define ForceFragmentSuffix ".forcefragment.all"
-#define OrderSuffix ".Order"
 
 // ======================================= FUNCTIONS ==========================================
@@ -84,4 +86,13 @@
 };
 
+// ======================================= CLASS HessianMatrix =============================
+
+class HessianMatrix : public MatrixContainer {
+  public: 
+    bool ParseIndices(char *name);
+    bool SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign);
+    bool ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns);
+};
+
 // ======================================= CLASS KeySetsContainer =============================
 
