source: src/analyzer.cpp@ 1907a7

/** \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 "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 Hcorrection;
        ForceMatrix Force;
        ForceMatrix Shielding;
        ForceMatrix ShieldingPAS;
        EnergyMatrix Time;
        EnergyMatrix EnergyFragments;
        EnergyMatrix HcorrectionFragments;
        ForceMatrix ForceFragments;
        ForceMatrix ShieldingFragments;
        ForceMatrix ShieldingPASFragments;
        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 Hcorrected = true;
        double norm;
        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 = (char *) Malloc(sizeof(char)*(strlen(argv[2])+2), "main: *dir");
                strcpy(dir, "/");
                strcat(dir, argv[2]);
        }

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

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

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

        // ------------- Parse through all Fragment subdirs --------
        if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
        Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0);
        if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
        if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1;
        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;
        }

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

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

        // ---------- Parse the shielding indices into an array ---------------
        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;
        }

        // ---------- 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 (Hcorrected) HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
        if (!ForceFragments.ParseFragmentMatrix(argv[1], dir, ForceFragmentSuffix,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;
        }

        // +++++++++++++++ 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 << endl;
        for(int j=0;j<Energy.RowCounter[Energy.MatrixCounter];j++) {
                for(int k=0;k<Energy.ColumnCounter;k++)
                        output << scientific << Energy.Matrix[ Energy.MatrixCounter ][j][k] << "\t";
                output << endl;
        }
        output << endl;

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

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

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

        output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
        for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
                for(int k=0;k<Time.ColumnCounter;k++) {
                        output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
                }
                output << endl;
        }
        output << endl;
        output.close();
        for(int k=0;k<Time.ColumnCounter;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 (!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;k--;) {
                        Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
                }
        counter = 0;
        output << "#Order\tFrag.No.\t" << Time.Header << endl;
        output2 << "#Order\tFrag.No.\t" << Time.Header << endl;
        for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
                for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
                        for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
                                for(int k=Time.ColumnCounter;k--;) {
                                        Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
                                }
                counter += KeySet.FragmentsPerOrder[BondOrder];
                output << BondOrder+1 << "\t" << counter;
                output2 << BondOrder+1 << "\t" << counter;
                for(int k=0;k<Time.ColumnCounter;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();

        // +++++++++++++++++++++++++++++++++++++++ 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;k++)
                                output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.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;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]");

        // +++++++++++++++++++++++++++++++++++++++ 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;
                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((void **)&dir, "main: *dir");
        cout << "done." << endl;
        return 0;
};

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