Index: src/parser.cpp
===================================================================
--- src/parser.cpp	(revision f731aed8e27c7b0f1fcfad0872e0742f3403b3b8)
+++ src/parser.cpp	(revision eeec8fb96069e10986806be23efb05adf344d003)
@@ -56,8 +56,9 @@
 MatrixContainer::MatrixContainer() {
   Indices = NULL;
-  Header = (char *) Malloc(sizeof(char)*1023, "MatrixContainer::MatrixContainer: *Header");
+  Header = (char **) Malloc(sizeof(char)*1, "MatrixContainer::MatrixContainer: **Header");
   Matrix = (double ***) Malloc(sizeof(double **)*(1), "MatrixContainer::MatrixContainer: ***Matrix"); // one more each for the total molecule
   RowCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *RowCounter");
   ColumnCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *ColumnCounter");
+  Header[0] = NULL;
   Matrix[0] = NULL;
   RowCounter[0] = -1;
@@ -90,5 +91,8 @@
   Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices"); 
   
-  Free((void **)&Header, "MatrixContainer::~MatrixContainer: *Header");
+  if (Header != NULL)
+    for(int i=MatrixCounter+1;i--;)
+      Free((void **)&Header[i], "MatrixContainer::~MatrixContainer: *Header[]");
+  Free((void **)&Header, "MatrixContainer::~MatrixContainer: **Header");
   Free((void **)&RowCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
   Free((void **)&ColumnCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
@@ -120,6 +124,9 @@
         return false;
       Indices[i] = (int *) Malloc(sizeof(int)*Matrix->RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
-      for(int j=Matrix->RowCounter[i];j--;)
+      for(int j=Matrix->RowCounter[i];j--;) {
         Indices[i][j] = Matrix->Indices[i][j];
+        //cout << Indices[i][j] << "\t";
+      }
+      //cout << endl;
     }
   }
@@ -153,13 +160,14 @@
     return false;
   }
-  
-  // skip some initial lines
+
+  // parse header
+  Header[MatrixNr] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseMatrix: *Header[]");
   for (int m=skiplines+1;m--;)
-    input.getline(Header, 1023);
+    input.getline(Header[MatrixNr], 1023);
   
   // scan header for number of columns
-  line.str(Header);
+  line.str(Header[MatrixNr]);
   for(int k=skipcolumns;k--;)
-    line >> Header;
+    line >> Header[MatrixNr];
   //cout << line.str() << endl;
   ColumnCounter[MatrixNr]=0;
@@ -169,5 +177,5 @@
   }
   //cout << line.str() << endl;
-  cout << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
+  //cout << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
   if (ColumnCounter[MatrixNr] == 0)
     cerr << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
@@ -183,5 +191,5 @@
     }
   }
-  cout << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
+  //cout << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
   if (RowCounter[MatrixNr] == 0)
     cerr << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
@@ -189,5 +197,5 @@
   // allocate matrix if it's not zero dimension in one direction
   if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
-    Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
+    Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseMatrix: **Matrix[]");
   
     // parse in each entry for this matrix
@@ -195,11 +203,11 @@
     input.seekg(ios::beg);
     for (int m=skiplines+1;m--;)
-      input.getline(Header, 1023);    // skip header
-    line.str(Header);
+      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, line.str().c_str(), 1023);  
+    strncpy(Header[MatrixNr], line.str().c_str(), 1023);  
     for(int j=0;j<RowCounter[MatrixNr];j++) {
-      Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+      Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[][]");
       input.getline(filename, 1023);
       stringstream lines(filename);
@@ -212,5 +220,5 @@
       }
       //cout << endl;
-      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseFragmentMatrix: *Matrix[RowCounter[MatrixNr]][]");
+      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[RowCounter[MatrixNr]][]");
       for(int j=ColumnCounter[MatrixNr];j--;)
         Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
@@ -266,4 +274,5 @@
   
   cout << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl;
