Index: src/analyzer.cpp
===================================================================
--- src/analyzer.cpp	(revision c685eb762cafcd3ac83229798d1b5348383442c3)
+++ src/analyzer.cpp	(revision 3902484c207d407476dea9c34dc0482f306011aa)
@@ -25,4 +25,5 @@
   periodentafel *periode = NULL; // and a period table of all elements
   EnergyMatrix Energy;
+  EnergyMatrix Hcorrection;
   ForceMatrix Force;
   ForceMatrix Shielding;
@@ -30,4 +31,5 @@
   MatrixContainer Time;
   EnergyMatrix EnergyFragments;
+  EnergyMatrix HcorrectionFragments;
   ForceMatrix ForceFragments;
   ForceMatrix ShieldingFragments;
@@ -45,5 +47,9 @@
   stringstream Orderxrange;
   stringstream Fragmentxrange;
-
+  stringstream yrange;
+  char *dir = NULL;
+  bool Hcorrected = true;
+  double norm;
+  
   cout << "ANOVA Analyzer" << endl;
   cout << "==============" << endl;
@@ -57,5 +63,10 @@
     cout << "[elementsdb]\tpath to elements database, needed for shieldings." << endl;
     return 1;
-  }
+  } else {
+    dir = (char *) Malloc(sizeof(char)*(strlen(argv[2])+2), "main: *dir");
+    strcpy(dir, "/");
+    strcat(dir, argv[2]);
+  }
+  
   if (argc > 4) {
     cout << "Loading periodentafel." << endl;
@@ -71,17 +82,21 @@
   
   // ------------- Parse through all Fragment subdirs --------
-  if (!Energy.ParseMatrix(argv[1], argv[2], EnergySuffix,0,0)) return 1;
-  if (!Force.ParseMatrix(argv[1], argv[2], ForcesSuffix,0,0)) return 1;
-  if (!Time.ParseMatrix(argv[1], argv[2], TimeSuffix, 10,1)) return 1;
+  if (!Energy.ParseMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
+  Hcorrected = Hcorrection.ParseMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0);
+  if (!Force.ParseMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
+  //if (!Time.ParseMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1;
   if (periode != NULL) { // also look for PAS values
-    if (!Shielding.ParseMatrix(argv[1], argv[2], ShieldingSuffix, 1, 0)) return 1;
-    if (!ShieldingPAS.ParseMatrix(argv[1], argv[2], ShieldingPASSuffix, 1, 0)) return 1;
+    if (!Shielding.ParseMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
+    if (!ShieldingPAS.ParseMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
   }
 
   // ---------- Parse the TE Factors into an array -----------------
   if (!Energy.ParseIndices()) return 1;
+  if (Hcorrected) Hcorrection.ParseIndices();
   
   // ---------- Parse the Force indices into an array ---------------
   if (!Force.ParseIndices(argv[1])) return 1;
+  if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
+  if (!ForceFragments.ParseIndices(argv[1])) return 1;
 
   // ---------- Parse the shielding indices into an array ---------------
@@ -96,9 +111,10 @@
   
   // ---------- Parse fragment files created by 'joiner' into an array -------------
-  if (!EnergyFragments.ParseMatrix(argv[1], argv[2], EnergyFragmentSuffix,0,0)) return 1;
-  if (!ForceFragments.ParseMatrix(argv[1], argv[2], ForceFragmentSuffix,0,0)) return 1;
+  if (!EnergyFragments.ParseMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
+  if (Hcorrected) HcorrectionFragments.ParseMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
+  if (!ForceFragments.ParseMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
   if (periode != NULL) { // also look for PAS values
-    if (!ShieldingFragments.ParseMatrix(argv[1], argv[2], ShieldingSuffix, 1, 0)) return 1;
-    if (!ShieldingPASFragments.ParseMatrix(argv[1], argv[2], ShieldingPASSuffix, 1, 0)) return 1;
+    if (!ShieldingFragments.ParseMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
+    if (!ShieldingPASFragments.ParseMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
   }
 
@@ -143,18 +159,18 @@
   }
   
-  output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
-  Time.SetLastMatrix(0., 0);
-  for (int BondOrder=KeySet.Order;BondOrder--;)
-    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
-      for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
-        for(int k=Time.ColumnCounter;k--;) {
-          Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
-        }
-  for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
-    for(int k=0;k<Time.ColumnCounter;k++)
-      output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
-    output << endl;
-  }
-  output << endl;
+//  output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
+//  Time.SetLastMatrix(0., 0);
+//  for (int BondOrder=KeySet.Order;BondOrder--;)
+//    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
+//      for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
+//        for(int k=Time.ColumnCounter;k--;) {
+//          Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
+//        }
+//  for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
+//    for(int k=0;k<Time.ColumnCounter;k++)
+//      output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
+//    output << endl;
+//  }
+//  output << endl;
   output.close();
 
@@ -165,150 +181,53 @@
   // ======================================= Creating the data files ==============================================================
 
-  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
-  if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
-  Time.SetLastMatrix(0., 0);
-  output << "#Order\tFrag.No.\t" << Time.Header << endl;
-  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
-    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
-      for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
-        for(int k=Time.ColumnCounter;k--;) {
-          Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
-        }
-    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
-    for(int k=0;k<Time.ColumnCounter;k++)
-      output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][1][k];
-    output << endl;
-  }
-  output.close();
+//  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
+//  if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
+//  Time.SetLastMatrix(0., 0);
+//  output << "#Order\tFrag.No.\t" << Time.Header << endl;
+//  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+//    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
+//      for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
+//        for(int k=Time.ColumnCounter;k--;) {
+//          Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
+//        }
+//    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
+//    for(int k=0;k<Time.ColumnCounter;k++)
+//      output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][1][k];
+//    output << endl;
+//  }
+//  output.close();
 
   // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
-  if (!EnergyFragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0)) return 1;
-  if (!CreateDataEnergyOrder(Energy, EnergyFragments, KeySet, argv[3], "DeltaEnergies-Order", "Plot of error between approximated and full energies energies versus the Bond Order", datum)) return 1;
+  if (!CreateDataDeltaEnergyOrder(Energy, EnergyFragments, KeySet, argv[3], "DeltaEnergies-Order", "Plot of error between approximated and full energies energies versus the Bond Order", datum)) return 1;
   
   // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
-  if (!EnergyFragments.SetLastMatrix(1.,0)) return 1;
-  if (!CreateDataEnergyOrder(Energy, EnergyFragments, KeySet, argv[3], "Energies-Order", "Plot of approximated energies versus the Bond Order", datum)) return 1;
+  if (!CreateDataEnergyOrder(EnergyFragments, KeySet, argv[3], "Energies-Order", "Plot of approximated energies versus the Bond Order", datum)) return 1;
   
   // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
-  if (!ForceFragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0)) return 1;
-  if (!OpenOutputFile(output, argv[3], "DeltaForces-Order.dat" )) return false;
-  if (!OpenOutputFile(output2, argv[3], "DeltaMinForces-Order.dat" )) return false;
-  if (!OpenOutputFile(output3, argv[3], "DeltaMeanForces-Order.dat" )) return false;
-  if (!OpenOutputFile(output4, argv[3], "DeltaMaxForces-Order.dat" )) return false;
-  cout << "Plot of per atom and min/mean/max error between approximated forces and full forces versus the Bond Order ... " << endl;
-  output << "# Plot of error between approximated forces and full forces versus the Bond Order, created on " << datum;
-  output << "# AtomNo" << Force.Header << endl;
-  output2 << "# Plot of min error between approximated forces and full forces versus the Bond Order, created on " << datum;
-  output2 << "# Order\tFrag.No.\t" << Force.Header << endl;
-  output3 << "# Plot of mean error between approximated forces and full forces versus the Bond Order, created on " << datum;
-  output3 << "# Order\tFrag.No.\t" << Force.Header << endl;
-  output4 << "# Plot of max error between approximated forces and full forces versus the Bond Order, created on " << datum;
-  output4 << "# Order\tFrag.No.\t" << Force.Header << endl;
-  int FragmentCounter = 0;
-  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
-    for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
-      for(int l=0;l<Force.RowCounter[ KeySet.OrderSet[BondOrder][i] ];l++) {
-        int j = Force.Indices[ KeySet.OrderSet[BondOrder][i] ][l];
-        if (j > Force.RowCounter[Force.MatrixCounter]) {
-          cerr << "Current force index " << j << " is greater than " << Force.RowCounter[Force.MatrixCounter] << "!" << endl;
-          return 1;
-        }
-        if (j != -1)
-          for(int k=0;k<Force.ColumnCounter;k++) {
-            Force.Matrix[Force.MatrixCounter][j][k] -= ForceFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][l][k];
-          }
-      }
-      FragmentCounter++;
-    }
-    // errors per atom
-    output << "#Order\t" << BondOrder+1 << endl;
-    for(int j=0;j<Force.RowCounter[ Force.MatrixCounter ];j++) { 
-      output << Force.Indices[Force.MatrixCounter][j] << "\t";
-      for (int l=0;l<Force.ColumnCounter;l++)
-        if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ j ][l]) < MYEPSILON)
-          output << scientific << Force.Matrix[Force.MatrixCounter][ j ][l] << "\t";
-        else
-          output << scientific << (Force.Matrix[Force.MatrixCounter][ j ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ j ][l]) << "\t";
-      output << endl;
-    }
-    output << endl;
-    
-    // min error
-    output2 << BondOrder+1 << "\t" << FragmentCounter;
-    CreateMinimumForce(Force, Force.MatrixCounter);
-    for (int l=0;l<Force.ColumnCounter;l++)
-      if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]) < MYEPSILON)
-        output2 << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
-      else
-        output2 << scientific << "\t" << (Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]);
-    output2 << endl;
-    
-    // mean error
-    output3 << BondOrder+1 << "\t" << FragmentCounter;
-    CreateMeanForce(Force, Force.MatrixCounter);
-    for (int l=0;l<Force.ColumnCounter;l++)
-      if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]) < MYEPSILON)
-        output3 << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
-      else
-        output3 << scientific << "\t" << (Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]);
-    output3 << endl;
-    
-    // max error
-    output4 << BondOrder+1 << "\t" << FragmentCounter;
-    CreateMaximumForce(Force, Force.MatrixCounter);
-    for (int l=0;l<Force.ColumnCounter;l++)
-      if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]) < MYEPSILON)
-        output4 << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
-      else
-        output4 << scientific << "\t" << (Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]);
-    output4 << endl;
-  }
-  output.close();
-  output2.close();
-  output3.close();
-  output4.close();
+  if (!CreateDataDeltaForcesOrderPerAtom(Force, ForceFragments, KeySet, argv[3], "DeltaForces-Order", "Plot of error between approximated forces and full forces versus the Bond Order", datum)) return 1;
+
+  // min force
+  if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMinForces-Order", "Plot of min error between approximated forces and full forces versus the Bond Order", datum, CreateMinimumForce)) return 1;
+
+  // mean force
+  if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMeanForces-Order", "Plot of mean error between approximated forces and full forces versus the Bond Order", datum, CreateMeanForce)) return 1;
+
+  // max force
+  if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMaxForces-Order", "Plot of max error between approximated forces and full forces versus the Bond Order", datum, CreateMaximumForce)) return 1;
 
   // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
