Index: src/Fragmentation/analyzer.cpp
===================================================================
--- src/Fragmentation/analyzer.cpp	(revision 1d17e57c9c64b6781a9d183fcf18872bae9cf9b9)
+++ 	(revision )
@@ -1,603 +1,0 @@
-/*
- * Project: MoleCuilder
- * Description: creates and alters molecular systems
- * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
- * 
- *
- *   This file is part of MoleCuilder.
- *
- *    MoleCuilder is free software: you can redistribute it and/or modify
- *    it under the terms of the GNU General Public License as published by
- *    the Free Software Foundation, either version 2 of the License, or
- *    (at your option) any later version.
- *
- *    MoleCuilder is distributed in the hope that it will be useful,
- *    but WITHOUT ANY WARRANTY; without even the implied warranty of
- *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *    GNU General Public License for more details.
- *
- *    You should have received a copy of the GNU General Public License
- *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-/** \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 "Element/periodentafel.hpp"
-#include "Fragmentation/defs.hpp"
-#include "Fragmentation/EnergyMatrix.hpp"
-#include "Fragmentation/ForceMatrix.hpp"
-#include "Fragmentation/helpers.hpp"
-#include "Fragmentation/HessianMatrix.hpp"
-#include "Fragmentation/KeySetsContainer.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;
-
-  LOG(0, "ANOVA Analyzer");
-  LOG(0, "==============");
-
-  // Get the command line options
-  if (argc < 4) {
-    LOG(0, "Usage: " << argv[0] << " <inputdir> <prefix> <outputdir> [elementsdb]");
-    LOG(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.");
-    LOG(0, "<prefix>\tprefix of energy and forces file.");
-    LOG(0, "<outputdir>\tcreated plotfiles and datafiles are placed into this directory ");
-    LOG(0, "[elementsdb]\tpath to elements database, needed for shieldings.");
-    return 1;
-  } else {
-    dir = new char[strlen(argv[2]) + 2];
-    strcpy(dir, "/");
-    strcat(dir, argv[2]);
-  }
-
-  if (argc > 4) {
-    LOG(0, "Loading periodentafel.");
-    periode = new periodentafel;
-    periode->LoadPeriodentafel(argv[4]);
-  }
-
-  // Test the given directory
-  if (!TestParams(argc, argv)) {
-    delete[] dir;
-    delete periode;
-    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;
-    ELOG(2, "No HCorrection file found, skipping these.");
-  }
-  
-  if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
-  if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix,0,0)) {
-    NoHessian = true;
-    ELOG(2, "No Hessian file found, skipping these.");
-  }
-  if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) {
-    NoTime = true;
-    ELOG(2, "No speed file found, skipping these.");
-  }
-  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 ---------------
-  {
-    std::stringstream filename;
-    filename << FRAGMENTPREFIX << KEYSETFILE;
-    if (!KeySet.ParseKeySets(argv[1], filename.str(), 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 ++++++++++++++++++++++++++++++
-
-  LOG(0, "Analyzing ...");
-
-  // ======================================= 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);
-  LOG(0, "done.");
-  return 0;
-};
-
-//============================ END ===========================
Index: src/Fragmentation/datacreator.cpp
===================================================================
--- src/Fragmentation/datacreator.cpp	(revision 1d17e57c9c64b6781a9d183fcf18872bae9cf9b9)
+++ 	(revision )
@@ -1,866 +1,0 @@
-/*
- * Project: MoleCuilder
- * Description: creates and alters molecular systems
- * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
- * 
- *
- *   This file is part of MoleCuilder.
- *
- *    MoleCuilder is free software: you can redistribute it and/or modify
- *    it under the terms of the GNU General Public License as published by
- *    the Free Software Foundation, either version 2 of the License, or
- *    (at your option) any later version.
- *
- *    MoleCuilder is distributed in the hope that it will be useful,
- *    but WITHOUT ANY WARRANTY; without even the implied warranty of
- *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *    GNU General Public License for more details.
- *
- *    You should have received a copy of the GNU General Public License
- *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-/** \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 "Fragmentation/EnergyMatrix.hpp"
-#include "Fragmentation/ForceMatrix.hpp"
-#include "Fragmentation/HessianMatrix.hpp"
-#include "Fragmentation/KeySetsContainer.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) {
-    LOG(0, "Unable to open " << name.str() << " for writing, is directory correct?");
-    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) {
-    LOG(0, "Unable to open " << name.str() << " for writing, is directory correct?");
-    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;
-  LOG(0, msg);
-  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;
-  LOG(0, msg);
-  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;
-  LOG(0, msg);
-  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;
-  LOG(0, msg);
-  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;
-  LOG(0, msg);
-  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(0, "Current order is " << BondOrder << ".");
-    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;
-  LOG(0, msg);
-  output << "# " << msg << ", created on " << datum;
-  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
-  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
-    //LOG(0, "Current order is " << BondOrder << ".");
-    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; 
-  LOG(0, msg);
-  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(0, "Current order is " << BondOrder << ".");
-    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; 
-  LOG(0, msg);
-  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(0, "Current order is " << BondOrder << ".");
-    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; 
-  LOG(0, msg);
-  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(0, "Current order is " << BondOrder << ".");
-    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;
-  LOG(0, msg);
-  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;
-  LOG(0, msg);
-  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 1d17e57c9c64b6781a9d183fcf18872bae9cf9b9)
+++ 	(revision )
@@ -1,64 +1,0 @@
-/** \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 1d17e57c9c64b6781a9d183fcf18872bae9cf9b9)
+++ 	(revision )
@@ -1,293 +1,0 @@
-/*
- * Project: MoleCuilder
- * Description: creates and alters molecular systems
- * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
- * 
- *
- *   This file is part of MoleCuilder.
- *
- *    MoleCuilder is free software: you can redistribute it and/or modify
- *    it under the terms of the GNU General Public License as published by
- *    the Free Software Foundation, either version 2 of the License, or
- *    (at your option) any later version.
- *
- *    MoleCuilder is distributed in the hope that it will be useful,
- *    but WITHOUT ANY WARRANTY; without even the implied warranty of
- *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *    GNU General Public License for more details.
- *
- *    You should have received a copy of the GNU General Public License
- *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-/** \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 "Element/periodentafel.hpp"
-#include "Fragmentation/defs.hpp"
-#include "Fragmentation/EnergyMatrix.hpp"
-#include "Fragmentation/ForceMatrix.hpp"
-#include "Fragmentation/helpers.hpp"
-#include "Fragmentation/HessianMatrix.hpp"
-#include "Fragmentation/KeySetsContainer.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;
-
-  LOG(0, "Joiner");
-  LOG(0, "======");
-
-  // Get the command line options
-  if (argc < 3) {
-    LOG(0, "Usage: " << argv[0] << " <inputdir> <prefix> [elementsdb]");
-    LOG(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.");
-    LOG(0, "<prefix>\tprefix of energy and forces file.");
-    LOG(0, "[elementsdb]\tpath to elements database, needed for shieldings.");
-    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;
-    LOG(0, "No HCorrection matrices found, skipping these.");
-  }
-  if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix, 0,0)) return 1;
-  if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix, 0,0)) {
-    NoHessian = true;
-    LOG(0, "No hessian matrices found, skipping these.");
-  }
-  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 ---------------
-  {
-    std::stringstream filename;
-    filename << FRAGMENTPREFIX << KEYSETFILE;
-    if (!KeySet.ParseKeySets(argv[1], filename.str(), 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 --------------------
-    LOG(0, "Summing energy of order " << BondOrder+1 << " ...");
-    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 --------------------
-    LOG(0, "Summing forces of order " << BondOrder+1 << " ...");
-    if (!ForceFragments.SumSubManyBodyTerms(Force, KeySet, BondOrder)) return 1;
-    if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return 1;
-    // --------- sum up Hessian --------------------
-    if (!NoHessian) {
-      LOG(0, "Summing Hessian of order " << BondOrder+1 << " ...");
-      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
-      LOG(0, "Summing shieldings and susceptibilities of order " << BondOrder+1 << " ...");
-      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);
-    LOG(0, "Writing files " << argv[1] << prefix.str() << ". ...");
-    // 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);
-  LOG(0, "done.");
-  return 0;
-};
-
-//============================ END ===========================