+  Header = (char **) ReAlloc(Header, sizeof(char *)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: **Header"); // one more each for the total molecule
   Matrix = (double ***) ReAlloc(Matrix, sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
   RowCounter = (int *) ReAlloc(RowCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
@@ -271,4 +280,5 @@
   for(int i=MatrixCounter+1;i--;) {
     Matrix[i] = NULL;
+    Header[i] = NULL;
     RowCounter[i] = -1;
     ColumnCounter[i] = -1;
@@ -287,25 +297,25 @@
 
 /** Allocates and resets the memory for a number \a MCounter of matrices.
- * \param *GivenHeader Header line
+ * \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 matrices
+ * \param *CCounter number of columns for each matrix
  * \return Allocation successful 
  */
-bool MatrixContainer::AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int *CCounter)
-{
-  Header = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseFragmentMatrix: *EnergyHeader");
-  strncpy(Header, GivenHeader, 1023);
-
+bool MatrixContainer::AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter)
+{
   MatrixCounter = MCounter;
-  Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
-  RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
-  ColumnCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *ColumnCounter");
+  Header = (char **) Malloc(sizeof(char *)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *Header");
+  Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: ***Matrix"); // one more each for the total molecule
+  RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *RowCounter");
+  ColumnCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *ColumnCounter");
   for(int i=MatrixCounter+1;i--;) {
+    Header[i] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::AllocateMatrix: *Header[i]");
+    strncpy(Header[i], GivenHeader[i], 1023);
     RowCounter[i] = RCounter[i];
     ColumnCounter[i] = CCounter[i];
-    Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");  
+    Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::AllocateMatrix: **Matrix[]");  
     for(int j=RowCounter[i]+1;j--;) {
-      Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter[i], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+      Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter[i], "MatrixContainer::AllocateMatrix: *Matrix[][]");
       for(int k=ColumnCounter[i];k--;)
         Matrix[i][j][k] = 0.;
@@ -385,5 +395,5 @@
 };
 
-/** Sums the energy with each factor and put into last element of \a ***Energies.
+/** 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.
@@ -462,5 +472,5 @@
       return false;
     }
-    output << Header << endl;
+    output << Header[i] << endl;
     for(int j=0;j<RowCounter[i];j++) {
       for(int k=0;k<ColumnCounter[i];k++)
@@ -491,5 +501,5 @@
     return false;
   }
-  output << Header << endl;
+  output << Header[MatrixCounter] << endl;
   for(int j=0;j<RowCounter[MatrixCounter];j++) {
     for(int k=0;k<ColumnCounter[MatrixCounter];k++)
@@ -765,5 +775,5 @@
       int j = Indices[ FragmentNr ][l];
       if (j > RowCounter[MatrixCounter]) {
-        cerr << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
+        cerr << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
         return false;
       }
@@ -772,5 +782,5 @@
           int k = Indices[ FragmentNr ][m];
           if (k > ColumnCounter[MatrixCounter]) {
-            cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << "!" << endl;
+            cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
             return false;
           }
@@ -785,4 +795,83 @@
 };
 
+/** 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 KeySet 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 &KeySet, int Order)
+{
+  // go through each order
+  for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
+    //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.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<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
+        if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
+          //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
+          // if the fragment's indices are all in the current fragment
+          for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
+            int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
+            //cout << "Current row index is " << k << "/" << m << "." << endl;
+            if (m != -1) { // if it's not an added hydrogen
+              for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
+                //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
+                if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
+                  m = l;
+                  break;  
+                }
+              }
+              //cout << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl;
+              if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
+                cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
+                return false;
+              }
+              
+              for(int l=0;l<ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l++) { 
+                int n = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][l];
+                //cout << "Current column index is " << l << "/" << n << "." << endl;
+                if (n != -1) { // if it's not an added hydrogen
+                  for (int p=0;p<ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment
+                    //cout << "Comparing " << n << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p] << "." << endl;
+                    if (n == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p]) {
+                      n = p;
+                      break;  
+                    }
+                  }
+                  //cout << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl;
+                  if (n > ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
+                    cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
+                    return false;
+                  }
+                  if (Order == SubOrder) { // equal order is always copy from Energies
+                    Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
+                  } else {
+                    Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
+                  }
+                }
+              }
+            }
+            //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
+             //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
+          }
+        } else {
+          //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
+        }
+      }
+    }
+   //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.AtomCounter[0]-1][1] << endl;
+  }
+  
+  return true;
+};
 
 /** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