-  if(!OpenOutputFile(output, argv[3], "Forces-Order.dat")) return 1;
-  cout << "Plot of approximated forces versus the Bond Order ... " << endl;
-  output << "# Plot of approximated forces versus the Bond Order, created on " << datum;
-  output << "# AtomNo" << Force.Header << endl;
-  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
-    for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
-      for(int l=0;l<Force.RowCounter[ KeySet.OrderSet[BondOrder][i] ];l++) {
-        int j = Force.Indices[ KeySet.OrderSet[BondOrder][i] ][l];
-        if (j > Force.RowCounter[Force.MatrixCounter]) {
-          cerr << "Current force index " << j << " is greater than " << Force.RowCounter[Force.MatrixCounter] << "!" << endl;
-          return 1;
-        }
-        if (j != -1)
-          for(int k=0;k<Force.ColumnCounter;k++) {
-            Force.Matrix[Force.MatrixCounter][j][k] -= ForceFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][l][k];
-          }
-      }
-    }
-    // errors per atom
-    output << "#Order\t" << BondOrder+1 << endl;
-    for(int j=0;j<Force.RowCounter[ Force.MatrixCounter ];j++) { 
-      output << Force.Indices[Force.MatrixCounter][j] << "\t";
-      for (int l=0;l<Force.ColumnCounter;l++)
-         output << Force.Matrix[Force.MatrixCounter][ j ][l] << "\t";
-      output << endl;
-    }
-    output << endl;
-  }
-  output.close();
-
-  // min force
-  if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "MinForces-Order", "Plot of min approximated forces versus the Bond Order", datum, CreateMinimumForce)) return 1;
-
-  // mean force
-  if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "MeanForces-Order", "Plot of mean approximated forces versus the Bond Order", datum, CreateMeanForce)) return 1;
-
-  // max force
-  if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "MaxForces-Order", "Plot of max approximated forces versus the Bond Order", datum, CreateMaximumForce)) return 1;
+  if (!CreateDataForcesOrderPerAtom(ForceFragments, KeySet, argv[3], "Forces-Order", "Plot of approximated forces versus the Bond Order", datum)) return 1;
+
+  // min force
+  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MinForces-Order", "Plot of min approximated forces versus the Bond Order", datum, CreateMinimumForce)) return 1;
+
+  // mean force
+  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MeanForces-Order", "Plot of mean approximated forces versus the Bond Order", datum, CreateMeanForce)) return 1;
+
+  // max force
+  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MaxForces-Order", "Plot of max approximated forces versus the Bond Order", datum, CreateMaximumForce)) return 1;
 
   // ++++++++++++++++++++++++++++++++++++++Plotting vector sum (should be 0) vs. bond order
