Changes in molecuilder/src/parser.cpp [25f5089:848729]

File:

• molecuilder/src/parser.cpp (modified) (24 diffs, 1 prop)

molecuilder/src/parser.cpp

@@ -56 +56 @@
 MatrixContainer::MatrixContainer() {
   Indices = NULL;
-  Header = (char **) Malloc(sizeof(char)*1, "MatrixContainer::MatrixContainer: **Header");
+  Header = (char *) Malloc(sizeof(char)*1023, "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;
   MatrixCounter = 0;
-  ColumnCounter[0] = -1;
+  ColumnCounter = 0;
 };
 
@@ -72 +70 @@
   if (Matrix != NULL) {
     for(int i=MatrixCounter;i--;) {
-      if ((ColumnCounter != NULL) && (RowCounter != NULL)) {
+      if (RowCounter != NULL) {
           for(int j=RowCounter[i]+1;j--;)
             Free((void **)&Matrix[i][j], "MatrixContainer::~MatrixContainer: *Matrix[][]");
@@ -78 +76 @@
       }
     }
-    if ((ColumnCounter != NULL) && (RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
+    if ((RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
       for(int j=RowCounter[MatrixCounter]+1;j--;)
         Free((void **)&Matrix[MatrixCounter][j], "MatrixContainer::~MatrixContainer: *Matrix[MatrixCounter][]");
@@ -91 +89 @@
   Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices"); 
   
-  if (Header != NULL)
-    for(int i=MatrixCounter+1;i--;)
-      Free((void **)&Header[i], "MatrixContainer::~MatrixContainer: *Header[]");
-  Free((void **)&Header, "MatrixContainer::~MatrixContainer: **Header");
+  Free((void **)&Header, "MatrixContainer::~MatrixContainer: *Header");
   Free((void **)&RowCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
-  Free((void **)&ColumnCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
-};
-
-/** 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) 
-{
-  cout << "Initialising indices";
-  if (Matrix == NULL) {
-    cout << " with trivial mapping." << endl;
-    Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "MatrixContainer::InitialiseIndices: **Indices");
-    for(int i=MatrixCounter+1;i--;) {
-      Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
-      for(int j=RowCounter[i];j--;)
-        Indices[i][j] = j;
-    }
-  } else {
-    cout << " from other MatrixContainer." << endl;
-    if (MatrixCounter != Matrix->MatrixCounter)
-      return false;
-    Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "MatrixContainer::InitialiseIndices: **Indices");
-    for(int i=MatrixCounter+1;i--;) {
-      if (RowCounter[i] != Matrix->RowCounter[i])
-        return false;
-      Indices[i] = (int *) Malloc(sizeof(int)*Matrix->RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
-      for(int j=Matrix->RowCounter[i];j--;) {
-        Indices[i][j] = Matrix->Indices[i][j];
-        //cout << Indices[i][j] << "\t";
-      }
-      //cout << endl;
-    }
-  }
-  return true; 
-};
+};
+
 
 /** Parsing a number of matrices.
@@ -160 +120 @@
     return false;
   }
-
-  // parse header
-  Header[MatrixNr] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseMatrix: *Header[]");
+  
+  // skip some initial lines
   for (int m=skiplines+1;m--;)
-    input.getline(Header[MatrixNr], 1023);
+    input.getline(Header, 1023);
   
   // scan header for number of columns
-  line.str(Header[MatrixNr]);
+  line.str(Header);
   for(int k=skipcolumns;k--;)
-    line >> Header[MatrixNr];
+    line >> Header;
   //cout << line.str() << endl;
-  ColumnCounter[MatrixNr]=0;
+  ColumnCounter=0;
   while ( getline(line,token, '\t') ) {
     if (token.length() > 0)
-      ColumnCounter[MatrixNr]++;
+      ColumnCounter++;
   }
   //cout << line.str() << 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;
+  //cout << "ColumnCounter: " << ColumnCounter << "." << endl;
+  if (ColumnCounter == 0)
+    cerr << "ColumnCounter: " << ColumnCounter << " from file " << name << ", this is probably an error!" << endl;
   
   // scan rest for number of rows/lines
@@ -196 +155 @@
   
   // 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::ParseMatrix: **Matrix[]");
+  if ((ColumnCounter > 0) && (RowCounter[MatrixNr] > -1)) {
+    Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
   
     // parse in each entry for this matrix
@@ -203 +162 @@
     input.seekg(ios::beg);
     for (int m=skiplines+1;m--;)
-      input.getline(Header[MatrixNr], 1023);    // skip header
-    line.str(Header[MatrixNr]);
+      input.getline(Header, 1023);    // skip header
+    line.str(Header);
     for(int k=skipcolumns;k--;)  // skip columns in header too
       line >> filename;
-    strncpy(Header[MatrixNr], line.str().c_str(), 1023);  
+    strncpy(Header, line.str().c_str(), 1023);  
     for(int j=0;j<RowCounter[MatrixNr];j++) {
-      Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[][]");
+      Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
       input.getline(filename, 1023);
       stringstream lines(filename);
@@ -215 +174 @@
       for(int k=skipcolumns;k--;)
         lines >> filename;
-      for(int k=0;(k<ColumnCounter[MatrixNr]) && (!lines.eof());k++) {
+      for(int k=0;(k<ColumnCounter) && (!lines.eof());k++) {
         lines >> Matrix[MatrixNr][j][k];
         //cout << " " << setprecision(2) << Matrix[MatrixNr][j][k];;
       }
       //cout << endl;
-      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] ] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[RowCounter[MatrixNr]][]");
+      for(int j=ColumnCounter;j--;)
         Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
     }
   } else {
-    cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl; 
+    cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl; 
   }
   input.close();
@@ -274 +233 @@
   
   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");
-  ColumnCounter = (int *) ReAlloc(ColumnCounter, sizeof(int)*(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++) {
@@ -297 +252 @@
 
 /** Allocates and resets the memory for a number \a MCounter of matrices.
- * \param **GivenHeader Header line for each matrix
+ * \param *GivenHeader Header line
  * \param MCounter number of matrices
  * \param *RCounter number of rows for each matrix
- * \param *CCounter number of columns for each matrix
+ * \param CCounter number of columns for all matrices
  * \return Allocation successful 
  */
-bool MatrixContainer::AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter)
-{
+bool MatrixContainer::AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int CCounter)
+{
+  Header = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseFragmentMatrix: *EnergyHeader");
+  strncpy(Header, GivenHeader, 1023);
+
   MatrixCounter = MCounter;
-  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");
+  ColumnCounter = CCounter;
+  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");
   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::AllocateMatrix: **Matrix[]");  
+    Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");  
     for(int j=RowCounter[i]+1;j--;) {
-      Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter[i], "MatrixContainer::AllocateMatrix: *Matrix[][]");
-      for(int k=ColumnCounter[i];k--;)
+      Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+      for(int k=ColumnCounter;k--;)
         Matrix[i][j][k] = 0.;
     }
@@ -332 +286 @@
   for(int i=MatrixCounter+1;i--;)
     for(int j=RowCounter[i]+1;j--;)
-      for(int k=ColumnCounter[i];k--;)
+      for(int k=ColumnCounter;k--;)
         Matrix[i][j][k] = 0.;
    return true;
@@ -345 +299 @@
   for(int i=MatrixCounter+1;i--;)
     for(int j=RowCounter[i]+1;j--;)
-      for(int k=ColumnCounter[i];k--;)
+      for(int k=ColumnCounter;k--;)
         if (fabs(Matrix[i][j][k]) > max)
           max = fabs(Matrix[i][j][k]);
@@ -361 +315 @@
   for(int i=MatrixCounter+1;i--;)
     for(int j=RowCounter[i]+1;j--;)
-      for(int k=ColumnCounter[i];k--;)
+      for(int k=ColumnCounter;k--;)
         if (fabs(Matrix[i][j][k]) < min)
           min = fabs(Matrix[i][j][k]);
@@ -377 +331 @@
 {
   for(int j=RowCounter[MatrixCounter]+1;j--;)
-    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
+    for(int k=skipcolumns;k<ColumnCounter;k++)
       Matrix[MatrixCounter][j][k] = value;
    return true;
@@ -390 +344 @@
 {
   for(int j=RowCounter[MatrixCounter]+1;j--;)
-    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
+    for(int k=skipcolumns;k<ColumnCounter;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 the energy with each factor and put into last element of \a ***Energies.
  * Sums over "E"-terms to create the "F"-terms
- * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
+ * \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
@@ -430 +384 @@
               }
               if (Order == SubOrder) { // equal order is always copy from Energies
-                for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;) // then adds/subtract each column 
+                for(int l=ColumnCounter;l--;) // then adds/subtract each column 
                   Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
               } else {
-                for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;)
+                for(int l=ColumnCounter;l--;)
                   Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
               }
             }
-            //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
+            //if ((ColumnCounter>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;
           }
@@ -472 +426 @@
       return false;
     }
-    output << Header[i] << endl;
+    output << Header << endl;
     for(int j=0;j<RowCounter[i];j++) {
-      for(int k=0;k<ColumnCounter[i];k++)
+      for(int k=0;k<ColumnCounter;k++)
         output << scientific << Matrix[i][j][k] << "\t";
       output << endl;
@@ -501 +455 @@
     return false;
   }
-  output << Header[MatrixCounter] << endl;
+  output << Header << endl;
   for(int j=0;j<RowCounter[MatrixCounter];j++) {
-    for(int k=0;k<ColumnCounter[MatrixCounter];k++)
+    for(int k=0;k<ColumnCounter;k++)
       output << scientific << Matrix[MatrixCounter][j][k] << "\t";
     output << endl;
@@ -513 +467 @@
 // ======================================= 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() 
+{
+  cout << "Parsing energy indices." << endl;
+  Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "EnergyMatrix::ParseIndices: **Indices");
+  for(int i=MatrixCounter+1;i--;) {
+    Indices[i] = (int *) Malloc(sizeof(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 Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
  * \param CorrectionFragments MatrixContainer with hydrogen saturation correction per fragments
  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
@@ -528 +498 @@
     for(int i=KeySet.FragmentsPerOrder[Order];i--;)
       for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
-        for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
+        for(int k=ColumnCounter;k--;)
           Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
   else
     for(int i=KeySet.FragmentsPerOrder[Order];i--;)
       for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
-        for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
+        for(int k=ColumnCounter;k--;)
           Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
   return true;
@@ -554 +524 @@
     // 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)
+    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
-    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
+    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[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     
     // try independently to parse global energysuffix file if present
@@ -623 +589 @@
 
 /** 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 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
@@ -643 +609 @@
       if (j != -1) {
         //if (j == 0) cout << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
-        for(int k=2;k<ColumnCounter[MatrixCounter];k++)
+        for(int k=2;k<ColumnCounter;k++)
           Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
       }
@@ -689 +655 @@
     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
-    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
+    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[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;
-  
-  cout << "Parsing hessian 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] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
-        return false;
-      }
-      if (j != -1) {
-        for(int m=0;m<ColumnCounter[ FragmentNr ];m++) {
-          int k = Indices[ FragmentNr ][m];
-          if (k > ColumnCounter[MatrixCounter]) {
-            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;
-          }
-          if (k != -1) {
-                //cout << "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 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
-                        //cout << "Adding " << MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
-                    Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
-                  } else {
-                        //cout << "Subtracting " << Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
-                    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.
- * \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;
-    ColumnCounter[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;
-          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
-    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " 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[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
 
     // try independently to parse global forcesuffix file if present