Index: src/Fragmentation/analyzer.cpp
===================================================================
--- src/Fragmentation/analyzer.cpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
+++ src/Fragmentation/analyzer.cpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
@@ -0,0 +1,576 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
+/** \file analyzer.cpp
+ *
+ * Takes evaluated fragments (energy and forces) and does evaluation of how sensible the BOSSANOVA
+ * approach was, e.g. in the decay of the many-body-contributions.
+ *
+ */
+
+//============================ INCLUDES ===========================
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include <cstring>
+#include <sstream>
+#include <fstream>
+#include <cmath>
+
+#include "datacreator.hpp"
+#include "parser.hpp"
+#include "periodentafel.hpp"
+#include "CodePatterns/Log.hpp"
+#include "CodePatterns/Verbose.hpp"
+
+//============================== MAIN =============================
+
+int main(int argc, char **argv)
+{
+  periodentafel *periode = NULL; // and a period table of all elements
+  EnergyMatrix Energy;
+  EnergyMatrix EnergyFragments;
+  ForceMatrix Force;
+  ForceMatrix ForceFragments;
+  HessianMatrix Hessian;
+  HessianMatrix HessianFragments;
+  EnergyMatrix Hcorrection;
+  EnergyMatrix HcorrectionFragments;
+  ForceMatrix Shielding;
+  ForceMatrix ShieldingPAS;
+  ForceMatrix Chi;
+  ForceMatrix ChiPAS;
+  EnergyMatrix Time;
+  ForceMatrix ShieldingFragments;
+  ForceMatrix ShieldingPASFragments;
+  ForceMatrix ChiFragments;
+  ForceMatrix ChiPASFragments;
+  KeySetsContainer KeySet;
+  std::ofstream output;
+  std::ofstream output2;
+  std::ofstream output3;
+  std::ofstream output4;
+  std::ifstream input;
+  std::stringstream filename;
+  time_t t = time(NULL);
+  struct tm *ts = localtime(&t);
+  char *datum = asctime(ts);
+  std::stringstream Orderxrange;
+  std::stringstream Fragmentxrange;
+  std::stringstream yrange;
+  char *dir = NULL;
+  bool NoHessian = false;
+  bool NoTime = false;
+  bool NoHCorrection = true;
+  int counter = 0;
+
+  DoLog(0) && (Log() << Verbose(0) << "ANOVA Analyzer" << endl);
+  DoLog(0) && (Log() << Verbose(0) << "==============" << endl);
+
+  // Get the command line options
+  if (argc < 4) {
+    DoLog(0) && (Log() << Verbose(0) << "Usage: " << argv[0] << " <inputdir> <prefix> <outputdir> [elementsdb]" << endl);
+    DoLog(0) && (Log() << Verbose(0) << "<inputdir>\ttherein the output of a molecuilder fragmentation is expected, each fragment with a subdir containing an energy.all and a forces.all file." << endl);
+    DoLog(0) && (Log() << Verbose(0) << "<prefix>\tprefix of energy and forces file." << endl);
+    DoLog(0) && (Log() << Verbose(0) << "<outputdir>\tcreated plotfiles and datafiles are placed into this directory " << endl);
+    DoLog(0) && (Log() << Verbose(0) << "[elementsdb]\tpath to elements database, needed for shieldings." << endl);
+    return 1;
+  } else {
+    dir = new char[strlen(argv[2]) + 2];
+    strcpy(dir, "/");
+    strcat(dir, argv[2]);
+  }
+
+  if (argc > 4) {
+    DoLog(0) && (Log() << Verbose(0) << "Loading periodentafel." << endl);
+    periode = new periodentafel;
+    periode->LoadPeriodentafel(argv[4]);
+  }
+
+  // Test the given directory
+  if (!TestParams(argc, argv))
+    return 1;
+
+  // +++++++++++++++++ PARSING +++++++++++++++++++++++++++++++
+
+  // ------------- Parse through all Fragment subdirs --------
+  if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
+  if (!Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0)) {
+    NoHCorrection = true;
+    DoeLog(2) && (eLog()<< Verbose(2) << "No HCorrection file found, skipping these." << endl);
+  }
+  
+  if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
+  if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix,0,0)) {
+    NoHessian = true;
+    DoeLog(2) && (eLog()<< Verbose(2) << "No Hessian file found, skipping these." << endl);
+  }
+  if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) {
+    NoTime = true;
+    DoeLog(2) && (eLog()<< Verbose(2) << "No speed file found, skipping these." << endl);
+  }
+  if (periode != NULL) { // also look for PAS values
+    if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
+    if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
+    if (!Chi.ParseFragmentMatrix(argv[1], dir, ChiSuffix, 1, 0)) return 1;
+    if (!ChiPAS.ParseFragmentMatrix(argv[1], dir, ChiPASSuffix, 1, 0)) return 1;
+  }
+
+  // ---------- Parse the TE Factors into an array -----------------
+  if (!Energy.ParseIndices()) return 1;
+  if (!NoHCorrection) 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.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
+
+  // ---------- Parse hessian indices into an array -----------------
+  if (!NoHessian) {
+    if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
+    if (!HessianFragments.AllocateMatrix(Hessian.Header, Hessian.MatrixCounter, Hessian.RowCounter, Hessian.ColumnCounter)) return 1;
+    if (!HessianFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
+  }
+
+  // ---------- Parse the shielding indices into an array ---------------
+  if (periode != NULL) { // also look for PAS values
+    if(!Shielding.ParseIndices(argv[1])) return 1;
+    if(!ShieldingPAS.ParseIndices(argv[1])) return 1;
+    if (!ShieldingFragments.AllocateMatrix(Shielding.Header, Shielding.MatrixCounter, Shielding.RowCounter, Shielding.ColumnCounter)) return 1;
+    if (!ShieldingPASFragments.AllocateMatrix(ShieldingPAS.Header, ShieldingPAS.MatrixCounter, ShieldingPAS.RowCounter, ShieldingPAS.ColumnCounter)) return 1;
+    if(!ShieldingFragments.ParseIndices(argv[1])) return 1;
+    if(!ShieldingPASFragments.ParseIndices(argv[1])) return 1;
+    if(!Chi.ParseIndices(argv[1])) return 1;
+    if(!ChiPAS.ParseIndices(argv[1])) return 1;
+    if (!ChiFragments.AllocateMatrix(Chi.Header, Chi.MatrixCounter, Chi.RowCounter, Chi.ColumnCounter)) return 1;
+    if (!ChiPASFragments.AllocateMatrix(ChiPAS.Header, ChiPAS.MatrixCounter, ChiPAS.RowCounter, ChiPAS.ColumnCounter)) return 1;
+    if(!ChiFragments.ParseIndices(argv[1])) return 1;
+    if(!ChiPASFragments.ParseIndices(argv[1])) return 1;
+  }
+
+  // ---------- Parse the KeySets into an array ---------------
+  if (!KeySet.ParseKeySets(argv[1], Force.RowCounter, Force.MatrixCounter)) return 1;
+  if (!KeySet.ParseManyBodyTerms()) return 1;
+
+  // ---------- Parse fragment files created by 'joiner' into an array -------------
+  if (!EnergyFragments.ParseFragmentMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
+  if (!NoHCorrection) 
+    HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
+  if (!ForceFragments.ParseFragmentMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
+  if (!NoHessian)
+    if (!HessianFragments.ParseFragmentMatrix(argv[1], dir, HessianFragmentSuffix,0,0)) return 1;
+  if (periode != NULL) { // also look for PAS values
+    if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingFragmentSuffix, 1, 0)) return 1;
+    if (!ShieldingPASFragments.ParseFragmentMatrix(argv[1], dir, ShieldingPASFragmentSuffix, 1, 0)) return 1;
+    if (!ChiFragments.ParseFragmentMatrix(argv[1], dir, ChiFragmentSuffix, 1, 0)) return 1;
+    if (!ChiPASFragments.ParseFragmentMatrix(argv[1], dir, ChiPASFragmentSuffix, 1, 0)) return 1;
+  }
+
+  // +++++++++++++++ TESTING ++++++++++++++++++++++++++++++
+
+  // print energy and forces to file
+  filename.str("");
+  filename << argv[3] << "/" << "energy-forces.all";
+  output.open(filename.str().c_str(), ios::out);
+  output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header[Energy.MatrixCounter] << endl;
+  for(int j=0;j<Energy.RowCounter[Energy.MatrixCounter];j++) {
+    for(int k=0;k<Energy.ColumnCounter[Energy.MatrixCounter];k++)
+      output << scientific << Energy.Matrix[ Energy.MatrixCounter ][j][k] << "\t";
+    output << endl;
+  }
+  output << endl;
+  
+  output << endl << "Total Forces" << endl << "===============" << endl << Force.Header[Force.MatrixCounter] << endl;
+  for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
+    for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
+      output << scientific << Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
+    output << endl;
+  }
+  output << endl;
+
+  if (!NoHessian) {
+    output << endl << "Total Hessian" << endl << "===============" << endl << Hessian.Header[Hessian.MatrixCounter] << endl;
+    for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
+      for(int k=0;k<Hessian.ColumnCounter[Hessian.MatrixCounter];k++)
+        output << scientific << Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
+      output << endl;
+    }
+    output << endl;
+  }
+
+  if (periode != NULL) { // also look for PAS values
+    output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header[Hessian.MatrixCounter] << endl;
+    for(int j=0;j<Shielding.RowCounter[Shielding.MatrixCounter];j++) {
+      for(int k=0;k<Shielding.ColumnCounter[Shielding.MatrixCounter];k++)
+        output << scientific << Shielding.Matrix[ Shielding.MatrixCounter ][j][k] << "\t";
+      output << endl;
+    }
+    output << endl;
+  
+    output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header[ShieldingPAS.MatrixCounter] << endl;
+    for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
+      for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
+        output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t";
+      output << endl;
+    }
+    output << endl;
+
+    output << endl << "Total Chis" << endl << "===============" << endl << Chi.Header[Chi.MatrixCounter] << endl;
+    for(int j=0;j<Chi.RowCounter[Chi.MatrixCounter];j++) {
+      for(int k=0;k<Chi.ColumnCounter[Chi.MatrixCounter];k++)
+        output << scientific << Chi.Matrix[ Chi.MatrixCounter ][j][k] << "\t";
+      output << endl;
+    }
+    output << endl;
+  
+    output << endl << "Total Chis PAS" << endl << "===============" << endl << ChiPAS.Header[ChiPAS.MatrixCounter] << endl;
+    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
+      for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
+        output << scientific << ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t";
+      output << endl;
+    }
+    output << endl;
+  }
+  
+  if (!NoTime) {
+    output << endl << "Total Times" << endl << "===============" << endl << Time.Header[Time.MatrixCounter] << endl;
+    for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
+      for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
+        output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
+      }
+      output << endl;
+    }
+    output << endl;
+  }
+  output.close();
+  if (!NoTime)
+    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];
+
+  // +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++
+
+  DoLog(0) && (Log() << Verbose(0) << "Analyzing ..." << endl);
+
+  // ======================================= Creating the data files ==============================================================
+
+  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
+  // +++++++++++++++++++++++++++++++++++++++ Plotting Delta Simtime vs Bond Order
+  if (!NoTime) {
+    if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
+    if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
+    for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
+      for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
+        Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
+      }
+    counter = 0;
+    output << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl;
+    output2 << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << 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[Time.MatrixCounter];k--;) {
+            Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
+          }
+      counter += KeySet.FragmentsPerOrder[BondOrder];
+      output << BondOrder+1 << "\t" << counter;
+      output2 << BondOrder+1 << "\t" << counter;
+      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)
+          output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
+        else
+          output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
+      }
+      output << endl;
+      output2 << endl;
+    }
+    output.close();
+    output2.close();
+  }
+
+  if (!NoHessian) {
+    // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in hessian to full QM
+    if (!CreateDataDeltaHessianOrderPerAtom(Hessian, HessianFragments, KeySet, argv[3], "DeltaHessian_xx-Order", "Plot of error between approximated hessian and full hessian versus the Bond Order", datum)) return 1;
+
+    if (!CreateDataDeltaFrobeniusOrderPerAtom(Hessian, HessianFragments, KeySet, argv[3], "DeltaFrobeniusHessian_xx-Order", "Plot of error between approximated hessian and full hessian in the frobenius norm versus the Bond Order", datum)) return 1;
+
+    // ++++++++++++++++++++++++++++++++++++++Plotting Hessian vs. Order
+    if (!CreateDataHessianOrderPerAtom(HessianFragments, KeySet, argv[3], "Hessian_xx-Order", "Plot of approximated hessian versus the Bond Order", datum)) return 1;
+    if (!AppendOutputFile(output, argv[3], "Hessian_xx-Order.dat" )) return false;
+    output << endl << "# Full" << endl;
+    for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
+      output << j << "\t";
+      for(int k=0;k<Hessian.ColumnCounter[Force.MatrixCounter];k++)
+        output << scientific <<  Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
+      output << endl;
+    }
+    output.close();
+  }
+
+  // +++++++++++++++++++++++++++++++++++++++ Plotting shieldings
+  if (periode != NULL) { // also look for PAS values
+    if (!CreateDataDeltaForcesOrderPerAtom(ShieldingPAS, ShieldingPASFragments, KeySet, argv[3], "DeltaShieldingsPAS-Order", "Plot of error between approximated shieldings and full shieldings versus the Bond Order", datum)) return 1;
+    if (!CreateDataForcesOrderPerAtom(ShieldingPASFragments, KeySet, argv[3], "ShieldingsPAS-Order", "Plot of approximated shieldings versus the Bond Order", datum)) return 1;
+    if (!AppendOutputFile(output, argv[3], "ShieldingsPAS-Order.dat" )) return false;
+    output << endl << "# Full" << endl;
+    for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
+      output << j << "\t";
+      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;
+    }
+    output.close();
+    if (!CreateDataDeltaForcesOrderPerAtom(ChiPAS, ChiPASFragments, KeySet, argv[3], "DeltaChisPAS-Order", "Plot of error between approximated Chis and full Chis versus the Bond Order", datum)) return 1;
+    if (!CreateDataForcesOrderPerAtom(ChiPASFragments, KeySet, argv[3], "ChisPAS-Order", "Plot of approximated Chis versus the Bond Order", datum)) return 1;
+    if (!AppendOutputFile(output, argv[3], "ChisPAS-Order.dat" )) return false;
+    output << endl << "# Full" << endl;
+    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
+      output << j << "\t";
+      for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
+        output << scientific <<  ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
+      output << endl;
+    }
+    output.close();
+  }
+
+
+  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
+  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 (!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 (!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 (!CreateDataForcesOrderPerAtom(ForceFragments, KeySet, argv[3], "Forces-Order", "Plot of approximated forces versus the Bond Order", datum)) return 1;
+  if (!AppendOutputFile(output, argv[3], "Forces-Order.dat" )) return false;
+  output << endl << "# Full" << endl;
+  for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
+    output << j << "\t";
+    for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
+      output << scientific <<  Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
+    output << endl;
+  }
+  output.close();
+  // 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(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
+  if (!CreateDataFragment(EnergyFragments, KeySet, argv[3], "Energies-Fragment", "Plot of fragment energy versus the Fragment No", datum, CreateEnergy)) return 1;
+  if (!CreateDataFragment(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", "Plot of fragment energy of each Fragment No vs. Bond Order", datum, CreateEnergy)) return 1;
+  if (!CreateDataFragmentOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", "Plot of maximum of fragment energy vs. Bond Order", datum, CreateMaxFragmentOrder)) return 1;
+  if (!CreateDataFragmentOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", "Plot of minimum of fragment energy vs. Bond Order", datum, CreateMinFragmentOrder)) return 1;
+
+  // +++++++++++++++++++++++++++++++Ploting min/mean/max forces for each fragment
+  // min force
+  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MinForces-Fragment", "Plot of min approximated forces versus the Fragment No", datum, CreateMinimumForce)) return 1;
+
+  // mean force
+  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MeanForces-Fragment", "Plot of mean approximated forces versus the Fragment No", datum, CreateMeanForce)) return 1;
+
+  // max force
+  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MaxForces-Fragment", "Plot of max approximated forces versus the Fragment No", datum, CreateMaximumForce)) return 1;
+
+  // +++++++++++++++++++++++++++++++Ploting min/mean/max forces for each fragment per order
+  // min force
+  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MinForces-FragmentOrder", "Plot of min approximated forces of each Fragment No vs. Bond Order", datum, CreateMinimumForce)) return 1;
+
+  // mean force
+  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MeanForces-FragmentOrder", "Plot of mean approximated forces of each Fragment No vs. Bond Order", datum, CreateMeanForce)) return 1;
+
+  // max force
+  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MaxForces-FragmentOrder", "Plot of max approximated forces of each Fragment No vs. Bond Order", datum, CreateMaximumForce)) return 1;
+
+  // ======================================= Creating the plot files ==============================================================
+
+  Orderxrange << "[1:" << KeySet.Order << "]";
+  Fragmentxrange << "[0:" << KeySet.FragmentCounter+1 << "]";
+  yrange.str("[1e-8:1e+1]");
+  
+  if (!NoTime) {
+    // +++++++++++++++++++++++++++++++++++++++ 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;
+  }
+  
+  // +++++++++++++++++++++++++++++++++++++++ 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(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
+
+  // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
+  if (!CreatePlotOrder(Energy, KeySet, argv[3], "Energies-Order", 1, "outside", "", "",  1, 1, "bond order k", "approximate energy [Ht]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1;
+
+  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
+  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
+  if(!OpenOutputFile(output, argv[3], "DeltaMinMeanMaxTotalForce-Order.pyx")) return 1;
+  CreatePlotHeader(output, "DeltaMinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL,  1, 1, "bond order k", "absolute error in total forces [Ht/a.u.]");
+  output << "plot " << Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;
+  output << "'DeltaMinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
+  output << "'DeltaMeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
+  output << "'DeltaMaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
+  output.close();
+
+  // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
+  // min force
+  if (!CreatePlotOrder(Force, KeySet, argv[3], "MinForces-Order", 2, "bottom right", "y", "",  1, 1, "bond order k", "absolute approximated min force [Ht/a.u.]", Orderxrange.str().c_str(), "", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
+
+  // mean force
+  if (!CreatePlotOrder(Force, KeySet, argv[3], "MeanForces-Order", 2, "bottom right", "y", "",  1, 1, "bond order k", "absolute approximated mean force [Ht/a.u.]", Orderxrange.str().c_str(), "", "1" , "with linespoints", AbsFirstForceValuePlotLine)) return 1;
+
+  // max force
+  if (!CreatePlotOrder(Force, KeySet, argv[3], "MaxForces-Order", 2, "bottom right", "y", "",  1, 1, "bond order k", "absolute approximated max force [Ht/a.u.]", Orderxrange.str().c_str(), "", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
+
+  // min/mean/max comparison for total force
+  if(!OpenOutputFile(output, argv[3],"MinMeanMaxTotalForce-Order.pyx")) return 1;
+  CreatePlotHeader(output, "MinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL, 1, 1, "bond order k", "absolute total force [Ht/a.u.]");
+  output << "plot "<< Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;
+  output << "'MinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
+  output << "'MeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
+  output << "'MaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
+  output.close();
+
+  // ++++++++++++++++++++++++++++++++++++++Plotting vector sum vs. Order
+
+  if (!CreatePlotOrder(Force, KeySet, argv[3], "VectorSum-Order", 2, "bottom right", "y" ,"", 1, 1, "bond order k", "vector sum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), "", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
+
+  // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
+  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(), 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(), 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(), 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(), 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(), 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(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
+
+  // +++++++++++++++++++++++++++++++=Ploting approximated and true shielding for each atom
+  if (periode != NULL) { // also look for PAS values
+    if(!OpenOutputFile(output, argv[3], "ShieldingsPAS-Order.pyx")) return 1;
+    if(!OpenOutputFile(output2, argv[3], "DeltaShieldingsPAS-Order.pyx")) return 1;
+    CreatePlotHeader(output, "ShieldingsPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical shielding value [ppm]");
+    CreatePlotHeader(output2, "DeltaShieldingsPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical shielding value [ppm]");
+    double step=0.8/KeySet.Order;
+    output << "set boxwidth " << step << endl;
+    output << "plot [0:" << ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter]+10 << "]\\" << endl;
+    output2 << "plot [0:" << ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter]+10 << "]\\" << endl;
+    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+      output << "'ShieldingsPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
+      output2 << "'DeltaShieldingsPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
+      if (BondOrder-1 != KeySet.Order)
+        output2 << ", \\" << endl;
+    }
+    output << "'ShieldingsPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
+    output2.close();  
+
+    if(!OpenOutputFile(output, argv[3], "ChisPAS-Order.pyx")) return 1;
+    if(!OpenOutputFile(output2, argv[3], "DeltaChisPAS-Order.pyx")) return 1;
+    CreatePlotHeader(output, "ChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
+    CreatePlotHeader(output2, "DeltaChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
+    output << "set boxwidth " << step << endl;
+    output << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
+    output2 << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
+    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+      output << "'ChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
+      output2 << "'DeltaChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
+      if (BondOrder-1 != KeySet.Order)
+        output2 << ", \\" << endl;
+    }
+    output << "'ChisPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
+    output.close();  
+    output2.close();  
+
+    if(!OpenOutputFile(output, argv[3], "ChisPAS-Order.pyx")) return 1;
+    if(!OpenOutputFile(output2, argv[3], "DeltaChisPAS-Order.pyx")) return 1;
+    CreatePlotHeader(output, "ChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
+    CreatePlotHeader(output2, "DeltaChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
+    output << "set boxwidth " << step << endl;
+    output << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
+    output2 << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
+    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+      output << "'ChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
+      output2 << "'DeltaChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
+      if (BondOrder-1 != KeySet.Order)
+        output2 << ", \\" << endl;
+    }
+    output << "'ChisPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
+    output.close();  
+    output2.close();  
+  }
+
+  // create Makefile
+  if(!OpenOutputFile(output, argv[3], "Makefile")) return 1;
+  output << "PYX = $(shell ls *.pyx)" << endl << endl;
+  output << "EPS = $(PYX:.pyx=.eps)" << endl << endl;
+  output << "%.eps: %.pyx" << endl;
+  output << "\t~/build/pyxplot/pyxplot $<" << endl << endl;
+  output << "all: $(EPS)" << endl << endl;
+  output << ".PHONY: clean" << endl;
+  output << "clean:" << endl;
+  output << "\trm -rf $(EPS)" << endl;
+  output.close();
+
+  // ++++++++++++++++ exit ++++++++++++++++++++++++++++++++++
+  delete(periode);
+  delete[](dir);
+  DoLog(0) && (Log() << Verbose(0) << "done." << endl);
+  return 0;
+};
+
+//============================ END ===========================
Index: src/Fragmentation/datacreator.cpp
===================================================================
--- src/Fragmentation/datacreator.cpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
+++ src/Fragmentation/datacreator.cpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
@@ -0,0 +1,848 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
+/** \file datacreator.cpp
+ *
+ * Declarations of assisting functions in creating data and plot files.
+ *
+ */
+
+//============================ INCLUDES ===========================
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+
+#include "CodePatterns/Log.hpp"
+#include "CodePatterns/Verbose.hpp"
+#include "LinearAlgebra/defs.hpp"
+#include "Helpers/defs.hpp"
+#include "datacreator.hpp"
+#include "parser.hpp"
+
+using namespace std;
+
+//=========================== FUNCTIONS============================
+
+/** Opens a file with \a *filename in \a *dir.
+ * \param output file handle on return
+ * \param *dir directory
+ * \param *filename name of file
+ * \return true if file has been opened
+ */
+bool OpenOutputFile(ofstream &output, const char *dir, const char *filename)
+{
+  stringstream name;
+  name << dir << "/" << filename;
+  output.open(name.str().c_str(), ios::out);
+  if (output == NULL) {
+    DoLog(0) && (Log() << Verbose(0) << "Unable to open " << name.str() << " for writing, is directory correct?" << endl);
+    return false;
+  }
+  return true;
+};
+
+/** Opens a file for appending with \a *filename in \a *dir.
+ * \param output file handle on return
+ * \param *dir directory
+ * \param *filename name of file
+ * \return true if file has been opened
+ */
+bool AppendOutputFile(ofstream &output, const char *dir, const char *filename)
+{
+  stringstream name;
+  name << dir << "/" << filename;
+  output.open(name.str().c_str(), ios::app);
+  if (output == NULL) {
+    DoLog(0) && (Log() << Verbose(0) << "Unable to open " << name.str() << " for writing, is directory correct?" << endl);
+    return false;
+  }
+  return true;
+};
+
+/** Plots an energy vs. order.
+ * \param &Fragments EnergyMatrix class containing matrix values
+ * \param KeySets 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 CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
+  DoLog(0) && (Log() << Verbose(0) << msg << endl);
+  output << "# " << msg << ", created on " << datum;
+  output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    for(int i=KeySets.FragmentsPerOrder[BondOrder];i--;) {
+      for(int j=Fragments.RowCounter[ KeySets.OrderSet[BondOrder][i] ];j--;)
+        for(int k=Fragments.ColumnCounter[ KeySets.OrderSet[BondOrder][i] ];k--;)
+          Fragments.Matrix[Fragments.MatrixCounter][j][k] += Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
+    }
+    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
+    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;
+  }
+  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 KeySets 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 &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
+  DoLog(0) && (Log() << Verbose(0) << msg << endl);
+  output << "# " << msg << ", created on " << datum;
+  output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
+  Fragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0);
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    for(int i=KeySets.FragmentsPerOrder[BondOrder];i--;) {
+      for(int j=Fragments.RowCounter[ KeySets.OrderSet[BondOrder][i] ];j--;)
+        for(int k=Fragments.ColumnCounter[ KeySets.OrderSet[BondOrder][i] ];k--;)
+          Fragments.Matrix[Fragments.MatrixCounter][j][k] -= Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
+    }
+    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
+    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];
+      else
+        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;
+  }
+  output.close();
+  return true;
+};
+
+/** Plot forces vs. order.
+ * \param &Fragments ForceMatrix class containing matrix values
+ * \param KeySets 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 &KeySets, const char *dir, const char *prefix,const  char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int))
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
+  DoLog(0) && (Log() << Verbose(0) << msg << endl);
+  output << "# " << msg << ", created on " << datum;
+  output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
+  Fragments.SetLastMatrix(0.,0);
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    Fragments.SumSubForces(Fragments, KeySets, BondOrder, 1.);
+    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
+    CreateForce(Fragments, Fragments.MatrixCounter);
+    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];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 KeySets 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 &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int))
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
+  DoLog(0) && (Log() << Verbose(0) << msg << endl);
+  output << "# " << msg << ", created on " << datum;
+  output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
+  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0);
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    Fragments.SumSubForces(Fragments, KeySets, BondOrder, -1.);
+    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
+    CreateForce(Fragments, Fragments.MatrixCounter);
+    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];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 KeySets 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 &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
+  stringstream filename;
+  ofstream output;
+  double norm = 0.;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
+  DoLog(0) && (Log() << Verbose(0) << msg << endl);
+  output << "# " << msg << ", created on " << datum;
+  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
+  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter], 0);
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    //Log() << Verbose(0) << "Current order is " << BondOrder << "." << endl;
+    Fragments.SumSubForces(Fragments, KeySets, BondOrder, -1.);
+    // errors per atom
+    output << endl << "#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[ Fragments.MatrixCounter ];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 << endl;
+    }
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+/** Plot forces error vs. vs atom vs. order.
+ * \param &Fragments ForceMatrix class containing matrix values
+ * \param KeySets 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 &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
+  DoLog(0) && (Log() << Verbose(0) << msg << endl);
+  output << "# " << msg << ", created on " << datum;
+  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    //Log() << Verbose(0) << "Current order is " << BondOrder << "." << endl;
+    Fragments.SumSubForces(Fragments, KeySets, BondOrder, 1.);
+    // errors per atom
+    output << endl << "#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[ Fragments.MatrixCounter ];l++)
+        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
+      output << endl;
+    }
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+
+/** Plot hessian error vs. vs atom vs. order.
+ * \param &Hessian HessianMatrix containing reference values (in MatrixCounter matrix)
+ * \param &Fragments HessianMatrix class containing matrix values
+ * \param KeySets 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 CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
+  DoLog(0) && (Log() << Verbose(0) << msg << endl);
+  output << "# " << msg << ", created on " << datum;
+  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
+  Fragments.SetLastMatrix(Hessian.Matrix[Hessian.MatrixCounter], 0);
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    //Log() << Verbose(0) << "Current order is " << BondOrder << "." << endl;
+    Fragments.SumSubHessians(Fragments, KeySets, BondOrder, -1.);
+    // errors per atom
+    output << endl << "#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[ Fragments.MatrixCounter ];l++) {
+        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
+      }
+      output << endl;
+    }
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+/** Plot hessian error vs. vs atom vs. order in the frobenius norm.
+ * \param &Hessian HessianMatrix containing reference values (in MatrixCounter matrix)
+ * \param &Fragments HessianMatrix class containing matrix values
+ * \param KeySets 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 CreateDataDeltaFrobeniusOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
+  stringstream filename;
+  ofstream output;
+  double norm = 0;
+  double tmp;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
+  DoLog(0) && (Log() << Verbose(0) << msg << endl);
+  output << "# " << msg << ", created on " << datum;
+  output << "# AtomNo\t";
+  Fragments.SetLastMatrix(Hessian.Matrix[Hessian.MatrixCounter], 0);
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    output << "Order" << BondOrder+1 << "\t";
+  }
+  output << endl;
+  output << Fragments.RowCounter[ Fragments.MatrixCounter ] << "\t";
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    //Log() << Verbose(0) << "Current order is " << BondOrder << "." << endl;
+    Fragments.SumSubHessians(Fragments, KeySets, BondOrder, -1.);
+    // frobenius norm of errors per atom
+    norm = 0.;
+    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) { 
+      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
+        tmp = Fragments.Matrix[Fragments.MatrixCounter][ j ][l];
+        norm += tmp*tmp;
+      }
+    }
+    output << scientific << sqrt(norm)/(Fragments.RowCounter[ Fragments.MatrixCounter ]*Fragments.ColumnCounter[ Fragments.MatrixCounter] ) << "\t";
+  }
+  output << endl;
+  output.close();
+  return true;
+};
+
+/** Plot hessian error vs. vs atom vs. order.
+ * \param &Fragments HessianMatrix class containing matrix values
+ * \param KeySets 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 CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
+  DoLog(0) && (Log() << Verbose(0) << msg << endl);
+  output << "# " << msg << ", created on " << datum;
+  output << "# AtomNo\t" << Fragments.Header[ Fragments.MatrixCounter ] << endl;
+  Fragments.SetLastMatrix(0., 0);
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    //Log() << Verbose(0) << "Current order is " << BondOrder << "." << endl;
+    Fragments.SumSubHessians(Fragments, KeySets, BondOrder, 1.);
+    // errors per atom
+    output << endl << "#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[ Fragments.MatrixCounter ];l++)
+        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
+      output << endl;
+    }
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+/** Plot matrix vs. fragment.
+ */
+bool CreateDataFragment(class MatrixContainer &Fragment, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragment)(class MatrixContainer &, int))
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
+  DoLog(0) && (Log() << Verbose(0) << msg << endl);
+  output << "# " << msg << ", created on " << datum << endl;
+  output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl;
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    for(int i=0;i<KeySets.FragmentsPerOrder[BondOrder];i++) {
+      output << BondOrder+1 << "\t" << KeySets.OrderSet[BondOrder][i]+1;
+      CreateFragment(Fragment, KeySets.OrderSet[BondOrder][i]);
+      for (int l=0;l<Fragment.ColumnCounter[ KeySets.OrderSet[BondOrder][i] ];l++)
+        output << scientific << "\t" << Fragment.Matrix[ KeySets.OrderSet[BondOrder][i] ][ Fragment.RowCounter[ KeySets.OrderSet[BondOrder][i] ] ][l];
+      output << endl;
+    }
+  }
+  output.close();
+  return true;
+};
+
+/** Copies fragment energy values into last matrix of \a Matrix with greatest total energy.
+ * \param &Matrix MatrixContainer with all fragment energy values
+ * \param &KeySets KeySetsContainer with associations of each fragment to a bond order
+ * \param BondOrder current bond order
+ */
+void CreateMaxFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySets, int BondOrder)
+{
+  for(int j=Fragments.RowCounter[ Fragments.MatrixCounter ];j--;) {
+    for(int i=KeySets.FragmentsPerOrder[BondOrder];i--;) {
+      if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < fabs(Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][1])) {
+        for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
+          Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
+      }
+    }
+  }
+};
+
+/** Copies fragment energy values into last matrix of \a Matrix with smallest total energy.
+ * \param &Matrix MatrixContainer with all fragment energy values
+ * \param &KeySets KeySetsContainer with associations of each fragment to a bond order
+ * \param BondOrder current bond order
+ */
+void CreateMinFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySets, int BondOrder)
+{
+  for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
+    int i=0;
+    do {  // first get a minimum value unequal to 0
+      for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
+        Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
+      i++;
+    } while ((fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < MYEPSILON) && (i<KeySets.FragmentsPerOrder[BondOrder]));
+    for(;i<KeySets.FragmentsPerOrder[BondOrder];i++) { // then find lowest
+      if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) > fabs(Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][1])) {
+        for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
+          Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
+      }
+    }
+  }
+};
+
+/** Plot matrix vs. fragment.
+ */
+bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int))
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
+  DoLog(0) && (Log() << Verbose(0) << msg << endl);
+  output << "# " << msg << ", created on " << datum;
+  output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl;
+  // max
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    Fragment.SetLastMatrix(0.,0);
+    CreateFragmentOrder(Fragment, KeySets, BondOrder);
+    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
+    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;
+  }
+  output.close();
+  return true;
+};
+
+/** Takes last but one row and copies into final row.
+ * \param Energy MatrixContainer with matrix values
+ * \param MatrixNumber the index for the ForceMatrix::matrix array
+ */
+void CreateEnergy(class MatrixContainer &Energy, int MatrixNumber)
+{
+  for(int k=0;k<Energy.ColumnCounter[MatrixNumber];k++)
+    Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber] ] [k] =  Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber]-1 ] [k];
+};
+
+/** Scans forces for the minimum in magnitude.
+ * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
+ * \param Force ForceMatrix class containing matrix values
+ * \param MatrixNumber the index for the ForceMatrix::matrix array
+ */
+void CreateMinimumForce(class MatrixContainer &Force, int MatrixNumber)
+{
+  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
+    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
+  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
+    double stored = 0;
+    int k=0;
+    do {
+      for (int m=NDIM;m--;) {
+        stored += Force.Matrix[MatrixNumber][ k ][l+m]
+              * Force.Matrix[MatrixNumber][ k ][l+m];
+        Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m]  = Force.Matrix[MatrixNumber][ k ][l+m];
+      }
+      k++;
+    } while ((fabs(stored) < MYEPSILON) && (k<Force.RowCounter[MatrixNumber]));
+    for (;k<Force.RowCounter[MatrixNumber];k++) {
+      double tmp = 0;
+      for (int m=NDIM;m--;)
+        tmp += Force.Matrix[MatrixNumber][ k ][l+m]
+              * Force.Matrix[MatrixNumber][ k ][l+m];
+      if ((fabs(tmp) > MYEPSILON) && (tmp < stored)) {  // current force is greater than stored
+        for (int m=NDIM;m--;)
+          Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m]  = Force.Matrix[MatrixNumber][ k ][l+m];
+        stored = tmp;
+      }
+    }
+  }
+};
+
+/** Scans forces for the mean in magnitude.
+ * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
+ * \param Force ForceMatrix class containing matrix values
+  * \param MatrixNumber the index for the ForceMatrix::matrix array
+ */
+void CreateMeanForce(class MatrixContainer &Force, int MatrixNumber)
+{
+  int divisor = 0;
+  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
+    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
+  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
+    double tmp = 0;
+    for (int k=Force.RowCounter[MatrixNumber];k--;) {
+      double norm = 0.;
+      for (int m=NDIM;m--;)
+        norm += Force.Matrix[MatrixNumber][ k ][l+m]
+              * Force.Matrix[MatrixNumber][ k ][l+m];
+      tmp += sqrt(norm);
+      if (fabs(norm) > MYEPSILON) divisor++;
+    }
+    tmp /= (double)divisor;
+    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = tmp;
+  }
+};
+
+/** Scans forces for the maximum in magnitude.
+ * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
+ * \param Force ForceMatrix class containing matrix values
+ * \param MatrixNumber the index for the ForceMatrix::matrix array
+ */
+void CreateMaximumForce(class MatrixContainer &Force, int MatrixNumber)
+{
+  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
+    double stored = 0;
+    for (int k=Force.RowCounter[MatrixNumber];k--;) {
+      double tmp = 0;
+      for (int m=NDIM;m--;)
+        tmp += Force.Matrix[MatrixNumber][ k ][l+m]
+              * Force.Matrix[MatrixNumber][ k ][l+m];
+      if (tmp > stored) {  // current force is greater than stored
+        for (int m=NDIM;m--;)
+          Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m]  = Force.Matrix[MatrixNumber][ k ][l+m];
+        stored = tmp;
+      }
+    }
+  }
+};
+
+/** 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.
+ * \param Force ForceMatrix class containing matrix values
+ * \param MatrixNumber the index for the ForceMatrix::matrix array
+ */
+void CreateVectorSumForce(class MatrixContainer &Force, int MatrixNumber)
+{
+  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
+    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
+  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];
+  }
+};
+
+/** Writes the standard pyxplot header info.
+ * \param keycolumns number of columns of the key
+ * \param *key position of key
+ * \param *logscale axis for logscale
+ * \param *extraline extra set lines if desired
+ * \param mxtics small tics at ...
+ * \param xtics large tics at ...
+ * \param *xlabel label for x axis
+ * \param *ylabel label for y axis
+ */
+void CreatePlotHeader(ofstream &output, const char *prefix, const int keycolumns, const char *key, const char *logscale, const char *extraline, const int mxtics, const int xtics, const char *xlabel, const char *ylabel)
+{
+  //output << "#!/home/heber/build/pyxplot/pyxplot" << endl << endl;
+  output << "reset" << endl;
+  output << "set keycolumns "<< keycolumns << endl;
+  output << "set key " << key << endl;
+  output << "set mxtics "<< mxtics << endl;
+  output << "set xtics "<< xtics << endl;
+  if (logscale != NULL)
+    output << "set logscale " << logscale << endl;
+  if (extraline != NULL)
+    output << extraline << endl;
+  output << "set xlabel '" << xlabel << "'" << endl;
+  output << "set ylabel '" << ylabel << "'" << endl;
+  output << "set terminal eps color" << endl;
+  output << "set output '"<< prefix << ".eps'" << endl;
+};
+
+/** Creates the pyxplotfile for energy data.
+ * \param Matrix MatrixContainer with matrix values
+ * \param KeySets contains bond order
+ * \param *dir directory
+ * \param *prefix prefix for all filenames (without ending)
+ * \param keycolumns number of columns of the key
+ * \param *key position of key
+ * \param logscale axis for logscale
+ * \param mxtics small tics at ...
+ * \param xtics large tics at ...
+ * \param xlabel label for x axis
+ * \param xlabel label for x axis
+ * \param *xrange xrange
+ * \param *yrange yrange
+ * \param *xargument number of column to plot against (x axis)
+ * \param uses using string for plot command
+ * \param (*CreatePlotLines) function reference that writes a single plot line
+ * \return true if file was written successfully
+ */
+bool CreatePlotOrder(class MatrixContainer &Matrix, const class KeySetsContainer &KeySets, const char *dir, const char *prefix, const int keycolumns, const char *key, const char *logscale, const char *extraline, const int mxtics, const int xtics, const char *xlabel, const char *ylabel, const char *xrange, const char *yrange, const char *xargument, const char *uses, void (*CreatePlotLines)(ofstream &, class MatrixContainer &, const char *, const char *, const char *))
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".pyx";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
+  CreatePlotHeader(output, prefix, keycolumns, key, logscale, extraline, mxtics, xtics, xlabel, ylabel);
+  output << "plot " << xrange << " " << yrange << " \\" << endl;
+  CreatePlotLines(output, Matrix, prefix, xargument, uses);
+  output.close();
+  return true;
+};
+
+/** Writes plot lines for absolute energies.
+ * \param output file handler
+ * \param Energy MatrixContainer with matrix values
+ * \param *prefix prefix of data file
+ * \param *xargument number of column to plot against (x axis)
+ * \param *uses uses command
+ */
+void AbsEnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses)
+{
+  stringstream line(Energy.Header[ Energy.MatrixCounter ]);
+  string token;
+
+  getline(line, token, '\t');
+  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[Energy.MatrixCounter]))
+      output << ", \\";
+    output << endl;
+  }
+};
+
+/** Writes plot lines for energies.
+ * \param output file handler
+ * \param Energy MatrixContainer with matrix values
+ * \param *prefix prefix of data file
+ * \param *xargument number of column to plot against (x axis)
+ * \param *uses uses command
+ */
+void EnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses)
+{
+  stringstream line(Energy.Header[Energy.MatrixCounter]);
+  string token;
+
+  getline(line, token, '\t');
+  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[Energy.MatrixCounter]))
+      output << ", \\";
+    output << endl;
+  }
+};
+
+/** Writes plot lines for absolute force magnitudes.
+ * \param output file handler
+ * \param Force MatrixContainer with matrix values
+ * \param *prefix prefix of data file
+ * \param *xargument number of column to plot against (x axis)
+ * \param *uses uses command
+ */
+void ForceMagnitudePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
+{
+  stringstream line(Force.Header[Force.MatrixCounter]);
+  string token;
+
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
+    getline(line, token, '\t');
+    while (token[0] == ' ') // remove leading white spaces
+      token.erase(0,1);
+    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[Force.MatrixCounter]-1))
+      output << ", \\";
+    output << endl;
+    getline(line, token, '\t');
+    getline(line, token, '\t');
+  }
+};
+
+/** Writes plot lines for first component of force vector.
+ * \param output file handler
+ * \param Force MatrixContainer with matrix values
+ * \param *prefix prefix of data file
+ * \param *xargument number of column to plot against (x axis)
+ * \param *uses uses command
+ */
+void AbsFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
+{
+  stringstream line(Force.Header[Force.MatrixCounter]);
+  string token;
+
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
+    getline(line, token, '\t');
+    while (token[0] == ' ') // remove leading white spaces
+      token.erase(0,1);
+    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[Force.MatrixCounter]-1))
+      output << ", \\";
+    output << endl;
+    getline(line, token, '\t');
+    getline(line, token, '\t');
+  }
+};
+
+/** Writes plot lines for force vector as boxes with a fillcolor.
+ * \param output file handler
+ * \param Force MatrixContainer with matrix values
+ * \param *prefix prefix of data file
+ * \param *xargument number of column to plot against (x axis)
+ * \param *uses uses command
+ */
+void BoxesForcePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
+{
+  stringstream line(Force.Header[Force.MatrixCounter]);
+  const char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
+  string token;
+
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
+    getline(line, token, '\t');
+    while (token[0] == ' ') // remove leading white spaces
+      token.erase(0,1);
+    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[Force.MatrixCounter]-1))
+      output << ", \\";
+    output << endl;
+    getline(line, token, '\t');
+    getline(line, token, '\t');
+  }
+};
+
+/** Writes plot lines for first force vector component as boxes with a fillcolor.
+ * \param output file handler
+ * \param Force MatrixContainer with matrix values
+ * \param *prefix prefix of data file
+ * \param *xargument number of column to plot against (x axis)
+ * \param *uses uses command
+ */
+void BoxesFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
+{
+  stringstream line(Force.Header[Force.MatrixCounter]);
+  const char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
+  string token;
+
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  getline(line, token, '\t');
+  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
+    getline(line, token, '\t');
+    while (token[0] == ' ') // remove leading white spaces
+      token.erase(0,1);
+    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[Force.MatrixCounter]-1))
+      output << ", \\";
+    output << endl;
+    getline(line, token, '\t');
+    getline(line, token, '\t');
+  }
+};
Index: src/Fragmentation/datacreator.hpp
===================================================================
--- src/Fragmentation/datacreator.hpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
+++ src/Fragmentation/datacreator.hpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
@@ -0,0 +1,64 @@
+/** \file datacreator.cpp
+ *
+ * Definitions of assisting functions in creating data and plot files.
+ *
+ */
+
+#ifndef DATACREATOR_HPP_
+#define DATACREATOR_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+using namespace std;
+
+#include <iosfwd>
+
+/****************************************** forward declarations *****************************/
+
+class EnergyMatrix;
+class ForceMatrix;
+class HessianMatrix;
+class KeySetsContainer;
+class MatrixContainer;
+
+/********************************************** declarations *******************************/
+
+bool OpenOutputFile(ofstream &output, const char *dir, const char *filename);
+bool AppendOutputFile(ofstream &output, const char *dir, const char *filename);
+
+bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum);
+bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum);
+bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int));
+bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum);
+bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int));
+bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 
+bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 
+bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 
+bool CreateDataDeltaFrobeniusOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 
+bool CreateDataFragment(class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int));
+bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int));
+
+void CreateEnergy(class MatrixContainer &Energy, int MatrixNumber);
+void CreateMaxFragmentOrder(class MatrixContainer &Matrix, class KeySetsContainer &KeySet, int BondOrder);
+void CreateMinFragmentOrder(class MatrixContainer &Matrix, class KeySetsContainer &KeySet, int BondOrder);
+void CreateMinimumForce(class MatrixContainer &Force, int MatrixNumber);
+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);
+
+void CreatePlotHeader(ofstream &output, const char *prefix, const int keycolumns, const char *key, const char *logscale, const char *extraline, const int mxtics, const int xtics, const char *xlabel, const char *ylabel);
+bool CreatePlotOrder(class MatrixContainer &Matrix, const class KeySetsContainer &KeySet, const char *dir, const char *prefix, const int keycolumns, const char *key, const char *logscale, const char *extraline, const int mxtics, const int xtics, const char *xlabel, const char *ylabel, const char *xrange, const char *yrange, const char *xargument, const char *uses, void (*CreatePlotLines)(ofstream &, class MatrixContainer &, const char *, const char *, const char *));
+
+void AbsEnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses);
+void EnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses);
+void ForceMagnitudePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses);
+void AbsFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses);
+void BoxesForcePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses);
+void BoxesFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses);
+
+#endif /*DATACREATOR_HPP_*/
Index: src/Fragmentation/joiner.cpp
===================================================================
--- src/Fragmentation/joiner.cpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
+++ src/Fragmentation/joiner.cpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
@@ -0,0 +1,269 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
+/** \file joiner.cpp
+ *
+ * Takes evaluated fragments (energy and forces) and by reading the factors files determines total energy
+ * and each force for the whole molecule.
+ *
+ */
+
+//============================ INCLUDES ===========================
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include <cstring>
+#include <cmath>
+
+#include "datacreator.hpp"
+#include "parser.hpp"
+#include "periodentafel.hpp"
+#include "CodePatterns/Log.hpp"
+#include "CodePatterns/Verbose.hpp"
+
+//============================== MAIN =============================
+
+int main(int argc, char **argv)
+{
+  periodentafel *periode = NULL; // and a period table of all elements
+  EnergyMatrix Energy;
+  EnergyMatrix EnergyFragments;
+  
+  EnergyMatrix Hcorrection;
+  EnergyMatrix HcorrectionFragments;
+  
+  ForceMatrix Force;
+  ForceMatrix ForceFragments;
+
+  HessianMatrix Hessian;
+  HessianMatrix HessianFragments;
+  
+  ForceMatrix Shielding;
+  ForceMatrix ShieldingPAS;
+  ForceMatrix ShieldingFragments;
+  ForceMatrix ShieldingPASFragments;
+  ForceMatrix Chi;
+  ForceMatrix ChiPAS;
+  ForceMatrix ChiFragments;
+  ForceMatrix ChiPASFragments;
+  KeySetsContainer KeySet;  
+  stringstream prefix;
+  char *dir = NULL;
+  bool NoHCorrection = false;
+  bool NoHessian = false;
+
+  DoLog(0) && (Log() << Verbose(0) << "Joiner" << endl);
+  DoLog(0) && (Log() << Verbose(0) << "======" << endl);
+
+  // Get the command line options
+  if (argc < 3) {
+    DoLog(0) && (Log() << Verbose(0) << "Usage: " << argv[0] << " <inputdir> <prefix> [elementsdb]" << endl);
+    DoLog(0) && (Log() << Verbose(0) << "<inputdir>\ttherein the output of a molecuilder fragmentation is expected, each fragment with a subdir containing an energy.all and a forces.all file." << endl);
+    DoLog(0) && (Log() << Verbose(0) << "<prefix>\tprefix of energy and forces file." << endl);
+    DoLog(0) && (Log() << Verbose(0) << "[elementsdb]\tpath to elements database, needed for shieldings." << endl);
+    return 1;
+  } else {
+    dir = new char[strlen(argv[2]) + 2];
+    strcpy(dir, "/");
+    strcat(dir, argv[2]);
+  }
+  if (argc > 3) {
+    periode = new periodentafel;
+    periode->LoadPeriodentafel(argv[3]);
+  }
+
+  // Test the given directory
+  if (!TestParams(argc, argv))
+    return 1;
+
+  // +++++++++++++++++ PARSING +++++++++++++++++++++++++++++++
+
+  // ------------- Parse through all Fragment subdirs --------
+  if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix, 0,0)) return 1;
+  if (!Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX, 0,0)) {
+    NoHCorrection = true;
+    DoLog(0) && (Log() << Verbose(0) << "No HCorrection matrices found, skipping these." << endl);
+  }
+  if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix, 0,0)) return 1;
+  if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix, 0,0)) {
+    NoHessian = true;
+    DoLog(0) && (Log() << Verbose(0) << "No hessian matrices found, skipping these." << endl);
+  }
+  if (periode != NULL) { // also look for PAS values
+    if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
+    if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
+    if (!Chi.ParseFragmentMatrix(argv[1], dir, ChiSuffix, 1, 0)) return 1;
+    if (!ChiPAS.ParseFragmentMatrix(argv[1], dir, ChiPASSuffix, 1, 0)) return 1;
+  }
+
+  // ---------- Parse the TE Factors into an array -----------------
+  if (!Energy.InitialiseIndices()) return 1;
+  if (!NoHCorrection) 
+    Hcorrection.InitialiseIndices();
+  
+  // ---------- Parse the Force indices into an array ---------------
+  if (!Force.ParseIndices(argv[1])) return 1;
+
+  // ---------- Parse the Hessian (=force) indices into an array ---------------
+  if (!NoHessian)
+    if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
+
+  // ---------- Parse the shielding indices into an array ---------------
+  if (periode != NULL) { // also look for PAS values
+    if(!Shielding.ParseIndices(argv[1])) return 1;
+    if(!ShieldingPAS.ParseIndices(argv[1])) return 1;
+    if(!Chi.ParseIndices(argv[1])) return 1;
+    if(!ChiPAS.ParseIndices(argv[1])) return 1;
+  }
+
+  // ---------- Parse the KeySets into an array ---------------
+  if (!KeySet.ParseKeySets(argv[1], Force.RowCounter, Force.MatrixCounter)) return 1;
+
+  if (!KeySet.ParseManyBodyTerms()) return 1;
+
+  if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return 1;
+  if (!NoHCorrection)  
+    HcorrectionFragments.AllocateMatrix(Hcorrection.Header, Hcorrection.MatrixCounter, Hcorrection.RowCounter, Hcorrection.ColumnCounter);
+  if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
+  if (!NoHessian)
+    if (!HessianFragments.AllocateMatrix(Hessian.Header, Hessian.MatrixCounter, Hessian.RowCounter, Hessian.ColumnCounter)) return 1;
+  if (periode != NULL) { // also look for PAS values
+    if (!ShieldingFragments.AllocateMatrix(Shielding.Header, Shielding.MatrixCounter, Shielding.RowCounter, Shielding.ColumnCounter)) return 1;
+    if (!ShieldingPASFragments.AllocateMatrix(ShieldingPAS.Header, ShieldingPAS.MatrixCounter, ShieldingPAS.RowCounter, ShieldingPAS.ColumnCounter)) return 1;
+    if (!ChiFragments.AllocateMatrix(Chi.Header, Chi.MatrixCounter, Chi.RowCounter, Chi.ColumnCounter)) return 1;
+    if (!ChiPASFragments.AllocateMatrix(ChiPAS.Header, ChiPAS.MatrixCounter, ChiPAS.RowCounter, ChiPAS.ColumnCounter)) return 1;
+  }
+
+  // ----------- Resetting last matrices (where full QM values are stored right now)
+  if(!Energy.SetLastMatrix(0., 0)) return 1;
+  if(!Force.SetLastMatrix(0., 2)) return 1;
+  if (!NoHessian)
+    if (!Hessian.SetLastMatrix(0., 0)) return 1;
+  if (periode != NULL) { // also look for PAS values
+    if(!Shielding.SetLastMatrix(0., 2)) return 1;
+    if(!ShieldingPAS.SetLastMatrix(0., 2)) return 1;
+    if(!Chi.SetLastMatrix(0., 2)) return 1;
+    if(!ChiPAS.SetLastMatrix(0., 2)) return 1;
+  }
+
+  // +++++++++++++++++ SUMMING +++++++++++++++++++++++++++++++
+
+  // --------- sum up and write for each order----------------
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    // --------- sum up energy --------------------
+    DoLog(0) && (Log() << Verbose(0) << "Summing energy of order " << BondOrder+1 << " ..." << endl);
+    if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return 1;
+    if (!NoHCorrection) { 
+      HcorrectionFragments.SumSubManyBodyTerms(Hcorrection, KeySet, BondOrder);
+      if (!Energy.SumSubEnergy(EnergyFragments, &HcorrectionFragments, KeySet, BondOrder, 1.)) return 1;
+      Hcorrection.SumSubEnergy(HcorrectionFragments, NULL, KeySet, BondOrder, 1.);
+    } else 
+      if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return 1;
+    // --------- sum up Forces --------------------
+    DoLog(0) && (Log() << Verbose(0) << "Summing forces of order " << BondOrder+1 << " ..." << endl);
+    if (!ForceFragments.SumSubManyBodyTerms(Force, KeySet, BondOrder)) return 1;
+    if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return 1;
+    // --------- sum up Hessian --------------------
+    if (!NoHessian) {
+      DoLog(0) && (Log() << Verbose(0) << "Summing Hessian of order " << BondOrder+1 << " ..." << endl);
+      if (!HessianFragments.SumSubManyBodyTerms(Hessian, KeySet, BondOrder)) return 1;
+      if (!Hessian.SumSubHessians(HessianFragments, KeySet, BondOrder, 1.)) return 1;
+    }
+    if (periode != NULL) { // also look for PAS values
+      DoLog(0) && (Log() << Verbose(0) << "Summing shieldings and susceptibilities of order " << BondOrder+1 << " ..." << endl);
+      if (!ShieldingFragments.SumSubManyBodyTerms(Shielding, 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, 1.)) return 1;
+      if (!ChiFragments.SumSubManyBodyTerms(Chi, KeySet, BondOrder)) return 1;
+      if (!Chi.SumSubForces(ChiFragments, KeySet, BondOrder, 1.)) return 1;
+      if (!ChiPASFragments.SumSubManyBodyTerms(ChiPAS, KeySet, BondOrder)) return 1;
+      if (!ChiPAS.SumSubForces(ChiPASFragments, KeySet, BondOrder, 1.)) return 1;
+    }
+
+    // --------- write the energy and forces file --------------------
+    prefix.str(" ");
+    prefix << dir << OrderSuffix << (BondOrder+1);
+    DoLog(0) && (Log() << Verbose(0) << "Writing files " << argv[1] << prefix.str() << ". ..." << endl);
+    // energy
+    if (!Energy.WriteLastMatrix(argv[1], (prefix.str()).c_str(), EnergySuffix)) return 1;
+    // forces
+    if (!Force.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ForcesSuffix)) return 1;
+    // hessian
+    if (!NoHessian)
+      if (!Hessian.WriteLastMatrix(argv[1], (prefix.str()).c_str(), HessianSuffix)) return 1;
+    // shieldings
+    if (periode != NULL) { // also look for PAS values
+      if (!Shielding.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ShieldingSuffix)) return 1;
+      if (!ShieldingPAS.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ShieldingPASSuffix)) return 1;
+      if (!Chi.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ChiSuffix)) return 1;
+      if (!ChiPAS.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ChiPASSuffix)) return 1;
+    }
+  }
+  // fragments
+  prefix.str(" ");
+  prefix << dir << EnergyFragmentSuffix;
+  if (!EnergyFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+  if (!NoHCorrection) {
+    prefix.str(" ");
+    prefix << dir << HcorrectionFragmentSuffix;
+    if (!HcorrectionFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+  }
+  prefix.str(" ");
+  prefix << dir << ForceFragmentSuffix;
+  if (!ForceFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+  {
+    std::string fileprefix(FRAGMENTPREFIX);
+    fileprefix += ENERGYPERFRAGMENT;
+    if (!CreateDataFragment(EnergyFragments, KeySet, argv[1], fileprefix.c_str(), "fragment energy versus the Fragment No", "today", CreateEnergy)) return 1;
+  }
+  if (!NoHessian) {
+    prefix.str(" ");
+    prefix << dir << HessianFragmentSuffix;
+    if (!HessianFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+  }
+  if (periode != NULL) { // also look for PAS values
+    prefix.str(" ");
+    prefix << dir << ShieldingFragmentSuffix;
+    if (!ShieldingFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+    prefix.str(" ");
+    prefix << dir << ShieldingPASFragmentSuffix;
+    if (!ShieldingPASFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+    prefix.str(" ");
+    prefix << dir << ChiFragmentSuffix;
+    if (!ChiFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+    prefix.str(" ");
+    prefix << dir << ChiPASFragmentSuffix;
+    if (!ChiPASFragments.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], dir, EnergyFragmentSuffix)) return 1;
+  if (!NoHCorrection) Hcorrection.WriteLastMatrix(argv[1], dir, HcorrectionFragmentSuffix);
+  if (!Force.WriteLastMatrix(argv[1], dir, ForceFragmentSuffix)) return 1;
+  if (!NoHessian)
+    if (!Hessian.WriteLastMatrix(argv[1], dir, HessianFragmentSuffix)) return 1;
+  if (periode != NULL) { // also look for PAS values
+    if (!Shielding.WriteLastMatrix(argv[1], dir, ShieldingFragmentSuffix)) return 1;
+    if (!ShieldingPAS.WriteLastMatrix(argv[1], dir, ShieldingPASFragmentSuffix)) return 1;
+    if (!Chi.WriteLastMatrix(argv[1], dir, ChiFragmentSuffix)) return 1;
+    if (!ChiPAS.WriteLastMatrix(argv[1], dir, ChiPASFragmentSuffix)) return 1;
+  }
+
+  // exit
+  delete(periode);
+  delete[](dir);
+  DoLog(0) && (Log() << Verbose(0) << "done." << endl);
+  return 0;
+};
+
+//============================ END ===========================
Index: src/Fragmentation/parser.cpp
===================================================================
--- src/Fragmentation/parser.cpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
+++ src/Fragmentation/parser.cpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
@@ -0,0 +1,1214 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
+/** \file parsing.cpp
+ *
+ * Declarations of various class functions for the parsing of value files.
+ *
+ */
+
+// ======================================= INCLUDES ==========================================
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include <cstring>
+#include <fstream>
+#include <iomanip>
+
+#include "parser.hpp"
+#include "CodePatterns/Log.hpp"
+#include "CodePatterns/Verbose.hpp"
+#include "Helpers/defs.hpp"
+#include "Helpers/helpers.hpp"
+
+using namespace std;
+
+// ======================================= FUNCTIONS ==========================================
+
+/** Test if given filename can be opened.
+ * \param filename name of file
+ * \param test true - no error message, false - print error
+ * \return given file exists
+ */
+bool FilePresent(const char *filename, bool test)
+{
+  ifstream input;
+
+  input.open(filename, ios::in);
+  if (input.fail()) {
+    if (!test)
+      DoLog(0) && (Log() << Verbose(0) << endl << "FilePresent: Unable to open " << filename << ", is the directory correct?" << endl);
+    return false;
+  }
+  input.close();
+  return true;
+};
+
+/** Test the given parameters.
+ * \param argc argument count
+ * \param **argv arguments array
+ * \return given inputdir is valid
+ */
+bool TestParams(int argc, char **argv)
+{
+  ifstream input;
+  stringstream line;
+
+  line << argv[1] << FRAGMENTPREFIX << KEYSETFILE;
+  return FilePresent(line.str().c_str(), false);
+};
+
+/** Returns a string with \a i prefixed with 0s to match order of total number of molecules in digits.
+ * \param FragmentNumber total number of fragments to determine necessary number of digits
+ * \param digits number to create with 0 prefixed
+ * \return allocated(!) char array with number in digits, ten base.
+ */
+inline char *FixedDigitNumber(const int FragmentNumber, const int digits)
+{
+  char *returnstring;
+  int number = FragmentNumber;
+  int order = 0;
+  while (number != 0) { // determine number of digits needed
+    number = (int)floor(((double)number / 10.));
+    order++;
+    //Log() << Verbose(0) << "Number is " << number << ", order is " << order << "." << endl;
+  }
+  // allocate string
+  returnstring = new char[order + 2];
+  // terminate  and fill string array from end backward
+  returnstring[order] = '\0';
+  number = digits;
+  for (int i=order;i--;) {
+    returnstring[i] = '0' + (char)(number % 10);
+    number = (int)floor(((double)number / 10.));
+  }
+  //Log() << Verbose(0) << returnstring << endl;
+  return returnstring;
+};
+
+// ======================================= CLASS MatrixContainer =============================
+
+/** Constructor of MatrixContainer class.
+ */
+MatrixContainer::MatrixContainer() :
+  Indices(NULL)
+{
+  Header = new char*[1];
+  Matrix = new double**[1]; // one more each for the total molecule
+  RowCounter = new int[1];
+  ColumnCounter = new int[1];
+  Header[0] = NULL;
+  Matrix[0] = NULL;
+  RowCounter[0] = -1;
+  MatrixCounter = 0;
+  ColumnCounter[0] = -1;
+};
+
+/** Destructor of MatrixContainer class.
+ */
+MatrixContainer::~MatrixContainer() {
+  if (Matrix != NULL) {
+    // free Matrix[MatrixCounter]
+    if ((ColumnCounter != NULL) && (RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
+      for(int j=RowCounter[MatrixCounter]+1;j--;)
+        delete[](Matrix[MatrixCounter][j]);
+    //if (MatrixCounter != 0)
+    delete[](Matrix[MatrixCounter]);
+    // free all matrices but ultimate one
+    for(int i=MatrixCounter;i--;) {
+      if ((ColumnCounter != NULL) && (RowCounter != NULL)) {
+          for(int j=RowCounter[i]+1;j--;)
+            delete[](Matrix[i][j]);
+          delete[](Matrix[i]);
+      }
+    }
+    delete[](Matrix);
+  }
+  // free indices
+  if (Indices != NULL)
+    for(int i=MatrixCounter+1;i--;) {
+      delete[](Indices[i]);
+    }
+  delete[](Indices);
+  
+  // free header and counters
+  if (Header != NULL)
+    for(int i=MatrixCounter+1;i--;)
+      delete[](Header[i]);
+  delete[](Header);
+  delete[](RowCounter);
+  delete[](ColumnCounter);
+};
+
+/** 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) 
+{
+  DoLog(0) && (Log() << Verbose(0) << "Initialising indices");
+  if (Matrix == NULL) {
+    DoLog(0) && (Log() << Verbose(0) << " with trivial mapping." << endl);
+    Indices = new int*[MatrixCounter + 1];
+    for(int i=MatrixCounter+1;i--;) {
+      Indices[i] = new int[RowCounter[i]];
+      for(int j=RowCounter[i];j--;)
+        Indices[i][j] = j;
+    }
+  } else {
+    DoLog(0) && (Log() << Verbose(0) << " from other MatrixContainer." << endl);
+    if (MatrixCounter != Matrix->MatrixCounter)
+      return false;
+    Indices = new int*[MatrixCounter + 1];
+    for(int i=MatrixCounter+1;i--;) {
+      if (RowCounter[i] != Matrix->RowCounter[i])
+        return false;
+      Indices[i] = new int[Matrix->RowCounter[i]];
+      for(int j=Matrix->RowCounter[i];j--;) {
+        Indices[i][j] = Matrix->Indices[i][j];
+        //Log() << Verbose(0) << Indices[i][j] << "\t";
+      }
+      //Log() << Verbose(0) << endl;
+    }
+  }
+  return true; 
+};
+
+/** Parsing a number of matrices.
+ *    -# open the matrix file
+ *    -# skip some lines (\a skiplines)
+ *    -# scan header lines for number of columns
+ *    -# scan lines for number of rows
+ *    -# allocate matrix
+ *    -# loop over found column and row counts and parse in each entry
+ * \param &input input stream
+ * \param skiplines number of inital lines to skip
+ * \param skiplines number of inital columns to skip
+ * \param MatrixNr index number in Matrix array to parse into
+ * \return parsing successful
+ */
+bool MatrixContainer::ParseMatrix(std::istream &input, int skiplines, int skipcolumns, int MatrixNr)
+{
+  stringstream line;
+  string token;
+  char filename[1023];
+
+  if (input.fail()) {
+    DoeLog(1) && (eLog()<< Verbose(1) << endl << "MatrixContainer::ParseMatrix: Unable to parse istream." << endl);
+    //performCriticalExit();
+    return false;
+  }
+
+  // parse header
+  if (Header[MatrixNr] != NULL)
+    delete[] Header[MatrixNr];
+  Header[MatrixNr] = new char[1024];
+  for (int m=skiplines+1;m--;)
+    input.getline(Header[MatrixNr], 1023);
+  
+  // scan header for number of columns
+  line.str(Header[MatrixNr]);
+  for(int k=skipcolumns;k--;)
+    line >> Header[MatrixNr];
+  Log() << Verbose(0) << line.str() << endl;
+  ColumnCounter[MatrixNr]=0;
+  while ( getline(line,token, '\t') ) {
+    if (token.length() > 0)
+      ColumnCounter[MatrixNr]++;
+  }
+  Log() << Verbose(0) << line.str() << endl;
+  Log() << Verbose(1) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
+  if (ColumnCounter[MatrixNr] == 0) {
+    DoeLog(0) && (eLog()<< Verbose(0) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from ostream." << endl);
+    performCriticalExit();
+  }
+  
+  // scan rest for number of rows/lines
+  RowCounter[MatrixNr]=-1;    // counts one line too much
+  while (!input.eof()) {
+    input.getline(filename, 1023);
+    Log() << Verbose(0) << "Comparing: " << strncmp(filename,"MeanForce",9) << endl;
+    RowCounter[MatrixNr]++; // then line was not last MeanForce
+    if (strncmp(filename,"MeanForce", 9) == 0) {
+      break;
+    }
+  }
+  Log() << Verbose(1) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from input stream." << endl;
+  if (RowCounter[MatrixNr] == 0) {
+    DoeLog(0) && (eLog()<< Verbose(0) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from input stream." << endl);
+    performCriticalExit();
+  }
+
+  // allocate matrix if it's not zero dimension in one direction
+  if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
+    if (Matrix[MatrixNr] != NULL)
+      delete[] Matrix[MatrixNr];
+    Matrix[MatrixNr] = new double*[RowCounter[MatrixNr] + 1];
+    for(int j=0;j<RowCounter[MatrixNr]+1;j++)
+      Matrix[MatrixNr][j] = 0;
+
+    // parse in each entry for this matrix
+    input.clear();
+    input.seekg(ios::beg);
+    for (int m=skiplines+1;m--;)
+      input.getline(Header[MatrixNr], 1023);    // skip header
+    line.str(Header[MatrixNr]);
+    Log() << Verbose(0) << "Header: " << line.str() << endl;
+    for(int k=skipcolumns;k--;)  // skip columns in header too
+      line >> filename;
+    strncpy(Header[MatrixNr], line.str().c_str(), 1023);  
+    for(int j=0;j<RowCounter[MatrixNr];j++) {
+      if (Matrix[MatrixNr][j] != NULL)
+        delete[] Matrix[MatrixNr][j];
+      Matrix[MatrixNr][j] = new double[ColumnCounter[MatrixNr]];
+      for(int k=0;k<ColumnCounter[MatrixNr];k++)
+        Matrix[MatrixNr][j][k] = 0;
+
+      input.getline(filename, 1023);
+      stringstream lines(filename);
+      Log() << Verbose(2) << "Matrix at level " << j << ":";// << filename << endl;
+      for(int k=skipcolumns;k--;)
+        lines >> filename;
+      for(int k=0;(k<ColumnCounter[MatrixNr]) && (!lines.eof());k++) {
+        lines >> Matrix[MatrixNr][j][k];
+        Log() << Verbose(1) << " " << setprecision(2) << Matrix[MatrixNr][j][k] << endl;
+      }
+      if (Matrix[MatrixNr][ RowCounter[MatrixNr] ] != NULL)
+        delete[] Matrix[MatrixNr][ RowCounter[MatrixNr] ];
+      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = new double[ColumnCounter[MatrixNr]];
+      for(int j=ColumnCounter[MatrixNr];j--;)
+        Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
+    }
+  } else {
+    DoeLog(1) && (eLog()<< Verbose(1) << "Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl);
+  }
+  return true;
+};
+
+/** Parsing a number of matrices.
+ * -# First, count the number of matrices by counting lines in KEYSETFILE
+ * -# Then,
+ *    -# construct the fragment number
+ *    -# open the matrix file
+ *    -# skip some lines (\a skiplines)
+ *    -# scan header lines for number of columns
+ *    -# scan lines for number of rows
+ *    -# allocate matrix
+ *    -# loop over found column and row counts and parse in each entry
+ * -# Finally, allocate one additional matrix (\a MatrixCounter) containing combined or temporary values
+ * \param *name directory with files
+ * \param *prefix prefix of each matrix file
+ * \param *suffix suffix of each matrix file
+ * \param skiplines number of inital lines to skip
+ * \param skiplines number of inital columns to skip
+ * \return parsing successful
+ */
+bool MatrixContainer::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
+{
+  char filename[1023];
+  ifstream input;
+  char *FragmentNumber = NULL;
+  stringstream file;
+  string token;
+
+  // count the number of matrices
+  MatrixCounter = -1; // we count one too much
+  file << name << FRAGMENTPREFIX << KEYSETFILE;
+  input.open(file.str().c_str(), ios::in);
+  if (input.bad()) {
+    DoLog(0) && (Log() << Verbose(0) << endl << "MatrixContainer::ParseFragmentMatrix: Unable to open " << file.str() << ", is the directory correct?" << endl);
+    return false;
+  }
+  while (!input.eof()) {
+    input.getline(filename, 1023);
+    stringstream zeile(filename);
+    MatrixCounter++;
+  }
+  input.close();
+  DoLog(0) && (Log() << Verbose(0) << "Determined " << MatrixCounter << " fragments." << endl);
+
+  DoLog(0) && (Log() << Verbose(0) << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl);
+  delete[](Header);
+  delete[](Matrix);
+  delete[](RowCounter);
+  delete[](ColumnCounter);
+  Header = new char*[MatrixCounter + 1]; // one more each for the total molecule
+  Matrix = new double**[MatrixCounter + 1]; // one more each for the total molecule
+  RowCounter = new int[MatrixCounter + 1];
+  ColumnCounter = new int[MatrixCounter + 1];
+  for(int i=MatrixCounter+1;i--;) {
+    Matrix[i] = NULL;
+    Header[i] = NULL;
+    RowCounter[i] = -1;
+    ColumnCounter[i] = -1;
+  }
+  for(int i=0; i < MatrixCounter;i++) {
+    // open matrix file
+    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
+    file.str(" ");
+    file << name << FRAGMENTPREFIX << FragmentNumber << prefix << suffix;
+    std::ifstream input(file.str().c_str());
+    DoLog(0) &&( Log() << Verbose(0) << "Opening " << file.str() << " ... " << endl);
+    if (!ParseMatrix(input, skiplines, skipcolumns, i)) {
+      input.close();
+      return false;
+    }
+    input.close();
+    delete[](FragmentNumber);
+  }
+  return true;
+};
+
+/** Allocates and resets the memory for a number \a MCounter of matrices.
+ * \param **GivenHeader Header line for each matrix
+ * \param MCounter number of matrices
+ * \param *RCounter number of rows for each matrix
+ * \param *CCounter number of columns for each matrix
+ * \return Allocation successful 
+ */
+bool MatrixContainer::AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter)
+{
+  MatrixCounter = MCounter;
+  Header = new char*[MatrixCounter + 1];
+  Matrix = new double**[MatrixCounter + 1]; // one more each for the total molecule
+  RowCounter = new int[MatrixCounter + 1];
+  ColumnCounter = new int[MatrixCounter + 1];
+  for(int i=MatrixCounter+1;i--;) {
+    Header[i] = new char[1024];
+    strncpy(Header[i], GivenHeader[i], 1023);
+    RowCounter[i] = RCounter[i];
+    ColumnCounter[i] = CCounter[i];
+    Matrix[i] = new double*[RowCounter[i] + 1];
+    for(int j=RowCounter[i]+1;j--;) {
+      Matrix[i][j] = new double[ColumnCounter[i]];
+      for(int k=ColumnCounter[i];k--;)
+        Matrix[i][j][k] = 0.;
+    }
+  }
+  return true;
+};
+
+/** Resets all values in MatrixContainer::Matrix.
+ * \return true if successful
+ */
+bool MatrixContainer::ResetMatrix()
+{
+  for(int i=MatrixCounter+1;i--;)
+    for(int j=RowCounter[i]+1;j--;)
+      for(int k=ColumnCounter[i];k--;)
+        Matrix[i][j][k] = 0.;
+   return true;
+};
+
+/** 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[i];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[i];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
+ * \param skipcolumns skip initial columns
+ * \return true if successful
+ */
+bool MatrixContainer::SetLastMatrix(double value, int skipcolumns)
+{
+  for(int j=RowCounter[MatrixCounter]+1;j--;)
+    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
+      Matrix[MatrixCounter][j][k] = value;
+   return true;
+};
+
+/** Sets all values in the last of MatrixContainer::Matrix to \a value.
+ * \param **values matrix with each value (must have at least same dimensions!)
+ * \param skipcolumns skip initial columns
+ * \return true if successful
+ */
+bool MatrixContainer::SetLastMatrix(double **values, int skipcolumns)
+{
+  for(int j=RowCounter[MatrixCounter]+1;j--;)
+    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
+      Matrix[MatrixCounter][j][k] = values[j][k];
+   return true;
+};
+
+/** Sums the entries with each factor and put into last element of \a ***Matrix.
+ * Sums over "E"-terms to create the "F"-terms
+ * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
+ * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
+ * \param Order bond order
+ * \return true if summing was successful
+ */
+bool MatrixContainer::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySets, int Order)
+{
+  // go through each order
+  for (int CurrentFragment=0;CurrentFragment<KeySets.FragmentsPerOrder[Order];CurrentFragment++) {
+    //Log() << Verbose(0) << "Current Fragment is " << CurrentFragment << "/" << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
+    // then go per order through each suborder and pick together all the terms that contain this fragment
+    for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
+      for (int j=0;j<KeySets.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
+        if (KeySets.Contains(KeySets.OrderSet[Order][CurrentFragment], KeySets.OrderSet[SubOrder][j])) {
+          //Log() << Verbose(0) << "Current other fragment is " << j << "/" << KeySets.OrderSet[SubOrder][j] << "." << endl;
+          // if the fragment's indices are all in the current fragment
+          for(int k=0;k<RowCounter[ KeySets.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
+            int m = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][k];
+            //Log() << Verbose(0) << "Current index is " << k << "/" << m << "." << endl;
+            if (m != -1) { // if it's not an added hydrogen
+              for (int l=0;l<RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
+                //Log() << Verbose(0) << "Comparing " << m << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
+                if (m == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l]) {
+                  m = l;
+                  break;
+                }
+              }
+              //Log() << Verbose(0) << "Corresponding index in CurrentFragment is " << m << "." << endl;
+              if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
+                DoeLog(0) && (eLog()<< Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current force index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl);
+                performCriticalExit();
+                return false;
+              }
+              if (Order == SubOrder) { // equal order is always copy from Energies
+                for(int l=ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l--;) // then adds/subtract each column
+                  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
+              } else {
+                for(int l=ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l--;)
+                  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
+              }
+            }
+            //if ((ColumnCounter[ KeySets.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
+             //Log() << Verbose(0) << "Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
+          }
+        } else {
+          //Log() << Verbose(0) << "Fragment " << KeySets.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
+        }
+      }
+    }
+   //Log() << Verbose(0) << "Final Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << KeySets.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][KeySets.AtomCounter[0]-1][1] << endl;
+  }
+
+  return true;
+};
+
+/** Writes the summed total fragment terms \f$F_{ij}\f$ to file.
+ * \param *name inputdir
+ * \param *prefix prefix before \a EnergySuffix
+ * \return file was written
+ */
+bool MatrixContainer::WriteTotalFragments(const char *name, const char *prefix)
+{
+  ofstream output;
+  char *FragmentNumber = NULL;
+
+  DoLog(0) && (Log() << Verbose(0) << "Writing fragment files." << endl);
+  for(int i=0;i<MatrixCounter;i++) {
+    stringstream line;
+    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
+    line << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix;
+    delete[](FragmentNumber);
+    output.open(line.str().c_str(), ios::out);
+    if (output == NULL) {
+      DoeLog(0) && (eLog()<< Verbose(0) << "MatrixContainer::WriteTotalFragments: Unable to open output energy file " << line.str() << "!" << endl);
+      performCriticalExit();
+      return false;
+    }
+    output << Header[i] << endl;
+    for(int j=0;j<RowCounter[i];j++) {
+      for(int k=0;k<ColumnCounter[i];k++)
+        output << scientific << Matrix[i][j][k] << "\t";
+      output << endl;
+    }
+    output.close();
+  }
+  return true;
+};
+
+/** Writes the summed total values in the last matrix to file.
+ * \param *name inputdir
+ * \param *prefix prefix
+ * \param *suffix suffix
+ * \return file was written
+ */
+bool MatrixContainer::WriteLastMatrix(const char *name, const char *prefix, const char *suffix)
+{
+  ofstream output;
+  stringstream line;
+
+  DoLog(0) && (Log() << Verbose(0) << "Writing matrix values of " << suffix << "." << endl);
+  line << name << prefix << suffix;
+  output.open(line.str().c_str(), ios::out);
+  if (output == NULL) {
+    DoeLog(0) && (eLog()<< Verbose(0) << "MatrixContainer::WriteLastMatrix: Unable to open output matrix file " << line.str() << "!" << endl);
+    performCriticalExit();
+    return false;
+  }
+  output << Header[MatrixCounter] << endl;
+  for(int j=0;j<RowCounter[MatrixCounter];j++) {
+    for(int k=0;k<ColumnCounter[MatrixCounter];k++)
+      output << scientific << Matrix[MatrixCounter][j][k] << "\t";
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+std::ostream & operator << (std::ostream &ost, const MatrixContainer &m)
+{
+  for (int i=0;i<=m.MatrixCounter;++i) {
+    ost << "Printing matrix " << i << ":" << std::endl;
+    for (int j=0;j<=m.RowCounter[i];++j) {
+      for (int k=0;k<m.ColumnCounter[i];++k) {
+        ost << m.Matrix[i][j][k] << " ";
+      }
+      ost << std::endl;
+    }
+  }
+  ost << std::endl;
+  return ost;
+}
+
+// ======================================= 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()
+{
+  DoLog(0) && (Log() << Verbose(0) << "Parsing energy indices." << endl);
+  Indices = new int*[MatrixCounter + 1];
+  for(int i=MatrixCounter+1;i--;) {
+    Indices[i] = new int[RowCounter[i]];
+    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 CorrectionFragments MatrixContainer with hydrogen saturation correction per fragments
+ * \param KeySets 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 EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySets, int Order, double sign)
+{
+  // sum energy
+  if (CorrectionFragments == NULL)
+    for(int i=KeySets.FragmentsPerOrder[Order];i--;)
+      for(int j=RowCounter[ KeySets.OrderSet[Order][i] ];j--;)
+        for(int k=ColumnCounter[ KeySets.OrderSet[Order][i] ];k--;)
+          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySets.OrderSet[Order][i] ][j][k];
+  else
+    for(int i=KeySets.FragmentsPerOrder[Order];i--;)
+      for(int j=RowCounter[ KeySets.OrderSet[Order][i] ];j--;)
+        for(int k=ColumnCounter[ KeySets.OrderSet[Order][i] ];k--;)
+          Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySets.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySets.OrderSet[Order][i] ][j][k]);
+  return true;
+};
+
+/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
+ * \param *name directory with files
+ * \param *prefix prefix of each matrix file
+ * \param *suffix suffix of each matrix file
+ * \param skiplines number of inital lines to skip
+ * \param skiplines number of inital columns to skip
+ * \return parsing successful
+ */
+bool EnergyMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
+{
+  char filename[1024];
+  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
+
+  if (status) {
+    // count maximum of columns
+    RowCounter[MatrixCounter] = 0;
+    ColumnCounter[MatrixCounter] = 0;
+    for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
+      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
+    DoLog(0) && (Log() << Verbose(0) << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl);
+    Matrix[MatrixCounter] = new double*[RowCounter[MatrixCounter] + 1];
+    for(int j=0;j<=RowCounter[MatrixCounter];j++)
+      Matrix[MatrixCounter][j] = new double[ColumnCounter[MatrixCounter]];
+    
+    // try independently to parse global energysuffix file if present
+    strncpy(filename, name, 1023);
+    strncat(filename, prefix, 1023-strlen(filename));
+    strncat(filename, suffix.c_str(), 1023-strlen(filename));
+    std::ifstream input(filename);
+    ParseMatrix(input, skiplines, skipcolumns, MatrixCounter);
+    input.close();
+  }
+  return status;
+};
+
+// ======================================= CLASS ForceMatrix =============================
+
+/** Parsing force Indices of each fragment
+ * \param *name directory with \a ForcesFile
+ * \return parsing successful
+ */
+bool ForceMatrix::ParseIndices(const char *name)
+{
+  ifstream input;
+  char *FragmentNumber = NULL;
+  char filename[1023];
+  stringstream line;
+
+  DoLog(0) && (Log() << Verbose(0) << "Parsing force indices for " << MatrixCounter << " matrices." << endl);
+  Indices = new int*[MatrixCounter + 1];
+  line << name << FRAGMENTPREFIX << FORCESFILE;
+  input.open(line.str().c_str(), ios::in);
+  //Log() << Verbose(0) << "Opening " << line.str() << " ... "  << input << endl;
+  if (input.fail()) {
+    DoLog(0) && (Log() << Verbose(0) << endl << "ForceMatrix::ParseIndices: Unable to open " << line.str() << ", is the directory correct?" << endl);
+    return false;
+  }
+  for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
+    // get the number of atoms for this fragment
+    input.getline(filename, 1023);
+    line.str(filename);
+    // parse the values
+    Indices[i] = new int[RowCounter[i]];
+    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
+    //Log() << Verbose(0) << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
+    delete[](FragmentNumber);
+    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
+      line >> Indices[i][j];
+      //Log() << Verbose(0) << " " << Indices[i][j];
+    }
+    //Log() << Verbose(0) << endl;
+  }
+  Indices[MatrixCounter] = new int[RowCounter[MatrixCounter]];
+  for(int j=RowCounter[MatrixCounter];j--;) {
+    Indices[MatrixCounter][j] = j;
+  }
+  input.close();
+  return true;
+};
+
+
+/** 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 KeySets 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 &KeySets, int Order, double sign)
+{
+  int FragmentNr;
+  // sum forces
+  for(int i=0;i<KeySets.FragmentsPerOrder[Order];i++) {
+    FragmentNr = KeySets.OrderSet[Order][i];
+    for(int l=0;l<RowCounter[ FragmentNr ];l++) {
+      int j = Indices[ FragmentNr ][l];
+      if (j > RowCounter[MatrixCounter]) {
+        DoeLog(0) && (eLog()<< Verbose(0) << "Current force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl);
+        performCriticalExit();
+        return false;
+      }
+      if (j != -1) {
+        //if (j == 0) Log() << Verbose(0) << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
+        for(int k=2;k<ColumnCounter[MatrixCounter];k++)
+          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
+      }
+    }
+  }
+  return true;
+};
+
+
+/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
+ * \param *name directory with files
+ * \param *prefix prefix of each matrix file
+ * \param *suffix suffix of each matrix file
+ * \param skiplines number of inital lines to skip
+ * \param skiplines number of inital columns to skip
+ * \return parsing successful
+ */
+bool ForceMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
+{
+  char filename[1023];
+  ifstream input;
+  stringstream file;
+  int nr;
+  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
+
+  if (status) {
+    // count number of atoms for last plus one matrix
+    file << name << FRAGMENTPREFIX << KEYSETFILE;
+    input.open(file.str().c_str(), ios::in);
+    if (input.fail()) {
+      DoLog(0) && (Log() << Verbose(0) << endl << "ForceMatrix::ParseFragmentMatrix: Unable to open " << file.str() << ", is the directory correct?" << endl);
+      return false;
+    }
+    RowCounter[MatrixCounter] = 0;
+    while (!input.eof()) {
+      input.getline(filename, 1023);
+      stringstream zeile(filename);
+      while (!zeile.eof()) {
+        zeile >> nr;
+        //Log() << Verbose(0) << "Current index: " << getNr() << "." << endl;
+        if (nr > RowCounter[MatrixCounter])
+          RowCounter[MatrixCounter] = nr;
+      }
+    }
+    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
+    DoLog(0) && (Log() << Verbose(0) << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl);
+    Matrix[MatrixCounter] = new double*[RowCounter[MatrixCounter] + 1];
+    for(int j=0;j<=RowCounter[MatrixCounter];j++)
+      Matrix[MatrixCounter][j] = new double[ColumnCounter[MatrixCounter]];
+
+    // try independently to parse global forcesuffix file if present
+    strncpy(filename, name, 1023);
+    strncat(filename, prefix, 1023-strlen(filename));
+    strncat(filename, suffix.c_str(), 1023-strlen(filename));
+    std::ifstream input(filename);
+    ParseMatrix(input, skiplines, skipcolumns, MatrixCounter);
+    input.close();
+  }
+ 
+
+  return status;
+};
+
+// ======================================= CLASS HessianMatrix =============================
+
+/** Parsing force Indices of each fragment
+ * \param *name directory with \a ForcesFile
+ * \return parsing successful 
+ */
+bool HessianMatrix::ParseIndices(char *name) 
+{
+  ifstream input;
+  char *FragmentNumber = NULL;
+  char filename[1023];
+  stringstream line;
+  
+  DoLog(0) && (Log() << Verbose(0) << "Parsing hessian indices for " << MatrixCounter << " matrices." << endl);
+  Indices = new int*[MatrixCounter + 1];
+  line << name << FRAGMENTPREFIX << FORCESFILE;
+  input.open(line.str().c_str(), ios::in);
+  //Log() << Verbose(0) << "Opening " << line.str() << " ... "  << input << endl;
+  if (input.fail()) {
+    DoLog(0) && (Log() << Verbose(0) << endl << "HessianMatrix::ParseIndices: Unable to open " << line.str() << ", is the directory correct?" << endl);
+    return false;
+  }
+  for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
+    // get the number of atoms for this fragment
+    input.getline(filename, 1023);
+    line.str(filename);
+    // parse the values
+    Indices[i] = new int[RowCounter[i]];
+    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
+    //Log() << Verbose(0) << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
+    delete[](FragmentNumber);
+    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
+      line >> Indices[i][j];
+      //Log() << Verbose(0) << " " << Indices[i][j];
+    }
+    //Log() << Verbose(0) << endl;
+  }
+  Indices[MatrixCounter] = new int[RowCounter[MatrixCounter]];
+  for(int j=RowCounter[MatrixCounter];j--;) {
+    Indices[MatrixCounter][j] = j;
+  }
+  input.close();
+  return true;
+};
+
+
+/** Sums the hessian entries and puts into last element of \a HessianMatrix::Matrix.
+ * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
+ * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
+ * \param Order bond order
+ *  \param sign +1 or -1
+ * \return true if summing was successful
+ */
+bool HessianMatrix::SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySets, int Order, double sign)
+{
+  int FragmentNr;
+  // sum forces
+  for(int i=0;i<KeySets.FragmentsPerOrder[Order];i++) {
+    FragmentNr = KeySets.OrderSet[Order][i];
+    for(int l=0;l<RowCounter[ FragmentNr ];l++) {
+      int j = Indices[ FragmentNr ][l];
+      if (j > RowCounter[MatrixCounter]) {
+        DoeLog(0) && (eLog()<< Verbose(0) << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl);
+        performCriticalExit();
+        return false;
+      }
+      if (j != -1) {
+        for(int m=0;m<ColumnCounter[ FragmentNr ];m++) {
+          int k = Indices[ FragmentNr ][m];
+          if (k > ColumnCounter[MatrixCounter]) {
+            DoeLog(0) && (eLog()<< Verbose(0) << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl);
+            performCriticalExit();
+            return false;
+          }
+          if (k != -1) {
+            //Log() << Verbose(0) << "Adding " << sign*Fragments.Matrix[ FragmentNr ][l][m] << " from [" << l << "][" << m << "] onto [" << j << "][" << k << "]." << endl;
+            Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][m];
+          }
+        }
+      }
+    }
+  }
+  return true;
+};
+
+/** Constructor for class HessianMatrix.
+ */
+HessianMatrix::HessianMatrix() :
+  MatrixContainer(),
+  IsSymmetric(true)
+{}
+
+/** Sums the hessian entries with each factor and put into last element of \a ***Matrix.
+ * Sums over "E"-terms to create the "F"-terms
+ * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
+ * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
+ * \param Order bond order
+ * \return true if summing was successful
+ */
+bool HessianMatrix::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySets, int Order)
+{
+  // go through each order
+  for (int CurrentFragment=0;CurrentFragment<KeySets.FragmentsPerOrder[Order];CurrentFragment++) {
+    //Log() << Verbose(0) << "Current Fragment is " << CurrentFragment << "/" << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
+    // then go per order through each suborder and pick together all the terms that contain this fragment
+    for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
+      for (int j=0;j<KeySets.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
+        if (KeySets.Contains(KeySets.OrderSet[Order][CurrentFragment], KeySets.OrderSet[SubOrder][j])) {
+          //Log() << Verbose(0) << "Current other fragment is " << j << "/" << KeySets.OrderSet[SubOrder][j] << "." << endl;
+          // if the fragment's indices are all in the current fragment
+          for(int k=0;k<RowCounter[ KeySets.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
+            int m = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][k];
+            //Log() << Verbose(0) << "Current row index is " << k << "/" << m << "." << endl;
+            if (m != -1) { // if it's not an added hydrogen
+              for (int l=0;l<RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
+                //Log() << Verbose(0) << "Comparing " << m << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
+                if (m == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l]) {
+                  m = l;
+                  break;  
+                }
+              }
+              //Log() << Verbose(0) << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl;
+              if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
+                DoeLog(0) && (eLog()<< Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl);
+                performCriticalExit();
+                return false;
+              }
+              
+              for(int l=0;l<ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l++) {
+                int n = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][l];
+                //Log() << Verbose(0) << "Current column index is " << l << "/" << n << "." << endl;
+                if (n != -1) { // if it's not an added hydrogen
+                  for (int p=0;p<ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment
+                    //Log() << Verbose(0) << "Comparing " << n << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][p] << "." << endl;
+                    if (n == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][p]) {
+                      n = p;
+                      break;  
+                    }
+                  }
+                  //Log() << Verbose(0) << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl;
+                  if (n > ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
+                    DoeLog(0) && (eLog()<< Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl);
+                    performCriticalExit();
+                    return false;
+                  }
+                  if (Order == SubOrder) { // equal order is always copy from Energies
+                    //Log() << Verbose(0) << "Adding " << MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
+                    Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
+                  } else {
+                    //Log() << Verbose(0) << "Subtracting " << Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
+                    Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
+                  }
+                }
+              }
+            }
+            //if ((ColumnCounter[ KeySets.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
+             //Log() << Verbose(0) << "Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
+          }
+        } else {
+          //Log() << Verbose(0) << "Fragment " << KeySets.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
+        }
+      }
+    }
+   //Log() << Verbose(0) << "Final Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << KeySets.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][KeySets.AtomCounter[0]-1][1] << endl;
+  }
+  
+  return true;
+};
+
+/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
+ * \param *name directory with files
+ * \param *prefix prefix of each matrix file
+ * \param *suffix suffix of each matrix file
+ * \param skiplines number of inital lines to skip
+ * \param skiplines number of inital columns to skip
+ * \return parsing successful
+ */ 
+bool HessianMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
+{
+  char filename[1023];
+  ifstream input;
+  stringstream file;
+  int nr;
+  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
+
+  if (status) {
+    // count number of atoms for last plus one matrix
+    file << name << FRAGMENTPREFIX << KEYSETFILE;
+    input.open(file.str().c_str(), ios::in);
+    if (input.fail()) {
+      DoLog(0) && (Log() << Verbose(0) << endl << "HessianMatrix::ParseFragmentMatrix: Unable to open " << file.str() << ", is the directory correct?" << endl);
+      return false;
+    }
+    RowCounter[MatrixCounter] = 0;
+    ColumnCounter[MatrixCounter] = 0;
+    while (!input.eof()) {
+      input.getline(filename, 1023);
+      stringstream zeile(filename);
+      while (!zeile.eof()) {
+        zeile >> nr;
+        //Log() << Verbose(0) << "Current index: " << getNr() << "." << endl;
+        if (nr > RowCounter[MatrixCounter]) {
+          RowCounter[MatrixCounter] = nr;
+          ColumnCounter[MatrixCounter] = nr;
+        }
+      }
+    }
+    RowCounter[MatrixCounter]++;    // Nr start at 0, count starts at 1
+    ColumnCounter[MatrixCounter]++;    // Nr start at 0, count starts at 1
+    input.close();
+  
+    // allocate last plus one matrix
+    DoLog(0) && (Log() << Verbose(0) << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl);
+    Matrix[MatrixCounter] = new double*[RowCounter[MatrixCounter] + 1];
+    for(int j=0;j<=RowCounter[MatrixCounter];j++)
+      Matrix[MatrixCounter][j] = new double[ColumnCounter[MatrixCounter]];
+
+    // try independently to parse global forcesuffix file if present
+    strncpy(filename, name, 1023);
+    strncat(filename, prefix, 1023-strlen(filename));
+    strncat(filename, suffix.c_str(), 1023-strlen(filename));
+    std::ifstream input(filename);
+    ParseMatrix(input, skiplines, skipcolumns, MatrixCounter);
+    input.close();
+  }
+
+
+  return status;
+};
+
+// ======================================= CLASS KeySetsContainer =============================
+
+/** Constructor of KeySetsContainer class.
+ */
+KeySetsContainer::KeySetsContainer() :
+  KeySets(NULL),
+  AtomCounter(NULL),
+  FragmentCounter(0),
+  Order(0),
+  FragmentsPerOrder(0),
+  OrderSet(NULL)
+{};
+
+/** Destructor of KeySetsContainer class.
+ */
+KeySetsContainer::~KeySetsContainer() {
+  for(int i=FragmentCounter;i--;)
+    delete[](KeySets[i]);
+  for(int i=Order;i--;)
+    delete[](OrderSet[i]);
+  delete[](KeySets);
+  delete[](OrderSet);
+  delete[](AtomCounter);
+  delete[](FragmentsPerOrder);
+};
+
+/** Parsing KeySets into array.
+ * \param *name directory with keyset file
+ * \param *ACounter number of atoms per fragment
+ * \param FCounter number of fragments
+ * \return parsing succesful
+ */
+bool KeySetsContainer::ParseKeySets(const char *name, const int *ACounter, const int FCounter) {
+  ifstream input;
+  char *FragmentNumber = NULL;
+  stringstream file;
+  char filename[1023];
+
+  FragmentCounter = FCounter;
+  DoLog(0) && (Log() << Verbose(0) << "Parsing key sets." << endl);
+  KeySets = new int*[FragmentCounter];
+  for(int i=FragmentCounter;i--;)
+    KeySets[i] = NULL;
+  file << name << FRAGMENTPREFIX << KEYSETFILE;
+  input.open(file.str().c_str(), ios::in);
+  if (input.fail()) {
+    DoLog(0) && (Log() << Verbose(0) << endl << "KeySetsContainer::ParseKeySets: Unable to open " << file.str() << ", is the directory correct?" << endl);
+    return false;
+  }
+
+  AtomCounter = new int[FragmentCounter];
+  for(int i=0;(i<FragmentCounter) && (!input.eof());i++) {
+    stringstream line;
+    AtomCounter[i] = ACounter[i];
+    // parse the values
+    KeySets[i] = new int[AtomCounter[i]];
+    for(int j=AtomCounter[i];j--;)
+      KeySets[i][j] = -1;
+    FragmentNumber = FixedDigitNumber(FragmentCounter, i);
+    //Log() << Verbose(0) << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
+    delete[](FragmentNumber);
+    input.getline(filename, 1023);
+    line.str(filename);
+    for(int j=0;(j<AtomCounter[i]) && (!line.eof());j++) {
+      line >> KeySets[i][j];
+      //Log() << Verbose(0) << " " << KeySets[i][j];
+    }
+    //Log() << Verbose(0) << endl;
+  }
+  input.close();
+  return true;
+};
+
+/** Parse many body terms, associating each fragment to a certain bond order.
+ * \return parsing succesful
+ */
+bool KeySetsContainer::ParseManyBodyTerms()
+{
+  int Counter;
+
+  DoLog(0) && (Log() << Verbose(0) << "Creating Fragment terms." << endl);
+  // scan through all to determine maximum order
+  Order=0;
+  for(int i=FragmentCounter;i--;) {
+    Counter=0;
+    for(int j=AtomCounter[i];j--;)
+      if (KeySets[i][j] != -1)
+        Counter++;
+    if (Counter > Order)
+      Order = Counter;
+  }
+  DoLog(0) && (Log() << Verbose(0) << "Found Order is " << Order << "." << endl);
+
+  // scan through all to determine fragments per order
+  FragmentsPerOrder = new int[Order];
+  for(int i=Order;i--;)
+    FragmentsPerOrder[i] = 0;
+  for(int i=FragmentCounter;i--;) {
+    Counter=0;
+    for(int j=AtomCounter[i];j--;)
+      if (KeySets[i][j] != -1)
+        Counter++;
+    FragmentsPerOrder[Counter-1]++;
+  }
+  for(int i=0;i<Order;i++)
+    DoLog(0) && (Log() << Verbose(0) << "Found No. of Fragments of Order " << i+1 << " is " << FragmentsPerOrder[i] << "." << endl);
+
+  // scan through all to gather indices to each order set
+  OrderSet = new int*[Order];
+  for(int i=Order;i--;)
+    OrderSet[i] = new int[FragmentsPerOrder[i]];
+  for(int i=Order;i--;)
+    FragmentsPerOrder[i] = 0;
+  for(int i=FragmentCounter;i--;) {
+    Counter=0;
+    for(int j=AtomCounter[i];j--;)
+      if (KeySets[i][j] != -1)
+        Counter++;
+    OrderSet[Counter-1][FragmentsPerOrder[Counter-1]] = i;
+    FragmentsPerOrder[Counter-1]++;
+  }
+  DoLog(0) && (Log() << Verbose(0) << "Printing OrderSet." << endl);
+  for(int i=0;i<Order;i++) {
+    for (int j=0;j<FragmentsPerOrder[i];j++) {
+      DoLog(0) && (Log() << Verbose(0) << " " << OrderSet[i][j]);
+    }
+    DoLog(0) && (Log() << Verbose(0) << endl);
+  }
+  DoLog(0) && (Log() << Verbose(0) << endl);
+
+
+  return true;
+};
+
+/** Compares each entry in \a *SmallerSet if it is containted in \a *GreaterSet.
+ * \param GreaterSet index to greater set
+ * \param SmallerSet index to smaller set
+ * \return true if all keys of SmallerSet contained in GreaterSet
+ */
+bool KeySetsContainer::Contains(const int GreaterSet, const int SmallerSet)
+{
+  bool result = true;
+  bool intermediate;
+  if ((GreaterSet < 0) || (SmallerSet < 0) || (GreaterSet > FragmentCounter) || (SmallerSet > FragmentCounter)) // index out of bounds
+    return false;
+  for(int i=AtomCounter[SmallerSet];i--;) {
+    intermediate = false;
+    for (int j=AtomCounter[GreaterSet];j--;)
+      intermediate = (intermediate || ((KeySets[SmallerSet][i] == KeySets[GreaterSet][j]) || (KeySets[SmallerSet][i] == -1)));
+    result = result && intermediate;
+  }
+
+  return result;
+};
+
+
+// ======================================= END =============================================
Index: src/Fragmentation/parser.hpp
===================================================================
--- src/Fragmentation/parser.hpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
+++ src/Fragmentation/parser.hpp	(revision bc28b0f3f8149d6034a044085ee96722cb9a75e3)
@@ -0,0 +1,140 @@
+/** \file parser.hpp
+ *
+ * Definitions of various class functions for the parsing of value files.
+ *
+ */
+
+
+#ifndef PARSER_HPP_
+#define PARSER_HPP_
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+/****************************************** forward declarations *****************************/
+
+class EnergyMatrix;
+class ForceMatrix;
+class HessianMatrix;
+class KeySetsContainer;
+class MatrixContainer;
+
+/********************************************** definitions *********************************/
+
+#define EnergySuffix ".energy.all"
+#define EnergyFragmentSuffix ".energyfragment.all"
+#define ForcesSuffix ".forces.all"
+#define ForceFragmentSuffix ".forcefragment.all"
+#define HcorrectionSuffix ".Hcorrection.all"
+#define HcorrectionFragmentSuffix ".Hcorrectionfragment.all"
+#define HessianSuffix ".hessian_xx.all"
+#define HessianFragmentSuffix ".hessianfragment_xx.all"
+#define OrderSuffix ".Order"
+#define ShieldingSuffix ".sigma_all.csv"
+#define ShieldingPASSuffix ".sigma_all_PAS.csv"
+#define ShieldingFragmentSuffix ".sigma_all_fragment.all"
+#define ShieldingPASFragmentSuffix ".sigma_all_PAS_fragment.all"
+#define ChiSuffix ".chi_all.csv"
+#define ChiPASSuffix ".chi_all_PAS.csv"
+#define ChiFragmentSuffix ".chi_all_fragment.all"
+#define ChiPASFragmentSuffix ".chi_all_PAS_fragment.all"
+#define TimeSuffix ".speed"
+
+// ======================================= FUNCTIONS ==========================================
+
+bool FilePresent(const char *filename, bool test);
+bool TestParams(int argc, char **argv);
+
+// ======================================= CLASS KeySetsContainer =============================
+
+class KeySetsContainer {
+  public:
+    int **KeySets;
+    int *AtomCounter;
+    int FragmentCounter;
+    int Order;
+    int *FragmentsPerOrder;
+    int **OrderSet;
+
+  KeySetsContainer();
+  ~KeySetsContainer();
+
+  bool ParseKeySets(const char *name, const int *ACounter, const int FCounter);
+  bool ParseManyBodyTerms();
+  bool Contains(const int GreaterSet, const int SmallerSet);
+};
+
+// ======================================= CLASS MatrixContainer =============================
+
+class MatrixContainer {
+  public:
+    double ***Matrix;
+    int **Indices;
+    char **Header;
+    int MatrixCounter;
+    int *RowCounter;
+    int *ColumnCounter;
+  
+  MatrixContainer();
+  virtual ~MatrixContainer();
+  
+  bool InitialiseIndices(class MatrixContainer *Matrix = NULL);
+  bool ParseMatrix(std::istream &input, int skiplines, int skipcolumns, int MatrixNr);
+  virtual bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns);
+  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);
+  //bool ParseIndices();
+  //bool SumSubValues();
+  bool SumSubManyBodyTerms(class MatrixContainer &Matrix, class KeySetsContainer &KeySets, int Order);
+  bool WriteTotalFragments(const char *name, const char *prefix);
+  bool WriteLastMatrix(const char *name, const char *prefix, const char *suffix);
+};
+
+std::ostream & operator << (std::ostream &ost, const MatrixContainer &m);
+
+// ======================================= CLASS EnergyMatrix =============================
+
+class EnergyMatrix : public MatrixContainer {
+  public:
+    bool ParseIndices();
+    bool SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySets, int Order, double sign);
+    bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns);
+};
+
+// ======================================= CLASS ForceMatrix =============================
+
+class ForceMatrix : public MatrixContainer {
+  public:
+    bool ParseIndices(const char *name);
+    bool SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySets, int Order, double sign);
+    bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns);
+};
+
+// ======================================= CLASS HessianMatrix =============================
+
+class HessianMatrix : public MatrixContainer { 
+  public: 
+    HessianMatrix();
+    //~HessianMatrix();
+    bool ParseIndices(char *name);
+    bool SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySets, int Order);
+    bool SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySets, int Order, double sign);
+    bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns);
+  private:
+    bool IsSymmetric;
+};
+
+
+// ======================================= END =============================================
+
+#endif /*PARSER_HPP_*/