-  if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "VectorSum-Order", "Plot of vector sum of the approximated forces versus the Bond Order", datum, CreateVectorSumForce)) return 1;
+  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "VectorSum-Order", "Plot of vector sum of the approximated forces versus the Bond Order", datum, CreateVectorSumForce)) return 1;
 
   // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
@@ -342,10 +261,11 @@
   Orderxrange << "[1:" << KeySet.Order << "]";
   Fragmentxrange << "[0:" << KeySet.FragmentCounter+1 << "]";
+  yrange.str("[1e-8:1e+1]");
   
   // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
-  if (!CreatePlotOrder(Time, KeySet, argv[3], "SimTime-Order", 1, "below", "y", "",  1, 1, "bond order k", "Evaluation time [s]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1;
+  //if (!CreatePlotOrder(Time, KeySet, argv[3], "SimTime-Order", 1, "below", "y", "",  1, 1, "bond order k", "Evaluation time [s]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1;
   
   // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
-  if (!CreatePlotOrder(Energy, KeySet, argv[3], "DeltaEnergies-Order", 1, "outside", "y", "",  1, 1, "bond order k", "absolute error in energy [Ht]", Orderxrange.str().c_str(), "[1e-8:1e+1]", "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
+  if (!CreatePlotOrder(Energy, KeySet, argv[3], "DeltaEnergies-Order", 1, "outside", "y", "",  1, 1, "bond order k", "absolute error in energy [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
   
   // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
@@ -353,12 +273,13 @@
 
   // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
-  // min force
-  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMinForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in min force [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-8:1e+0]", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
-
-  // mean force
-  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMeanForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in mean force [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-8:1e+0]", "1" , "with linespoints", AbsFirstForceValuePlotLine)) return 1;
-
-  // max force
-  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMaxForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in max force [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-8:1e+0]", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
+  yrange.str("[1e-8:1e+0]");
+  // min force
+  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMinForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in min force [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
+
+  // mean force
+  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMeanForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in mean force [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsFirstForceValuePlotLine)) return 1;
+
+  // max force
+  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMaxForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in max force [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
 
   // min/mean/max comparison for total force
@@ -395,28 +316,32 @@
 
   // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
-  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-Fragment", 5, "below", "y", "", 1, 5, "fragment number", "Energies of each fragment [Ht]", Fragmentxrange.str().c_str(), "[1e-16:1e+1]", "2" , "with points", AbsEnergyPlotLine)) return 1;
-  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Energies of each fragment [Ht]", Orderxrange.str().c_str(), "[1e-10:1e+1]", "1" , "with points", AbsEnergyPlotLine)) return 1;
-  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Maximum fragment energy [Ht]", Orderxrange.str().c_str(), "[1e-6:1e+1]", "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
-  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Minimum fragment energy [Ht]", Orderxrange.str().c_str(), "[1e-6:1e+1]", "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
+  yrange.str("");
+  yrange << "[" << EnergyFragments.FindMinValue() << ":" << EnergyFragments.FindMaxValue() << "]";
+  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-Fragment", 5, "below", "y", "", 1, 5, "fragment number", "Energies of each fragment [Ht]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with points", AbsEnergyPlotLine)) return 1;
+  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Energies of each fragment [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with points", AbsEnergyPlotLine)) return 1;
+  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Maximum fragment energy [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
+  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Minimum fragment energy [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
 
   // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment
+  yrange.str("");
+  yrange << "[" << ForceFragments.FindMinValue() << ":" << ForceFragments.FindMaxValue()<< "]";
   // min
-  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "minimum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), "[1e-16:1e+1]", "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
+  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "minimum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
     
   // mean
-  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "mean of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), "[1e-16:1e+1]", "2" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
+  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "mean of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
     
   // max
-  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "maximum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), "[1e-16:1e+1]", "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
+  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "maximum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
 
   // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment per bond order
   // min
-  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "minimum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-16:1e+1]", "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
+  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "minimum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
     
   // mean
-  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "mean of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-16:1e+1]", "1" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
+  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "mean of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
     
   // max
-  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "maximum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-16:1e+1]", "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
+  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "maximum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
  
   // create Makefile
@@ -434,4 +359,5 @@
   // ++++++++++++++++ exit ++++++++++++++++++++++++++++++++++
   delete(periode);
+  Free((void **)&dir, "main: *dir");
   cout << "done." << endl;
   return 0;
Index: src/datacreator.cpp
===================================================================
--- src/datacreator.cpp	(revision c685eb762cafcd3ac83229798d1b5348383442c3)
+++ src/datacreator.cpp	(revision 3902484c207d407476dea9c34dc0482f306011aa)
@@ -30,6 +30,5 @@
 
 /** Plots an energy vs. order.
- * \param Energy EnergyMatrix class containing matrix values (last matrix is created by summing over Fragments)
- * \param EnergyFragments EnergyMatrix class containing matrix values
+ * \param &Fragments EnergyMatrix class containing matrix values
  * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
  * \param *prefix prefix in filename (without ending)
@@ -37,5 +36,5 @@
  * \return true if file was written successfully
  */
-bool CreateDataEnergyOrder(class MatrixContainer &Energy, class MatrixContainer &EnergyFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 
+bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 
 {
   stringstream filename;
@@ -46,17 +45,51 @@
   cout << msg << endl;
   output << "# " << msg << ", created on " << datum;
-  output << "#Order\tFrag.No.\t" << Energy.Header << endl;
+  output << "#Order\tFrag.No.\t" << Fragments.Header << endl;
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
-      for(int j=Energy.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
-        for(int k=Energy.ColumnCounter;k--;)
-          Energy.Matrix[Energy.MatrixCounter][j][k] -= EnergyFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
+      for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
+        for(int k=Fragments.ColumnCounter;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<Energy.ColumnCounter;l++)
-      if (fabs(EnergyFragments.Matrix[EnergyFragments.MatrixCounter][ EnergyFragments.RowCounter[EnergyFragments.MatrixCounter]-1 ][l]) < MYEPSILON)
-        output << scientific << "\t" << Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l];
+    for (int l=0;l<Fragments.ColumnCounter;l++)
+      output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+/** Plots an energy error vs. order.
+ * \param &Energy EnergyMatrix class containing reference values (in MatrixCounter matrix)
+ * \param &Fragments EnergyMatrix class containing matrix values
+ * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param *prefix prefix in filename (without ending)
+ * \param *msg message to be place in first line as a comment
+ * \return true if file was written successfully
+ */
+bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
+  cout << msg << endl;
+  output << "# " << msg << ", created on " << datum;
+  output << "#Order\tFrag.No.\t" << Fragments.Header << endl;
+  Fragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0);
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
+      for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
+        for(int k=Fragments.ColumnCounter;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++)
+      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];
       else
-        output << scientific << "\t" << (Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l] / EnergyFragments.Matrix[EnergyFragments.MatrixCounter][ EnergyFragments.RowCounter[EnergyFragments.MatrixCounter]-1 ][l]);
+        output << scientific << "\t" << (Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l] / Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]);
     output << endl;
   }
@@ -66,6 +99,12 @@
 
 /** Plot forces vs. order.
- */
-bool CreateDataForcesOrder(class MatrixContainer &Force, class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
+ * \param &Fragments ForceMatrix class containing matrix values
+ * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param *prefix prefix in filename (without ending)
+ * \param *msg message to be place in first line as a comment
+ * \param *datum current date and time
+ * \return true if file was written successfully
+ */
+bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
 {
   stringstream filename;
@@ -76,26 +115,127 @@
   cout << msg << endl;
   output << "# " << msg << ", created on " << datum;
-  output << "# Order\tFrag.No.\t" << Force.Header << endl;
-  for(int j=0;j<=Force.RowCounter[ Force.MatrixCounter ];j++)
-    for(int k=0;k<Force.ColumnCounter;k++)
-       Force.Matrix[Force.MatrixCounter][j][k] = 0.;
-  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
-    for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
-      for(int l=0;l<Force.RowCounter[ KeySet.OrderSet[BondOrder][i] ];l++) {
-        int j = Force.Indices[ KeySet.OrderSet[BondOrder][i] ][l];
-        if (j > Force.RowCounter[Force.MatrixCounter]) {
-          cerr << "Current force index " << j << " is greater than " << Force.RowCounter[Force.MatrixCounter] << "!" << endl;
-          return 1;
-        }
-        if (j != -1)
-          for(int k=Force.ColumnCounter;k--;) {
-            Force.Matrix[Force.MatrixCounter][j][k] += ForceFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][l][k];
-          }
+  output << "# Order\tFrag.No.\t" << Fragments.Header << endl;
+  Fragments.SetLastMatrix(0.,0);
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
+    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
+    CreateForce(Fragments, Fragments.MatrixCounter);
+    for (int l=0;l<Fragments.ColumnCounter;l++)
+       output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+/** Plot forces error vs. order.
+ * \param &Force ForceMatrix containing reference values (in MatrixCounter matrix)
+ * \param &Fragments ForceMatrix class containing matrix values
+ * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param *prefix prefix in filename (without ending)
+ * \param *msg message to be place in first line as a comment
+ * \param *datum current date and time
+ * \return true if file was written successfully
+ */
+bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
+  cout << msg << endl;
+  output << "# " << msg << ", created on " << datum;
+  output << "# Order\tFrag.No.\t" << Fragments.Header << endl;
+  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0);
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
+    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
+    CreateForce(Fragments, Fragments.MatrixCounter);
+    for (int l=0;l<Fragments.ColumnCounter;l++)
+       output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+/** Plot forces error vs. vs atom vs. order.
+ * \param &Force ForceMatrix containing reference values (in MatrixCounter matrix)
+ * \param &Fragments ForceMatrix class containing matrix values
+ * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param *prefix prefix in filename (without ending)
+ * \param *msg message to be place in first line as a comment
+ * \param *datum current date and time
+ * \return true if file was written successfully
+ */
+bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 
+{
+  stringstream filename;
+  ofstream output;
+  double norm = 0.;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
+  cout << msg << endl;
+  output << "# " << msg << ", created on " << datum;
+  output << "# AtomNo\t" << Fragments.Header << endl;
+  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter], 0);
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    //cout << "Current order is " << BondOrder << "." << endl;
+    Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
+    // errors per atom
+    output << "#Order\t" << BondOrder+1 << endl;
+    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++) {
+        if (((l+1) % 3) == 0) {
+          norm = 0.;
+          for (int m=0;m<NDIM;m++)
+            norm += Force.Matrix[Force.MatrixCounter][ j ][l+m]*Force.Matrix[Force.MatrixCounter][ j ][l+m];
+          norm = sqrt(norm);
+        }                                                                                                           
+//        if (norm < MYEPSILON)
+          output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
+//        else
+//          output << scientific << (Fragments.Matrix[Fragments.MatrixCounter][ j ][l] / norm) << "\t";
       }
-    }
-    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
-    CreateForce(Force, Force.MatrixCounter);
-    for (int l=0;l<Force.ColumnCounter;l++)
-       output << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
+      output << endl;
+    }
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+/** Plot forces error vs. vs atom vs. order.
+ * \param &Fragments ForceMatrix class containing matrix values
+ * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param *prefix prefix in filename (without ending)
+ * \param *msg message to be place in first line as a comment
+ * \param *datum current date and time
+ * \return true if file was written successfully
+ */
+bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
+  cout << msg << endl;
+  output << "# " << msg << ", created on " << datum;
+  output << "# AtomNo\t" << Fragments.Header << endl;
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    //cout << "Current order is " << BondOrder << "." << endl;
+    Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
+    // errors per atom
+    output << "#Order\t" << BondOrder+1 << endl;
+    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++)
+        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
+      output << endl;
+    }
     output << endl;
   }
