/** \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 "datacreator.hpp"
#include "helpers.hpp"
#include "memoryallocator.hpp"
#include "parser.hpp"
#include "periodentafel.hpp"

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif


//============================== 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;
  ofstream output;
  ofstream output2;
  ofstream output3;
  ofstream output4;
  ifstream input;
  stringstream filename;
  time_t t = time(NULL);
  struct tm *ts = localtime(&t);
  char *datum = asctime(ts);
  stringstream Orderxrange;
  stringstream Fragmentxrange;
  stringstream yrange;
  char *dir = NULL;
  bool NoHessian = false;
  bool NoTime = false;
  bool NoHCorrection = true;
  int counter;

  cout << "ANOVA Analyzer" << endl;
  cout << "==============" << endl;

  // Get the command line options
  if (argc < 4) {
    cout << "Usage: " << argv[0] << " <inputdir> <prefix> <outputdir> [elementsdb]" << endl;
    cout << "<inputdir>\ttherein the output of a molecuilder fragmentation is expected, each fragment with a subdir containing an energy.all and a forces.all file." << endl;
    cout << "<prefix>\tprefix of energy and forces file." << endl;
    cout << "<outputdir>\tcreated plotfiles and datafiles are placed into this directory " << endl;
    cout << "[elementsdb]\tpath to elements database, needed for shieldings." << endl;
    return 1;
  } else {
    dir = Malloc<char>(strlen(argv[2]) + 2, "main: *dir");
    strcpy(dir, "/");
    strcat(dir, argv[2]);
  }

  if (argc > 4) {
    cout << "Loading periodentafel." << endl;
    periode = Malloc<periodentafel>(1, "main - periode");
    periode->LoadPeriodentafel(argv[4]);
  }

  // Test the given directory
  if (!TestParams(argc, argv))
    return 1;

  // +++++++++++++++++ PARSING +++++++++++++++++++++++++++++++

  // ------------- Parse through all Fragment subdirs --------
  if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
  if (!Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0)) {
    NoHCorrection = true;
    cout << "No HCorrection file found, skipping these." << endl;
  }
  
  if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
  if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix,0,0)) {
    NoHessian = true;
    cout << "No Hessian file found, skipping these." << endl;
  }
  if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) {
    NoTime = true;
    cout << "No speed file found, skipping these." << endl;
  }
  if (periode != NULL) { // also look for PAS values
    if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
    if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
    if (!Chi.ParseFragmentMatrix(argv[1], dir, ChiSuffix, 1, 0)) return 1;
    if (!ChiPAS.ParseFragmentMatrix(argv[1], dir, ChiPASSuffix, 1, 0)) return 1;
  }

  // ---------- Parse the TE Factors into an array -----------------
  if (!Energy.ParseIndices()) return 1;
  if (!NoHCorrection) Hcorrection.ParseIndices();

  // ---------- Parse the Force indices into an array ---------------
  if (!Force.ParseIndices(argv[1])) return 1;
  if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
  if (!ForceFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1;

  // ---------- Parse hessian indices into an array -----------------
  if (!NoHessian) {
    if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
    if (!HessianFragments.AllocateMatrix(Hessian.Header, Hessian.MatrixCounter, Hessian.RowCounter, Hessian.ColumnCounter)) return 1;
    if (!HessianFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
  }

  // ---------- Parse the shielding indices into an array ---------------
  if (periode != NULL) { // also look for PAS values
    if(!Shielding.ParseIndices(argv[1])) return 1;
    if(!ShieldingPAS.ParseIndices(argv[1])) return 1;
    if (!ShieldingFragments.AllocateMatrix(Shielding.Header, Shielding.MatrixCounter, Shielding.RowCounter, Shielding.ColumnCounter)) return 1;
    if (!ShieldingPASFragments.AllocateMatrix(ShieldingPAS.Header, ShieldingPAS.MatrixCounter, ShieldingPAS.RowCounter, ShieldingPAS.ColumnCounter)) return 1;
    if(!ShieldingFragments.ParseIndices(argv[1])) return 1;
    if(!ShieldingPASFragments.ParseIndices(argv[1])) return 1;
    if(!Chi.ParseIndices(argv[1])) return 1;
    if(!ChiPAS.ParseIndices(argv[1])) return 1;
    if (!ChiFragments.AllocateMatrix(Chi.Header, Chi.MatrixCounter, Chi.RowCounter, Chi.ColumnCounter)) return 1;
    if (!ChiPASFragments.AllocateMatrix(ChiPAS.Header, ChiPAS.MatrixCounter, ChiPAS.RowCounter, ChiPAS.ColumnCounter)) return 1;
    if(!ChiFragments.ParseIndices(argv[1])) return 1;
    if(!ChiPASFragments.ParseIndices(argv[1])) return 1;
  }

  // ---------- Parse the KeySets into an array ---------------
  if (!KeySet.ParseKeySets(argv[1], Force.RowCounter, Force.MatrixCounter)) return 1;
  if (!KeySet.ParseManyBodyTerms()) return 1;

  // ---------- Parse fragment files created by 'joiner' into an array -------------
  if (!EnergyFragments.ParseFragmentMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
  if (!NoHCorrection) 
    HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
  if (!ForceFragments.ParseFragmentMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
  if (!NoHessian)
    if (!HessianFragments.ParseFragmentMatrix(argv[1], dir, HessianFragmentSuffix,0,0)) return 1;
  if (periode != NULL) { // also look for PAS values
    if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingFragmentSuffix, 1, 0)) return 1;
    if (!ShieldingPASFragments.ParseFragmentMatrix(argv[1], dir, ShieldingPASFragmentSuffix, 1, 0)) return 1;
    if (!ChiFragments.ParseFragmentMatrix(argv[1], dir, ChiFragmentSuffix, 1, 0)) return 1;
    if (!ChiPASFragments.ParseFragmentMatrix(argv[1], dir, ChiPASFragmentSuffix, 1, 0)) return 1;
  }

  // +++++++++++++++ TESTING ++++++++++++++++++++++++++++++

  // print energy and forces to file
  filename.str("");
  filename << argv[3] << "/" << "energy-forces.all";
  output.open(filename.str().c_str(), ios::out);
  output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header[Energy.MatrixCounter] << endl;
  for(int j=0;j<Energy.RowCounter[Energy.MatrixCounter];j++) {
    for(int k=0;k<Energy.ColumnCounter[Energy.MatrixCounter];k++)
      output << scientific << Energy.Matrix[ Energy.MatrixCounter ][j][k] << "\t";
    output << endl;
  }
  output << endl;
  
  output << endl << "Total Forces" << endl << "===============" << endl << Force.Header[Force.MatrixCounter] << endl;
  for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
    for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
      output << scientific << Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
    output << endl;
  }
  output << endl;

  if (!NoHessian) {
    output << endl << "Total Hessian" << endl << "===============" << endl << Hessian.Header[Hessian.MatrixCounter] << endl;
    for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
      for(int k=0;k<Hessian.ColumnCounter[Hessian.MatrixCounter];k++)
        output << scientific << Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
      output << endl;
    }
    output << endl;
  }

  if (periode != NULL) { // also look for PAS values
    output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header[Hessian.MatrixCounter] << endl;
    for(int j=0;j<Shielding.RowCounter[Shielding.MatrixCounter];j++) {
      for(int k=0;k<Shielding.ColumnCounter[Shielding.MatrixCounter];k++)
        output << scientific << Shielding.Matrix[ Shielding.MatrixCounter ][j][k] << "\t";
      output << endl;
    }
    output << endl;
  
    output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header[ShieldingPAS.MatrixCounter] << endl;
    for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
      for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
        output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t";
      output << endl;
    }
    output << endl;

    output << endl << "Total Chis" << endl << "===============" << endl << Chi.Header[Chi.MatrixCounter] << endl;
    for(int j=0;j<Chi.RowCounter[Chi.MatrixCounter];j++) {
      for(int k=0;k<Chi.ColumnCounter[Chi.MatrixCounter];k++)
        output << scientific << Chi.Matrix[ Chi.MatrixCounter ][j][k] << "\t";
      output << endl;
    }
    output << endl;
  
    output << endl << "Total Chis PAS" << endl << "===============" << endl << ChiPAS.Header[ChiPAS.MatrixCounter] << endl;
    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
      for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
        output << scientific << ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t";
      output << endl;
    }
    output << endl;
  }
  
  if (!NoTime) {
    output << endl << "Total Times" << endl << "===============" << endl << Time.Header[Time.MatrixCounter] << endl;
    for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
      for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
        output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
      }
      output << endl;
    }
    output << endl;
  }
  output.close();
  if (!NoTime)
    for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++)
      Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];

  // +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++

  cout << "Analyzing ..." << endl;

  // ======================================= Creating the data files ==============================================================

  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
  // +++++++++++++++++++++++++++++++++++++++ Plotting Delta Simtime vs Bond Order
  if (!NoTime) {
    if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
    if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
    for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
      for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
        Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
      }
    counter = 0;
    output << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl;
    output2 << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl;
    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
      for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
        for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
          for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
            Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
          }
      counter += KeySet.FragmentsPerOrder[BondOrder];
      output << BondOrder+1 << "\t" << counter;
      output2 << BondOrder+1 << "\t" << counter;
      for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
        output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
        if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
          output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
        else
          output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
      }
      output << endl;
      output2 << endl;
    }
    output.close();
    output2.close();
  }

  if (!NoHessian) {
    // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in hessian to full QM
    if (!CreateDataDeltaHessianOrderPerAtom(Hessian, HessianFragments, KeySet, argv[3], "DeltaHessian_xx-Order", "Plot of error between approximated hessian and full hessian versus the Bond Order", datum)) return 1;

    if (!CreateDataDeltaFrobeniusOrderPerAtom(Hessian, HessianFragments, KeySet, argv[3], "DeltaFrobeniusHessian_xx-Order", "Plot of error between approximated hessian and full hessian in the frobenius norm versus the Bond Order", datum)) return 1;

    // ++++++++++++++++++++++++++++++++++++++Plotting Hessian vs. Order
    if (!CreateDataHessianOrderPerAtom(HessianFragments, KeySet, argv[3], "Hessian_xx-Order", "Plot of approximated hessian versus the Bond Order", datum)) return 1;
    if (!AppendOutputFile(output, argv[3], "Hessian_xx-Order.dat" )) return false;
    output << endl << "# Full" << endl;
    for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
      output << j << "\t";
      for(int k=0;k<Hessian.ColumnCounter[Force.MatrixCounter];k++)
        output << scientific <<  Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
      output << endl;
    }
    output.close();
  }

  // +++++++++++++++++++++++++++++++++++++++ Plotting shieldings
  if (periode != NULL) { // also look for PAS values
    if (!CreateDataDeltaForcesOrderPerAtom(ShieldingPAS, ShieldingPASFragments, KeySet, argv[3], "DeltaShieldingsPAS-Order", "Plot of error between approximated shieldings and full shieldings versus the Bond Order", datum)) return 1;
    if (!CreateDataForcesOrderPerAtom(ShieldingPASFragments, KeySet, argv[3], "ShieldingsPAS-Order", "Plot of approximated shieldings versus the Bond Order", datum)) return 1;
    if (!AppendOutputFile(output, argv[3], "ShieldingsPAS-Order.dat" )) return false;
    output << endl << "# Full" << endl;
    for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
      output << j << "\t";
      for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
        output << scientific <<  ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
      output << endl;
    }
    output.close();
    if (!CreateDataDeltaForcesOrderPerAtom(ChiPAS, ChiPASFragments, KeySet, argv[3], "DeltaChisPAS-Order", "Plot of error between approximated Chis and full Chis versus the Bond Order", datum)) return 1;
    if (!CreateDataForcesOrderPerAtom(ChiPASFragments, KeySet, argv[3], "ChisPAS-Order", "Plot of approximated Chis versus the Bond Order", datum)) return 1;
    if (!AppendOutputFile(output, argv[3], "ChisPAS-Order.dat" )) return false;
    output << endl << "# Full" << endl;
    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
      output << j << "\t";
      for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
        output << scientific <<  ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
      output << endl;
    }
    output.close();
  }


  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
  if (!CreateDataDeltaEnergyOrder(Energy, EnergyFragments, KeySet, argv[3], "DeltaEnergies-Order", "Plot of error between approximated and full energies energies versus the Bond Order", datum)) return 1;

  // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
  if (!CreateDataEnergyOrder(EnergyFragments, KeySet, argv[3], "Energies-Order", "Plot of approximated energies versus the Bond Order", datum)) return 1;

  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
  if (!CreateDataDeltaForcesOrderPerAtom(Force, ForceFragments, KeySet, argv[3], "DeltaForces-Order", "Plot of error between approximated forces and full forces versus the Bond Order", datum)) return 1;

  // min force
  if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMinForces-Order", "Plot of min error between approximated forces and full forces versus the Bond Order", datum, CreateMinimumForce)) return 1;

  // mean force
  if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMeanForces-Order", "Plot of mean error between approximated forces and full forces versus the Bond Order", datum, CreateMeanForce)) return 1;

  // max force
  if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMaxForces-Order", "Plot of max error between approximated forces and full forces versus the Bond Order", datum, CreateMaximumForce)) return 1;

  // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
  if (!CreateDataForcesOrderPerAtom(ForceFragments, KeySet, argv[3], "Forces-Order", "Plot of approximated forces versus the Bond Order", datum)) return 1;
  if (!AppendOutputFile(output, argv[3], "Forces-Order.dat" )) return false;
  output << endl << "# Full" << endl;
  for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
    output << j << "\t";
    for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
      output << scientific <<  Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
    output << endl;
  }
  output.close();
  // min force
  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MinForces-Order", "Plot of min approximated forces versus the Bond Order", datum, CreateMinimumForce)) return 1;

  // mean force
  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MeanForces-Order", "Plot of mean approximated forces versus the Bond Order", datum, CreateMeanForce)) return 1;

  // max force
  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MaxForces-Order", "Plot of max approximated forces versus the Bond Order", datum, CreateMaximumForce)) return 1;

  // ++++++++++++++++++++++++++++++++++++++Plotting vector sum (should be 0) vs. bond order
  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "VectorSum-Order", "Plot of vector sum of the approximated forces versus the Bond Order", datum, CreateVectorSumForce)) return 1;

  // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
  if (!CreateDataFragment(EnergyFragments, KeySet, argv[3], "Energies-Fragment", "Plot of fragment energy versus the Fragment No", datum, CreateEnergy)) return 1;
  if (!CreateDataFragment(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", "Plot of fragment energy of each Fragment No vs. Bond Order", datum, CreateEnergy)) return 1;
  if (!CreateDataFragmentOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", "Plot of maximum of fragment energy vs. Bond Order", datum, CreateMaxFragmentOrder)) return 1;
  if (!CreateDataFragmentOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", "Plot of minimum of fragment energy vs. Bond Order", datum, CreateMinFragmentOrder)) return 1;

  // +++++++++++++++++++++++++++++++Ploting min/mean/max forces for each fragment
  // min force
  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MinForces-Fragment", "Plot of min approximated forces versus the Fragment No", datum, CreateMinimumForce)) return 1;

  // mean force
  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MeanForces-Fragment", "Plot of mean approximated forces versus the Fragment No", datum, CreateMeanForce)) return 1;

  // max force
  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MaxForces-Fragment", "Plot of max approximated forces versus the Fragment No", datum, CreateMaximumForce)) return 1;

  // +++++++++++++++++++++++++++++++Ploting min/mean/max forces for each fragment per order
  // min force
  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MinForces-FragmentOrder", "Plot of min approximated forces of each Fragment No vs. Bond Order", datum, CreateMinimumForce)) return 1;

  // mean force
  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MeanForces-FragmentOrder", "Plot of mean approximated forces of each Fragment No vs. Bond Order", datum, CreateMeanForce)) return 1;

  // max force
  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MaxForces-FragmentOrder", "Plot of max approximated forces of each Fragment No vs. Bond Order", datum, CreateMaximumForce)) return 1;

  // ======================================= Creating the plot files ==============================================================

  Orderxrange << "[1:" << KeySet.Order << "]";
  Fragmentxrange << "[0:" << KeySet.FragmentCounter+1 << "]";
  yrange.str("[1e-8:1e+1]");
  
  if (!NoTime) {
    // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
    if (!CreatePlotOrder(Time, KeySet, argv[3], "SimTime-Order", 1, "below", "y", "",  1, 1, "bond order k", "Evaluation time [s]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1;
  }
  
  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
  if (!CreatePlotOrder(Energy, KeySet, argv[3], "DeltaEnergies-Order", 1, "outside", "y", "",  1, 1, "bond order k", "absolute error in energy [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;

  // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
  if (!CreatePlotOrder(Energy, KeySet, argv[3], "Energies-Order", 1, "outside", "", "",  1, 1, "bond order k", "approximate energy [Ht]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1;

  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
  yrange.str("[1e-8:1e+0]");
  // min force
  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMinForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in min force [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;

  // mean force
  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMeanForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in mean force [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsFirstForceValuePlotLine)) return 1;

  // max force
  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMaxForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in max force [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;

  // min/mean/max comparison for total force
  if(!OpenOutputFile(output, argv[3], "DeltaMinMeanMaxTotalForce-Order.pyx")) return 1;
  CreatePlotHeader(output, "DeltaMinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL,  1, 1, "bond order k", "absolute error in total forces [Ht/a.u.]");
  output << "plot " << Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;
  output << "'DeltaMinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
  output << "'DeltaMeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
  output << "'DeltaMaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
  output.close();

  // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
  // min force
  if (!CreatePlotOrder(Force, KeySet, argv[3], "MinForces-Order", 2, "bottom right", "y", "",  1, 1, "bond order k", "absolute approximated min force [Ht/a.u.]", Orderxrange.str().c_str(), "", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;

  // mean force
  if (!CreatePlotOrder(Force, KeySet, argv[3], "MeanForces-Order", 2, "bottom right", "y", "",  1, 1, "bond order k", "absolute approximated mean force [Ht/a.u.]", Orderxrange.str().c_str(), "", "1" , "with linespoints", AbsFirstForceValuePlotLine)) return 1;

  // max force
  if (!CreatePlotOrder(Force, KeySet, argv[3], "MaxForces-Order", 2, "bottom right", "y", "",  1, 1, "bond order k", "absolute approximated max force [Ht/a.u.]", Orderxrange.str().c_str(), "", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;

  // min/mean/max comparison for total force
  if(!OpenOutputFile(output, argv[3],"MinMeanMaxTotalForce-Order.pyx")) return 1;
  CreatePlotHeader(output, "MinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL, 1, 1, "bond order k", "absolute total force [Ht/a.u.]");
  output << "plot "<< Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;
  output << "'MinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
  output << "'MeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
  output << "'MaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
  output.close();

  // ++++++++++++++++++++++++++++++++++++++Plotting vector sum vs. Order

  if (!CreatePlotOrder(Force, KeySet, argv[3], "VectorSum-Order", 2, "bottom right", "y" ,"", 1, 1, "bond order k", "vector sum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), "", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;

  // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
  yrange.str("");
  yrange << "[" << EnergyFragments.FindMinValue() << ":" << EnergyFragments.FindMaxValue() << "]";
  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-Fragment", 5, "below", "y", "", 1, 5, "fragment number", "Energies of each fragment [Ht]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with points", AbsEnergyPlotLine)) return 1;
  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Energies of each fragment [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with points", AbsEnergyPlotLine)) return 1;
  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Maximum fragment energy [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Minimum fragment energy [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;

  // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment
  yrange.str("");
  yrange << "[" << ForceFragments.FindMinValue() << ":" << ForceFragments.FindMaxValue()<< "]";
  // min
  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "minimum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;

  // mean
  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "mean of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;

  // max
  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "maximum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;

  // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment per bond order
  // min
  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "minimum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;

  // mean
  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "mean of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;

  // max
  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "maximum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;

  // +++++++++++++++++++++++++++++++=Ploting approximated and true shielding for each atom
  if (periode != NULL) { // also look for PAS values
    if(!OpenOutputFile(output, argv[3], "ShieldingsPAS-Order.pyx")) return 1;
    if(!OpenOutputFile(output2, argv[3], "DeltaShieldingsPAS-Order.pyx")) return 1;
    CreatePlotHeader(output, "ShieldingsPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical shielding value [ppm]");
    CreatePlotHeader(output2, "DeltaShieldingsPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical shielding value [ppm]");
    double step=0.8/KeySet.Order;
    output << "set boxwidth " << step << endl;
    output << "plot [0:" << ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter]+10 << "]\\" << endl;
    output2 << "plot [0:" << ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter]+10 << "]\\" << endl;
    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
      output << "'ShieldingsPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
      output2 << "'DeltaShieldingsPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
      if (BondOrder-1 != KeySet.Order)
        output2 << ", \\" << endl;
    }
    output << "'ShieldingsPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
    output2.close();  

    if(!OpenOutputFile(output, argv[3], "ChisPAS-Order.pyx")) return 1;
    if(!OpenOutputFile(output2, argv[3], "DeltaChisPAS-Order.pyx")) return 1;
    CreatePlotHeader(output, "ChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
    CreatePlotHeader(output2, "DeltaChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
    output << "set boxwidth " << step << endl;
    output << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
    output2 << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
      output << "'ChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
      output2 << "'DeltaChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
      if (BondOrder-1 != KeySet.Order)
        output2 << ", \\" << endl;
    }
    output << "'ChisPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
    output.close();  
    output2.close();  

    if(!OpenOutputFile(output, argv[3], "ChisPAS-Order.pyx")) return 1;
    if(!OpenOutputFile(output2, argv[3], "DeltaChisPAS-Order.pyx")) return 1;
    CreatePlotHeader(output, "ChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
    CreatePlotHeader(output2, "DeltaChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
    output << "set boxwidth " << step << endl;
    output << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
    output2 << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
      output << "'ChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
      output2 << "'DeltaChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
      if (BondOrder-1 != KeySet.Order)
        output2 << ", \\" << endl;
    }
    output << "'ChisPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
    output.close();  
    output2.close();  
  }

  // create Makefile
  if(!OpenOutputFile(output, argv[3], "Makefile")) return 1;
  output << "PYX = $(shell ls *.pyx)" << endl << endl;
  output << "EPS = $(PYX:.pyx=.eps)" << endl << endl;
  output << "%.eps: %.pyx" << endl;
  output << "\t~/build/pyxplot/pyxplot $<" << endl << endl;
  output << "all: $(EPS)" << endl << endl;
  output << ".PHONY: clean" << endl;
  output << "clean:" << endl;
  output << "\trm -rf $(EPS)" << endl;
  output.close();

  // ++++++++++++++++ exit ++++++++++++++++++++++++++++++++++
  delete(periode);
  Free(&dir);
  cout << "done." << endl;
  return 0;
};

//============================ END ===========================