Index: src/analyzer.cpp
===================================================================
--- src/analyzer.cpp	(revision 85d278634d82c0d87206785b8cb85c4c6f794e06)
+++ src/analyzer.cpp	(revision f731aed8e27c7b0f1fcfad0872e0742f3403b3b8)
@@ -93,6 +93,6 @@
 
   // ---------- Parse the TE Factors into an array -----------------
-  if (!Energy.ParseIndices()) return 1;
-  if (Hcorrected) Hcorrection.ParseIndices();
+  if (!Energy.InitialiseIndices()) return 1;
+  if (Hcorrected) Hcorrection.InitialiseIndices();
   
   // ---------- Parse the Force indices into an array ---------------
@@ -132,5 +132,5 @@
   output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header << endl;
   for(int j=0;j<Energy.RowCounter[Energy.MatrixCounter];j++) {
-    for(int k=0;k<Energy.ColumnCounter;k++)
+    for(int k=0;k<Energy.ColumnCounter[Energy.MatrixCounter];k++)
       output << scientific << Energy.Matrix[ Energy.MatrixCounter ][j][k] << "\t";
     output << endl;
@@ -140,5 +140,5 @@
   output << endl << "Total Forces" << endl << "===============" << endl << Force.Header << endl;
   for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
-    for(int k=0;k<Force.ColumnCounter;k++)
+    for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
       output << scientific << Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
     output << endl;
@@ -149,5 +149,5 @@
     output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header << endl;
     for(int j=0;j<Shielding.RowCounter[Shielding.MatrixCounter];j++) {
-      for(int k=0;k<Shielding.ColumnCounter;k++)
+      for(int k=0;k<Shielding.ColumnCounter[Shielding.MatrixCounter];k++)
         output << scientific << Shielding.Matrix[ Shielding.MatrixCounter ][j][k] << "\t";
       output << endl;
@@ -157,5 +157,5 @@
     output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header << endl;
     for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
-      for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
+      for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
         output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t";
       output << endl;
@@ -166,5 +166,5 @@
   output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
   for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
-    for(int k=0;k<Time.ColumnCounter;k++) {
+    for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
       output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
     }
@@ -173,5 +173,5 @@
   output << endl;
   output.close();
-  for(int k=0;k<Time.ColumnCounter;k++)
+  for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++)
     Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];
 
@@ -187,5 +187,5 @@
   if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
   for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
-    for(int k=Time.ColumnCounter;k--;) {
+    for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
       Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
     }
@@ -196,5 +196,5 @@
     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
       for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
-        for(int k=Time.ColumnCounter;k--;) {
+        for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
           Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
         }
@@ -202,5 +202,5 @@
     output << BondOrder+1 << "\t" << counter;
     output2 << BondOrder+1 << "\t" << counter;
-    for(int k=0;k<Time.ColumnCounter;k++) {
+    for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
       output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
       if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
@@ -224,5 +224,5 @@
     for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
       output << j << "\t";
-      for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
+      for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
         output << scientific <<  ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
       output << endl;
@@ -256,5 +256,5 @@
   for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
     output << j << "\t";
-    for(int k=0;k<Force.ColumnCounter;k++)
+    for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
       output << scientific <<  Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
     output << endl;
Index: src/datacreator.cpp
===================================================================
--- src/datacreator.cpp	(revision 85d278634d82c0d87206785b8cb85c4c6f794e06)
+++ src/datacreator.cpp	(revision f731aed8e27c7b0f1fcfad0872e0742f3403b3b8)
@@ -67,9 +67,9 @@
     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
       for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
-        for(int k=Fragments.ColumnCounter;k--;)
+        for(int k=Fragments.ColumnCounter[ KeySet.OrderSet[BondOrder][i] ];k--;)
           Fragments.Matrix[Fragments.MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     }
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
-    for (int l=0;l<Fragments.ColumnCounter;l++)
+    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
       output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
     output << endl;
@@ -101,9 +101,9 @@
     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
       for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
-        for(int k=Fragments.ColumnCounter;k--;)
+        for(int k=Fragments.ColumnCounter[ KeySet.OrderSet[BondOrder][i] ];k--;)
           Fragments.Matrix[Fragments.MatrixCounter][j][k] -= Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     }
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
-    for (int l=0;l<Fragments.ColumnCounter;l++)
+    for (int l=0;l<Fragments.ColumnCounter[Energy.MatrixCounter];l++)
       if (fabs(Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]) < MYEPSILON)
         output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