@@ -286,4 +426,13 @@
 };
 
+/** Leaves the Force.Matrix as it is.
+ * \param Force ForceMatrix class containing matrix values
+ * \param MatrixNumber the index for the ForceMatrix::matrix array
+*/
+void CreateSameForce(class MatrixContainer &Force, int MatrixNumber)
+{
+  // does nothing
+};
+
 /** Adds vectorwise all forces.
  * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
@@ -295,5 +444,5 @@
   for (int l=Force.ColumnCounter;l--;)
     Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
-  for (int l=5;l<Force.ColumnCounter;l++) {
+  for (int l=0;l<Force.ColumnCounter;l++) {
     for (int k=Force.RowCounter[MatrixNumber];k--;)
       Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] += Force.Matrix[MatrixNumber][k][l];
@@ -375,8 +524,8 @@
 
   line >> item;
-  for (int i=3; i< Energy.ColumnCounter;i++) {
-    line >> item;
-    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":(abs($" << i+1 << ")) " << uses;
-    if (i != (Energy.ColumnCounter-1))
+  for (int i=2; i<= Energy.ColumnCounter;i++) {
+    line >> item;
+    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":(abs($" << i+2 << ")) " << uses;
+    if (i != (Energy.ColumnCounter))
       output << ", \\";
     output << endl;
@@ -397,8 +546,8 @@
 
   line >> item;
-  for (int i=3; i< Energy.ColumnCounter;i++) {
-    line >> item;
-    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":" << i+1 << " " << uses;
-    if (i != (Energy.ColumnCounter-1))
+  for (int i=2; i<= Energy.ColumnCounter;i++) {
+    line >> item;
+    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":" << i+2 << " " << uses;
+    if (i != (Energy.ColumnCounter))
       output << ", \\";
     output << endl;
Index: src/datacreator.hpp
===================================================================
--- src/datacreator.hpp	(revision c685eb762cafcd3ac83229798d1b5348383442c3)
+++ src/datacreator.hpp	(revision 3902484c207d407476dea9c34dc0482f306011aa)
@@ -19,6 +19,10 @@
 bool OpenOutputFile(ofstream &output, const char *dir, const char *filename);
 
-bool CreateDataEnergyOrder(class MatrixContainer &Energy, class MatrixContainer &EnergyFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum); 
-bool CreateDataForcesOrder(class MatrixContainer &Force, class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int));
+bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum); 
+bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum); 
+bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int));
+bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum);
+bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int));
+bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum); 
 bool CreateDataFragment(class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int));
 bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int));
@@ -30,4 +34,5 @@
 void CreateMeanForce(class MatrixContainer &Force, int MatrixNumber);
 void CreateMaximumForce(class MatrixContainer &Force, int MatrixNumber);
+void CreateSameForce(class MatrixContainer &Force, int MatrixNumber);
 void CreateVectorSumForce(class MatrixContainer &Force, int MatrixNumber);
 
Index: src/defs.hpp
===================================================================
--- src/defs.hpp	(revision c685eb762cafcd3ac83229798d1b5348383442c3)
+++ src/defs.hpp	(revision 3902484c207d407476dea9c34dc0482f306011aa)
@@ -49,4 +49,6 @@
 #define TEFACTORSFILE "TE-Factors.dat"    //!< default filename of BOSSANOVA total energy factors file
 #define FORCESFILE "Forces-Factors.dat"    //!< default filename of BOSSANOVA force factors file
+#define HCORRECTIONSUFFIX "Hcorrection.dat"    //!< default filename of BOSSANOVA H correction file (unwanted saturation interaction)
+#define FITCONSTANTSUFFIX "FitConstant.dat"   //!< suffix of default filename of BOSSANOVA fit constants file (unwanted saturation interaction)
 #define SHIELDINGSUFFIX "sigma_all.csv"                //!< default filename of BOSSANOVA shieldings file
 #define SHIELDINGPASSUFFIX "sigma_all_PAS.csv"                 //!< default filename of BOSSANOVA shieldings PAS file
Index: src/joiner.cpp
===================================================================
--- src/joiner.cpp	(revision c685eb762cafcd3ac83229798d1b5348383442c3)
+++ src/joiner.cpp	(revision 3902484c207d407476dea9c34dc0482f306011aa)
@@ -19,6 +19,8 @@
   periodentafel *periode = NULL; // and a period table of all elements
   EnergyMatrix Energy;
+  EnergyMatrix Hcorrection;
   ForceMatrix Force;
   EnergyMatrix EnergyFragments;
+  EnergyMatrix HcorrectionFragments;
   ForceMatrix ForceFragments;
   ForceMatrix Shielding;
@@ -28,4 +30,6 @@
   KeySetsContainer KeySet;  
   stringstream prefix;
+  char *dir = NULL;
+  bool Hcorrected = true;
 
   cout << "Joiner" << endl;
@@ -39,4 +43,8 @@
     cout << "[elementsdb]\tpath to elements database, needed for shieldings." << endl;
     return 1;
+  } else {
+    dir = (char *) Malloc(sizeof(char)*(strlen(argv[2])+2), "main: *dir");
+    strcpy(dir, "/");
+    strcat(dir, argv[2]);
   }
   if (argc > 3) {
@@ -52,13 +60,15 @@
   
   // ------------- Parse through all Fragment subdirs --------
-  if (!Energy.ParseMatrix(argv[1], argv[2], EnergySuffix, 0,0)) return 1;
-  if (!Force.ParseMatrix(argv[1], argv[2], ForcesSuffix, 0,0)) return 1;
+  if (!Energy.ParseMatrix(argv[1], dir, EnergySuffix, 0,0)) return 1;
+  Hcorrected = Hcorrection.ParseMatrix(argv[1], "", HCORRECTIONSUFFIX, 0,0);
+  if (!Force.ParseMatrix(argv[1], dir, ForcesSuffix, 0,0)) return 1;
   if (periode != NULL) { // also look for PAS values
-    if (!Shielding.ParseMatrix(argv[1], argv[2], ShieldingSuffix, 1, 0)) return 1;
-    if (!ShieldingPAS.ParseMatrix(argv[1], argv[2], ShieldingPASSuffix, 1, 0)) return 1;
+    if (!Shielding.ParseMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
+    if (!ShieldingPAS.ParseMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
   }
 
   // ---------- Parse the TE Factors into an array -----------------
   if (!Energy.ParseIndices()) return 1;
+  if (Hcorrected) Hcorrection.ParseIndices();
   
   // ---------- Parse the Force indices into an array ---------------
@@ -76,4 +86,5 @@
   if (!KeySet.ParseManyBodyTerms()) return 1;
   if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return 1;
+  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 (periode != NULL) { // also look for PAS values
@@ -86,5 +97,5 @@
   if(!Force.SetLastMatrix(0., 2)) return 1;
   if (periode != NULL) { // also look for PAS values
-    if(!Shielding.SetLastMatrix(0., 0)) return 1;
+    if(!Shielding.SetLastMatrix(0., 2)) return 1;
     if(!ShieldingPAS.SetLastMatrix(0., 2)) return 1;
   }
@@ -97,20 +108,25 @@
     cout << "Summing energy of order " << BondOrder+1 << " ..." << endl;
     if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return 1;
-    if (!Energy.SumSubEnergy(EnergyFragments, KeySet, BondOrder)) return 1;
+    if (Hcorrected) { 
+      HcorrectionFragments.SumSubManyBodyTerms(Hcorrection, KeySet, BondOrder);
+      if (!Energy.SumSubEnergy(EnergyFragments, &HcorrectionFragments, KeySet, BondOrder, 1.)) return 1;
+      if (Hcorrected) Hcorrection.SumSubEnergy(HcorrectionFragments, NULL, KeySet, BondOrder, 1.);
+    } else 
+      if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return 1;
     // --------- sum up Forces --------------------
     cout << "Summing forces of order " << BondOrder+1 << " ..." << endl;
     if (!ForceFragments.SumSubManyBodyTerms(Force, KeySet, BondOrder)) return 1;
-    if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder)) return 1;
+    if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return 1;
     if (periode != NULL) { // also look for PAS values
       cout << "Summing shieldings of order " << BondOrder+1 << " ..." << endl;
       if (!ShieldingFragments.SumSubManyBodyTerms(Shielding, KeySet, BondOrder)) return 1;
-      if (!Shielding.SumSubForces(ShieldingFragments, KeySet, BondOrder)) return 1;
+      if (!Shielding.SumSubForces(ShieldingFragments, KeySet, BondOrder, 1.)) return 1;
       if (!ShieldingPASFragments.SumSubManyBodyTerms(ShieldingPAS, KeySet, BondOrder)) return 1;
-      if (!ShieldingPAS.SumSubForces(ShieldingPASFragments, KeySet, BondOrder)) return 1;
+      if (!ShieldingPAS.SumSubForces(ShieldingPASFragments, KeySet, BondOrder, 1.)) return 1;
     }
 
     // --------- write the energy and forces file --------------------
     prefix.str(" ");
-    prefix << argv[2] << OrderSuffix << (BondOrder+1);
+    prefix << dir << OrderSuffix << (BondOrder+1);
     cout << "Writing files " << argv[1] << prefix.str() << ". ..." << endl;
     // energy
@@ -126,29 +142,36 @@
   // fragments
   prefix.str(" ");
-  prefix << argv[2] << EnergyFragmentSuffix;
+  prefix << dir << EnergyFragmentSuffix;
   if (!EnergyFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+  if (Hcorrected) {
+    prefix.str(" ");
+    prefix << dir << HcorrectionFragmentSuffix;
+    if (!HcorrectionFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+  }
   prefix.str(" ");
-  prefix << argv[2] << ForceFragmentSuffix;
+  prefix << dir << ForceFragmentSuffix;
   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 (periode != NULL) { // also look for PAS values
     prefix.str(" ");
-    prefix << argv[2] << ShieldingFragmentSuffix;
+    prefix << dir << ShieldingFragmentSuffix;
     if (!ShieldingFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
     prefix.str(" ");
-    prefix << argv[2] << ShieldingPASFragmentSuffix;
+    prefix << dir << ShieldingPASFragmentSuffix;
     if (!ShieldingPASFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
   }
 
   // write last matrices as fragments into central dir (not subdir as above), for analyzer to know index bounds
-  if (!Energy.WriteLastMatrix(argv[1], argv[2], EnergyFragmentSuffix)) return 1;
-  if (!Force.WriteLastMatrix(argv[1], argv[2], ForceFragmentSuffix)) return 1;
+  if (!Energy.WriteLastMatrix(argv[1], dir, EnergyFragmentSuffix)) return 1;
+  if (Hcorrected) Hcorrection.WriteLastMatrix(argv[1], dir, HcorrectionFragmentSuffix);
+  if (!Force.WriteLastMatrix(argv[1], dir, ForceFragmentSuffix)) return 1;
   if (periode != NULL) { // also look for PAS values
-    if (!Shielding.WriteLastMatrix(argv[1], argv[2], ShieldingFragmentSuffix)) return 1;
-    if (!ShieldingPAS.WriteLastMatrix(argv[1], argv[2], ShieldingPASFragmentSuffix)) return 1;
+    if (!Shielding.WriteLastMatrix(argv[1], dir, ShieldingFragmentSuffix)) return 1;
+    if (!ShieldingPAS.WriteLastMatrix(argv[1], dir, ShieldingPASFragmentSuffix)) return 1;
   }
 
   // exit  
   delete(periode);
+  Free((void **)&dir, "main: *dir");
   cout << "done." << endl;
   return 0;
Index: src/moleculelist.cpp
===================================================================
--- src/moleculelist.cpp	(revision c685eb762cafcd3ac83229798d1b5348383442c3)
+++ src/moleculelist.cpp	(revision 3902484c207d407476dea9c34dc0482f306011aa)
@@ -133,4 +133,191 @@
 };
 
+/** Calculates necessary hydrogen correction due to unwanted interaction between saturated ones.
+ * If for a pair of two hydrogen atoms a and b, at least is a saturated one, and a and b are not
+ * bonded to the same atom, then we add for this pair a correction term constructed from a Morse
+ * potential function fit to QM calculations with respecting to the interatomic hydrogen distance.
+ * \param *out output stream for debugging
+ * \param *path path to file
+ */
+bool MoleculeListClass::AddHydrogenCorrection(ofstream *out, char *path)
+{
+  atom *Walker = NULL;
+  atom *Runner = NULL;
+  double ***FitConstant = NULL, **correction = NULL;
+  int a,b;
+  ofstream output;
+  ifstream input;
+  string line;
+  stringstream zeile;
+  double distance;
+  char ParsedLine[1023];
+  double tmp;
+  char *FragmentNumber = NULL;
+
+  cout << Verbose(1) << "Saving hydrogen saturation correction ... ";
+  // 0. parse in fit constant files that should have the same dimension as the final energy files
+  // 0a. find dimension of matrices with constants
+  line = path;
+  line.append("/");
+  line += FRAGMENTPREFIX;
+  line += "1";
+  line += FITCONSTANTSUFFIX;
+  input.open(line.c_str());
+  if (input == NULL) {
+    cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
+    return false;
+  }
+  a=0;
+  b=-1; // we overcount by one 
+  while (!input.eof()) {
+    input.getline(ParsedLine, 1023);
+    zeile.str(ParsedLine);
+    int i=0;
+    while (!zeile.eof()) {
+      zeile >> distance;
+      i++;
+    }
+    if (i > a) 
+      a = i; 
+    b++;
+  }
+  cout << "I recognized " << a << " columns and " << b << " rows, ";
+  input.close();
+  
+  // 0b. allocate memory for constants
+  FitConstant = (double ***) Malloc(sizeof(double **)*3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
+  for (int k=0;k<3;k++) {
+    FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
+    for (int i=a;i--;) {
+      FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
+    }
+  }
+  // 0c. parse in constants
+  for (int i=0;i<3;i++) {
+    line = path;
+    line.append("/");
+    line += FRAGMENTPREFIX;
+    sprintf(ParsedLine, "%d", i+1);
+    line += ParsedLine;
+    line += FITCONSTANTSUFFIX;
+    input.open(line.c_str());
+    if (input == NULL) {
+      cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
+      return false;
+    }
+    int k = 0,l;
+    while ((!input.eof()) && (k < b)) {
+      input.getline(ParsedLine, 1023);
+      //cout << "Current Line: " << ParsedLine << endl;
+      zeile.str(ParsedLine);
+      zeile.clear();
+      l = 0;
+      while ((!zeile.eof()) && (l < a)) {
+        zeile >> FitConstant[i][l][k];
+        //cout << FitConstant[i][l][k] << "\t";
+        l++;
+      }
+      //cout << endl;
+      k++;
+    }
+    input.close();
+  }
+  for(int k=0;k<3;k++) {
+    cout << "Constants " << k << ":" << endl;
+    for (int j=0;j<b;j++) {
+      for (int i=0;i<a;i++) {
+        cout << FitConstant[k][i][j] << "\t";
+      }
+      cout << endl;
+    }
+    cout << endl;
+  }
+  
+  // 0d. allocate final correction matrix
+  correction = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **correction");
+  for (int i=a;i--;)
+    correction[i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
+        
+  // 1a. go through every molecule in the list
+  for(int i=NumberOfMolecules;i--;) {
+    // 1b. zero final correction matrix
+    for (int k=a;k--;)
+      for (int j=b;j--;)
+        correction[k][j] = 0.;
+    // 2. take every hydrogen that is a saturated one
+    Walker = ListOfMolecules[i]->start;
+    while (Walker->next != ListOfMolecules[i]->end) {
+      Walker = Walker->next;
+      //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
+      if ((Walker->type->Z == 1) && ((Walker->father == NULL) || (Walker->father->type->Z != 1))) { // if it's a hydrogen
+        Runner = ListOfMolecules[i]->start;
+        while (Runner->next != ListOfMolecules[i]->end) {
+          Runner = Runner->next;
+          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
+          // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
+          if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && (ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) {  // (hydrogens have only one bonding partner!)
+            // 4. evaluate the morse potential for each matrix component and add up
+            distance = sqrt(Runner->x.Distance(&Walker->x));
+            //cout << "Fragment " << i << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
+            for(int k=0;k<a;k++) {
+              for (int j=0;j<b;j++) {
+                switch(k) {
+                case 1:
+                case 7:
+                case 11:
+                  tmp = pow(FitConstant[0][k][j] * ( 1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]) ) ),2);
+                  break;
+                default:
+                  tmp = FitConstant[0][k][j] * pow( distance,  FitConstant[1][k][j]) + FitConstant[2][k][j];
+                };
+                correction[k][j] -= tmp;    // ground state is actually lower (disturbed by additional interaction)
+                //cout << tmp << "\t";
+              }
+              //cout << endl;
+            }
+            //cout << endl;
+          }
+        }
+      }
+    }
+    // 5. write final matrix to file
+    line = path;
+    line.append("/");
+    line += FRAGMENTPREFIX;
+    FragmentNumber = FixedDigitNumber(NumberOfMolecules, i);
+    line += FragmentNumber;
+    delete(FragmentNumber);
+    line += HCORRECTIONSUFFIX;
+    output.open(line.c_str());
+    output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
+    for (int j=0;j<b;j++) {
+      for(int i=0;i<a;i++)
+        output << correction[i][j] << "\t";
+      output << endl;
+    }
+    output.close();
+  }
+  line = path;
+  line.append("/");
+  line += HCORRECTIONSUFFIX;
+  output.open(line.c_str());
+  output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
+  for (int j=0;j<b;j++) {
+    for(int i=0;i<a;i++)
+      output << 0 << "\t";
+    output << endl;
+  }
+  output.close();
+  // 6. free memory of parsed matrices
+  FitConstant = (double ***) Malloc(sizeof(double **)*a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
+  for (int k=0;k<3;k++) {
+    FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
+    for (int i=a;i--;) {
+      FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
+    }
+  }
+  cout << "done." << endl;
+  return true;
+};
 
 /** Store force indices, i.e. the connection between the nuclear index in the total molecule config and the respective atom in fragment config.
@@ -224,5 +411,5 @@
     outputFragment.open(FragmentName, ios::out);
     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
-    if (intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment))
+    if ((intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment)))
       *out << " done." << endl;
     else
@@ -263,5 +450,5 @@
     configuration->SetDefaultPath(FragmentName);
     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
-    if (intermediateResult = configuration->Save(&outputFragment, ListOfMolecules[i]->elemente, ListOfMolecules[i]))
+    if ((intermediateResult = configuration->Save(&outputFragment, ListOfMolecules[i]->elemente, ListOfMolecules[i])))
       *out << " done." << endl;
     else
Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision c685eb762cafcd3ac83229798d1b5348383442c3)
+++ src/molecules.cpp	(revision 3902484c207d407476dea9c34dc0482f306011aa)
@@ -2707,4 +2707,7 @@
     StoreAdjacencyToFile(out, configuration->configpath);
 
+    // store Hydrogen saturation correction file
+    BondFragments->AddHydrogenCorrection(out, configuration->configpath);
+    
     // store adaptive orders into file
     StoreOrderAtSiteFile(out, configuration->configpath);
Index: src/molecules.hpp
===================================================================
--- src/molecules.hpp	(revision c685eb762cafcd3ac83229798d1b5348383442c3)
+++ src/molecules.hpp	(revision 3902484c207d407476dea9c34dc0482f306011aa)
@@ -289,4 +289,5 @@
 
   /// Output configs.
+  bool AddHydrogenCorrection(ofstream *out, char *path);
   bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
   bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
Index: src/parser.cpp
===================================================================
--- src/parser.cpp	(revision c685eb762cafcd3ac83229798d1b5348383442c3)
+++ src/parser.cpp	(revision 3902484c207d407476dea9c34dc0482f306011aa)
@@ -146,5 +146,5 @@
     file.str(" ");
     if (i != MatrixCounter)  
-      file << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix << suffix;
+      file << name << FRAGMENTPREFIX << FragmentNumber << prefix << suffix;
     else
       file << name << prefix << suffix;
@@ -265,4 +265,36 @@
 };
 
+/** Scans all elements of MatrixContainer::Matrix for greatest absolute value.
+ * \return greatest value of MatrixContainer::Matrix
+ */
+double MatrixContainer::FindMaxValue()
+{
+  double max = Matrix[0][0][0];
+  for(int i=MatrixCounter+1;i--;)
+    for(int j=RowCounter[i]+1;j--;)
+      for(int k=ColumnCounter;k--;)
+        if (fabs(Matrix[i][j][k]) > max)
+          max = fabs(Matrix[i][j][k]);
+  if (fabs(max) < MYEPSILON)
+    max += MYEPSILON;
+ return max;
+};
+
+/** Scans all elements of MatrixContainer::Matrix for smallest absolute value.
+ * \return smallest value of MatrixContainer::Matrix
+ */
+double MatrixContainer::FindMinValue()
+{
+  double min = Matrix[0][0][0];
+  for(int i=MatrixCounter+1;i--;)
+    for(int j=RowCounter[i]+1;j--;)
+      for(int k=ColumnCounter;k--;)
+        if (fabs(Matrix[i][j][k]) < min)
+          min = fabs(Matrix[i][j][k]);
+  if (fabs(min) < MYEPSILON)
+    min += MYEPSILON;
+  return min;
+};
+
 /** Sets all values in the last of MatrixContainer::Matrix to \a value.
  * \param value reset value
@@ -333,6 +365,6 @@
               }
             }
-          //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;
+            //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;
           }
         } else {
@@ -341,5 +373,5 @@
       }
     }
-    //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;
+   //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;
   }
   
@@ -428,15 +460,23 @@
  * Sums over the "F"-terms in ANOVA decomposition.
  * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
+ * \param CorrectionFragments MatrixContainer with hydrogen saturation correction per fragments
  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
  * \param Order bond order
+ * \parsm sign +1 or -1
  * \return true if summing was successful
  */
-bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, int Order)
+bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySet, int Order, double sign)
 {
   // sum energy
-  for(int i=KeySet.FragmentsPerOrder[Order];i--;)
-    for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
-      for(int k=ColumnCounter;k--;)
-        Matrix[MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
+  if (CorrectionFragments == NULL)
+    for(int i=KeySet.FragmentsPerOrder[Order];i--;)
+      for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
+        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;k--;)
+          Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
   return true;
 };
@@ -492,20 +532,26 @@
  * \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 ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order)
-{
+bool ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign)
+{
+  int FragmentNr;
   // sum forces
-  for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++)
-    for(int l=0;l<RowCounter[ KeySet.OrderSet[Order][i] ];l++) {
-      int j = Indices[ KeySet.OrderSet[Order][i] ][l];
+  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 force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
         return false;
       }
-      if (j != -1)
+      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++)
-          Matrix[MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[Order][i] ][l][k];
-    }
+          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
+      }
+    }
+  }
   return true;
 };
@@ -571,5 +617,5 @@
       KeySets[i][j] = -1;
     FragmentNumber = FixedDigitNumber(FragmentCounter, i);
-    cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
+    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
     Free((void **)&FragmentNumber, "KeySetsContainer::ParseKeySets: *FragmentNumber");
     input.getline(filename, 1023);
@@ -577,7 +623,7 @@
     for(int j=0;(j<AtomCounter[i]) && (!line.eof());j++) {
       line >> KeySets[i][j];
-      cout << " " << KeySets[i][j];
-    }
-    cout << endl;
+      //cout << " " << KeySets[i][j];
+    }
+    //cout << endl;
   }
   input.close();
Index: src/parser.hpp
===================================================================
--- src/parser.hpp	(revision c685eb762cafcd3ac83229798d1b5348383442c3)
+++ src/parser.hpp	(revision 3902484c207d407476dea9c34dc0482f306011aa)
@@ -19,4 +19,5 @@
 
 #define EnergySuffix ".energy.all"
+#define HcorrectionSuffix ".Hcorrection.all"
 #define ForcesSuffix ".forces.all"
 #define ShieldingSuffix ".sigma_all.csv"
@@ -26,4 +27,5 @@
 #define TimeSuffix ".speed"
 #define EnergyFragmentSuffix ".energyfragment.all"
+#define HcorrectionFragmentSuffix ".Hcorrectionfragment.all"
 #define ForceFragmentSuffix ".forcefragment.all"
 #define OrderSuffix ".Order"
@@ -52,4 +54,6 @@
   bool AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int CCounter);
   bool ResetMatrix();
+  double FindMinValue();
+  double FindMaxValue();
   bool SetLastMatrix(double value, int skipcolumns);
   bool SetLastMatrix(double **values, int skipcolumns);
@@ -66,5 +70,5 @@
   public:
     bool ParseIndices();
-    bool SumSubEnergy(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, int Order);
+    bool SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySet, int Order, double sign);
 };
 
@@ -74,5 +78,5 @@
   public: 
     bool ParseIndices(char *name);
-    bool SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order);
+    bool SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign);
 };
 