@@ -139,5 +139,5 @@
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
     CreateForce(Fragments, Fragments.MatrixCounter);
-    for (int l=0;l<Fragments.ColumnCounter;l++)
+    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
        output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
     output << endl;
@@ -171,5 +171,5 @@
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
     CreateForce(Fragments, Fragments.MatrixCounter);
-    for (int l=0;l<Fragments.ColumnCounter;l++)
+    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
        output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
     output << endl;
@@ -207,5 +207,5 @@
     for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) { 
       output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
-      for (int l=0;l<Fragments.ColumnCounter;l++) {
+      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
         if (((l+1) % 3) == 0) {
           norm = 0.;
@@ -252,5 +252,5 @@
     for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) { 
       output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
-      for (int l=0;l<Fragments.ColumnCounter;l++)
+      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++)
         output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
       output << endl;
@@ -278,5 +278,5 @@
       output << BondOrder+1 << "\t" << KeySet.OrderSet[BondOrder][i]+1;
       CreateFragment(Fragment, KeySet.OrderSet[BondOrder][i]);
-      for (int l=0;l<Fragment.ColumnCounter;l++)
+      for (int l=0;l<Fragment.ColumnCounter[ KeySet.OrderSet[BondOrder][i] ];l++)
         output << scientific << "\t" << Fragment.Matrix[ KeySet.OrderSet[BondOrder][i] ][ Fragment.RowCounter[ KeySet.OrderSet[BondOrder][i] ] ][l];
       output << endl;
@@ -297,5 +297,5 @@
     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
       if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
-        for (int k=Fragments.ColumnCounter;k--;)
+        for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
           Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
       }
@@ -314,5 +314,5 @@
     int i=0;
     do {  // first get a minimum value unequal to 0
-      for (int k=Fragments.ColumnCounter;k--;)
+      for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
         Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
       i++;
@@ -320,5 +320,5 @@
     for(;i<KeySet.FragmentsPerOrder[BondOrder];i++) { // then find lowest
       if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) > fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
-        for (int k=Fragments.ColumnCounter;k--;)
+        for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
           Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
       }
@@ -344,5 +344,5 @@
     CreateFragmentOrder(Fragment, KeySet, BondOrder);
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
-    for (int l=0;l<Fragment.ColumnCounter;l++)
+    for (int l=0;l<Fragment.ColumnCounter[ Fragment.MatrixCounter ];l++)
       output << scientific << "\t" << Fragment.Matrix[ Fragment.MatrixCounter ][ Fragment.RowCounter[ Fragment.MatrixCounter ]-1 ][l];
     output << endl;
@@ -358,5 +358,5 @@
 void CreateEnergy(class MatrixContainer &Energy, int MatrixNumber)
 {
-  for(int k=0;k<Energy.ColumnCounter;k++)
+  for(int k=0;k<Energy.ColumnCounter[MatrixNumber];k++)
     Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber] ] [k] =  Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber]-1 ] [k];
 };
@@ -369,7 +369,7 @@
 void CreateMinimumForce(class MatrixContainer &Force, int MatrixNumber)
 {
-  for (int l=Force.ColumnCounter;l--;)
+  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
     Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
-  for (int l=5;l<Force.ColumnCounter;l+=3) {
+  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
     double stored = 0;
     int k=0;
@@ -404,7 +404,7 @@
 {
   int divisor = 0;
-  for (int l=Force.ColumnCounter;l--;)
+  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
     Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
-  for (int l=5;l<Force.ColumnCounter;l+=3) {
+  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
     double tmp = 0;
     for (int k=Force.RowCounter[MatrixNumber];k--;) {
@@ -428,5 +428,5 @@
 void CreateMaximumForce(class MatrixContainer &Force, int MatrixNumber)
 {
-  for (int l=5;l<Force.ColumnCounter;l+=3) {
+  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
     double stored = 0;
     for (int k=Force.RowCounter[MatrixNumber];k--;) {
@@ -460,7 +460,7 @@
 void CreateVectorSumForce(class MatrixContainer &Force, int MatrixNumber)
 {
-  for (int l=Force.ColumnCounter;l--;)
+  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
     Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
-  for (int l=0;l<Force.ColumnCounter;l++) {
+  for (int l=0;l<Force.ColumnCounter[MatrixNumber];l++) {
     for (int k=Force.RowCounter[MatrixNumber];k--;)
       Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] += Force.Matrix[MatrixNumber][k][l];
@@ -542,10 +542,10 @@
 
   getline(line, token, '\t');
-  for (int i=2; i<= Energy.ColumnCounter;i++) {
+  for (int i=2; i<= Energy.ColumnCounter[Energy.MatrixCounter];i++) {
     getline(line, token, '\t');
     while (token[0] == ' ') // remove leading white spaces
       token.erase(0,1);
     output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(abs($" << i+2 << ")) " << uses;
-    if (i != (Energy.ColumnCounter))
+    if (i != (Energy.ColumnCounter[Energy.MatrixCounter]))
       output << ", \\";
     output << endl;
@@ -566,10 +566,10 @@
 
   getline(line, token, '\t');
-  for (int i=1; i<= Energy.ColumnCounter;i++) {
+  for (int i=1; i<= Energy.ColumnCounter[Energy.MatrixCounter];i++) {
     getline(line, token, '\t');
     while (token[0] == ' ') // remove leading white spaces
       token.erase(0,1);
     output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":" << i+2 << " " << uses;
-    if (i != (Energy.ColumnCounter))
+    if (i != (Energy.ColumnCounter[Energy.MatrixCounter]))
       output << ", \\";
     output << endl;
@@ -594,5 +594,5 @@
   getline(line, token, '\t');
   getline(line, token, '\t');
-  for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
+  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
     getline(line, token, '\t');
     while (token[0] == ' ') // remove leading white spaces
@@ -600,5 +600,5 @@
     token.erase(token.length(), 1);  // kill residual index char (the '0')
     output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(sqrt($" << i+1 << "*$" << i+1 << "+$" << i+2 << "*$" << i+2 << "+$" << i+3 << "*$" << i+3 << ")) " << uses;
-    if (i != (Force.ColumnCounter-1))
+    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
       output << ", \\";
     output << endl;
@@ -625,5 +625,5 @@
   getline(line, token, '\t');
   getline(line, token, '\t');
-  for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
+  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
     getline(line, token, '\t');
     while (token[0] == ' ') // remove leading white spaces
@@ -631,5 +631,5 @@
     token.erase(token.length(), 1);  // kill residual index char (the '0')
     output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(abs($" << i+1 << ")) " << uses;
-    if (i != (Force.ColumnCounter-1))
+    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
       output << ", \\";
     output << endl;
@@ -657,5 +657,5 @@
   getline(line, token, '\t');
   getline(line, token, '\t');
-  for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
+  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
     getline(line, token, '\t');
     while (token[0] == ' ') // remove leading white spaces
@@ -663,5 +663,5 @@
     token.erase(token.length(), 1);  // kill residual index char (the '0')
     output << "'" << prefix << ".dat' title '" << token << "' using ($" << xargument << "+" << fixed << setprecision(1) << (double)((i-7)/3)*0.2 << "):(sqrt($" << i+1 << "*$" << i+1 << "+$" << i+2 << "*$" << i+2 << "+$" << i+3 << "*$" << i+3 << ")) " << uses << " " << fillcolor[(i-7)/3];
-    if (i != (Force.ColumnCounter-1))
+    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
       output << ", \\";
     output << endl;
@@ -689,5 +689,5 @@
   getline(line, token, '\t');
   getline(line, token, '\t');
-  for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
+  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
     getline(line, token, '\t');
     while (token[0] == ' ') // remove leading white spaces
@@ -695,5 +695,5 @@
     token.erase(token.length(), 1);  // kill residual index char (the '0')
     output << "'" << prefix << ".dat' title '" << token << "' using ($" << xargument << "+" << fixed << setprecision(1) << (double)((i-7)/3)*0.2 << "):" << i+1 << " " << uses << " " << fillcolor[(i-7)/3];
-    if (i != (Force.ColumnCounter-1))
+    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
       output << ", \\";
     output << endl;
Index: src/joiner.cpp
===================================================================
--- src/joiner.cpp	(revision 85d278634d82c0d87206785b8cb85c4c6f794e06)
+++ src/joiner.cpp	(revision f731aed8e27c7b0f1fcfad0872e0742f3403b3b8)
@@ -80,6 +80,6 @@
 
   // ---------- Parse the TE Factors into an array -----------------
-  if (!Energy.ParseIndices()) return 1;
-  if (Hcorrected) Hcorrection.ParseIndices();
+  if (!Energy.InitialiseIndices()) return 1;
+  if (Hcorrected) Hcorrection.InitialiseIndices();
   
   // ---------- Parse the Force indices into an array ---------------
@@ -87,5 +87,5 @@
 
   // ---------- Parse the Hessian (=force) indices into an array ---------------
-  if (!Hessian.ParseIndices(argv[1])) return 1;
+  if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
 
   // ---------- Parse the shielding indices into an array ---------------
Index: src/parser.cpp
===================================================================
--- src/parser.cpp	(revision 85d278634d82c0d87206785b8cb85c4c6f794e06)
+++ src/parser.cpp	(revision f731aed8e27c7b0f1fcfad0872e0742f3403b3b8)
@@ -59,8 +59,9 @@
   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");
   Matrix[0] = NULL;
   RowCounter[0] = -1;
   MatrixCounter = 0;
-  ColumnCounter = 0;
+  ColumnCounter[0] = -1;
 };
 
@@ -70,5 +71,5 @@
   if (Matrix != NULL) {
     for(int i=MatrixCounter;i--;) {
-      if (RowCounter != NULL) {
+      if ((ColumnCounter != NULL) && (RowCounter != NULL)) {
           for(int j=RowCounter[i]+1;j--;)
             Free((void **)&Matrix[i][j], "MatrixContainer::~MatrixContainer: *Matrix[][]");
@@ -76,5 +77,5 @@
       }
     }
-    if ((RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
+    if ((ColumnCounter != NULL) && (RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
       for(int j=RowCounter[MatrixCounter]+1;j--;)
         Free((void **)&Matrix[MatrixCounter][j], "MatrixContainer::~MatrixContainer: *Matrix[MatrixCounter][]");
@@ -91,6 +92,38 @@
   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];
+    }
+  }
+  return true; 
+};
 
 /** Parsing a number of matrices.
@@ -130,13 +163,13 @@
     line >> Header;
   //cout << line.str() << endl;
-  ColumnCounter=0;
+  ColumnCounter[MatrixNr]=0;
   while ( getline(line,token, '\t') ) {
     if (token.length() > 0)
-      ColumnCounter++;
+      ColumnCounter[MatrixNr]++;
   }
   //cout << line.str() << endl;
-  //cout << "ColumnCounter: " << ColumnCounter << "." << endl;
-  if (ColumnCounter == 0)
-    cerr << "ColumnCounter: " << ColumnCounter << " from file " << name << ", this is probably an error!" << 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;
   
   // scan rest for number of rows/lines
@@ -150,10 +183,10 @@
     }
   }
-  //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;
   
   // allocate matrix if it's not zero dimension in one direction
-  if ((ColumnCounter > 0) && (RowCounter[MatrixNr] > -1)) {
+  if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
     Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
   
@@ -168,5 +201,5 @@
     strncpy(Header, line.str().c_str(), 1023);  
     for(int j=0;j<RowCounter[MatrixNr];j++) {
-      Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+      Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
       input.getline(filename, 1023);
       stringstream lines(filename);
@@ -174,15 +207,15 @@
       for(int k=skipcolumns;k--;)
         lines >> filename;
-      for(int k=0;(k<ColumnCounter) && (!lines.eof());k++) {
+      for(int k=0;(k<ColumnCounter[MatrixNr]) && (!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, "MatrixContainer::ParseFragmentMatrix: *Matrix[RowCounter[MatrixNr]][]");
-      for(int j=ColumnCounter;j--;)
+      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseFragmentMatrix: *Matrix[RowCounter[MatrixNr]][]");
+      for(int j=ColumnCounter[MatrixNr];j--;)
         Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
     }
   } else {
-    cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl; 
+    cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl; 
   }
   input.close();
@@ -235,7 +268,9 @@
   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;
     RowCounter[i] = -1;
+    ColumnCounter[i] = -1;
   }
   for(int i=0; i < MatrixCounter;i++) {
@@ -255,8 +290,8 @@
  * \param MCounter number of matrices
  * \param *RCounter number of rows for each matrix
- * \param CCounter number of columns for all matrices
+ * \param *CCounter number of columns for each 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");
@@ -264,13 +299,14 @@
 
   MatrixCounter = MCounter;
-  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");
+  ColumnCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *ColumnCounter");
   for(int i=MatrixCounter+1;i--;) {
     RowCounter[i] = RCounter[i];
+    ColumnCounter[i] = CCounter[i];
     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, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
-      for(int k=ColumnCounter;k--;)
+      Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter[i], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+      for(int k=ColumnCounter[i];k--;)
         Matrix[i][j][k] = 0.;
     }
@@ -286,5 +322,5 @@
   for(int i=MatrixCounter+1;i--;)
     for(int j=RowCounter[i]+1;j--;)
-      for(int k=ColumnCounter;k--;)
+      for(int k=ColumnCounter[i];k--;)
         Matrix[i][j][k] = 0.;
    return true;
@@ -299,5 +335,5 @@
   for(int i=MatrixCounter+1;i--;)
     for(int j=RowCounter[i]+1;j--;)
-      for(int k=ColumnCounter;k--;)
+      for(int k=ColumnCounter[i];k--;)
         if (fabs(Matrix[i][j][k]) > max)
           max = fabs(Matrix[i][j][k]);
@@ -315,5 +351,5 @@
   for(int i=MatrixCounter+1;i--;)
     for(int j=RowCounter[i]+1;j--;)
-      for(int k=ColumnCounter;k--;)
+      for(int k=ColumnCounter[i];k--;)
         if (fabs(Matrix[i][j][k]) < min)
           min = fabs(Matrix[i][j][k]);
@@ -331,5 +367,5 @@
 {
   for(int j=RowCounter[MatrixCounter]+1;j--;)
-    for(int k=skipcolumns;k<ColumnCounter;k++)
+    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
       Matrix[MatrixCounter][j][k] = value;
    return true;
@@ -344,5 +380,5 @@
 {
   for(int j=RowCounter[MatrixCounter]+1;j--;)
-    for(int k=skipcolumns;k<ColumnCounter;k++)
+    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
       Matrix[MatrixCounter][j][k] = values[j][k];
    return true;
@@ -351,5 +387,5 @@
 /** 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
@@ -384,12 +420,12 @@
               }
               if (Order == SubOrder) { // equal order is always copy from Energies
-                for(int l=ColumnCounter;l--;) // then adds/subtract each column 
+                for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];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;l--;)
+                for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;)
                   Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
               }
             }
-            //if ((ColumnCounter>1) && (RowCounter[0]-1 >= 1))
+            //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;
           }
@@ -428,5 +464,5 @@
     output << Header << endl;
     for(int j=0;j<RowCounter[i];j++) {
-      for(int k=0;k<ColumnCounter;k++)
+      for(int k=0;k<ColumnCounter[i];k++)
         output << scientific << Matrix[i][j][k] << "\t";
       output << endl;
@@ -457,5 +493,5 @@
   output << Header << endl;
   for(int j=0;j<RowCounter[MatrixCounter];j++) {
-    for(int k=0;k<ColumnCounter;k++)
+    for(int k=0;k<ColumnCounter[MatrixCounter];k++)
       output << scientific << Matrix[MatrixCounter][j][k] << "\t";
     output << endl;
@@ -467,23 +503,7 @@
 // ======================================= 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
@@ -498,10 +518,10 @@
     for(int i=KeySet.FragmentsPerOrder[Order];i--;)
       for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
-        for(int k=ColumnCounter;k--;)
+        for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];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;k--;)
+        for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
           Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
   return true;
@@ -524,12 +544,16 @@
     // count maximum of columns
     RowCounter[MatrixCounter] = 0;
-    for(int j=0; j < MatrixCounter;j++) // (energy matrix might be bigger than number of atoms in terms of rows)
+    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
-    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter << " columns." << endl;
+    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, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     
     // try independently to parse global energysuffix file if present
@@ -589,5 +613,5 @@
 
 /** 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
@@ -609,5 +633,5 @@
       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;k++)
+        for(int k=2;k<ColumnCounter[MatrixCounter];k++)
           Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
       }
@@ -655,10 +679,16 @@
     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 << " columns." << endl;
+    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, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
 
     // try independently to parse global forcesuffix file if present
@@ -686,5 +716,5 @@
   stringstream line;
   
-  cout << "Parsing force indices for " << MatrixCounter << " matrices." << endl;
+  cout << "Parsing hessian indices for " << MatrixCounter << " matrices." << endl;
   Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "HessianMatrix::ParseIndices: **Indices");
   line << name << FRAGMENTPREFIX << FORCESFILE;
@@ -720,5 +750,5 @@
 
 /** 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 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
@@ -739,11 +769,8 @@
       }
       if (j != -1) {
-        //for(int m=0;m<ColumnCounter[ FragmentNr ];m++) {
-        for(int m=0;m<ColumnCounter;m++) {
+        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] << "!" << endl;
-          if (k > ColumnCounter) {
-            cerr << "Current hessian index " << k << " is greater than " << ColumnCounter << "!" << endl;
+          if (k > ColumnCounter[MatrixCounter]) {
+            cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << "!" << endl;
             return false;
           }
@@ -784,4 +811,5 @@
     }
     RowCounter[MatrixCounter] = 0;
+    ColumnCounter[MatrixCounter] = 0;
     while (!input.eof()) {
       input.getline(filename, 1023);
@@ -790,16 +818,19 @@
         zeile >> nr;
         //cout << "Current index: " << nr << "." << endl;
-        if (nr > RowCounter[MatrixCounter])
+        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 << " columns." << endl;
+    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, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
 
     // try independently to parse global forcesuffix file if present
Index: src/parser.hpp
===================================================================
--- src/parser.hpp	(revision 85d278634d82c0d87206785b8cb85c4c6f794e06)
+++ src/parser.hpp	(revision f731aed8e27c7b0f1fcfad0872e0742f3403b3b8)
@@ -48,12 +48,13 @@
     int MatrixCounter;
     int *RowCounter;
-    int ColumnCounter;
+    int *ColumnCounter;
   
   MatrixContainer();
   ~MatrixContainer();
   
+  bool InitialiseIndices(class MatrixContainer *Matrix = NULL);
   bool ParseMatrix(const char *name, int skiplines, int skipcolumns, int MatrixNr);
   virtual bool ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns);
-  bool AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int CCounter);
+  bool AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int *CCounter);
   bool ResetMatrix();
   double FindMinValue();
@@ -72,5 +73,4 @@
 class EnergyMatrix : public MatrixContainer {
   public:
-    bool ParseIndices();
     bool SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySet, int Order, double sign);
     bool ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns);
