Ignore:
Timestamp:
Jul 23, 2009, 2:21:57 PM (16 years ago)
Author:
Frederik Heber <heber@…>
Children:
c3a303
Parents:
c95b69
Message:

fixed indentation from tabs to two spaces.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/analyzer.cpp

    rc95b69 ra048fa  
    2323int main(int argc, char **argv)
    2424{
    25         periodentafel *periode = NULL; // and a period table of all elements
    26         EnergyMatrix Energy;
    27         EnergyMatrix Hcorrection;
    28         ForceMatrix Force;
    29         ForceMatrix Shielding;
    30         ForceMatrix ShieldingPAS;
    31         EnergyMatrix Time;
    32         EnergyMatrix EnergyFragments;
    33         EnergyMatrix HcorrectionFragments;
    34         ForceMatrix ForceFragments;
    35         ForceMatrix ShieldingFragments;
    36         ForceMatrix ShieldingPASFragments;
    37         KeySetsContainer KeySet;
    38         ofstream output;
    39         ofstream output2;
    40         ofstream output3;
    41         ofstream output4;
    42         ifstream input;
    43         stringstream filename;
    44         time_t t = time(NULL);
    45         struct tm *ts = localtime(&t);
    46         char *datum = asctime(ts);
    47         stringstream Orderxrange;
    48         stringstream Fragmentxrange;
    49         stringstream yrange;
    50         char *dir = NULL;
    51         bool Hcorrected = true;
    52         int counter;
    53 
    54         cout << "ANOVA Analyzer" << endl;
    55         cout << "==============" << endl;
    56 
    57         // Get the command line options
    58         if (argc < 4) {
    59                 cout << "Usage: " << argv[0] << " <inputdir> <prefix> <outputdir> [elementsdb]" << endl;
    60                 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;
    61                 cout << "<prefix>\tprefix of energy and forces file." << endl;
    62                 cout << "<outputdir>\tcreated plotfiles and datafiles are placed into this directory " << endl;
    63                 cout << "[elementsdb]\tpath to elements database, needed for shieldings." << endl;
    64                 return 1;
    65         } else {
    66                 dir = (char *) Malloc(sizeof(char)*(strlen(argv[2])+2), "main: *dir");
    67                 strcpy(dir, "/");
    68                 strcat(dir, argv[2]);
    69         }
    70 
    71         if (argc > 4) {
    72                 cout << "Loading periodentafel." << endl;
    73                 periode = new periodentafel;
    74                 periode->LoadPeriodentafel(argv[4]);
    75         }
    76 
    77         // Test the given directory
    78         if (!TestParams(argc, argv))
    79                 return 1;
    80 
    81         // +++++++++++++++++ PARSING +++++++++++++++++++++++++++++++
    82 
    83         // ------------- Parse through all Fragment subdirs --------
    84         if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
    85         Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0);
    86         if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
    87         if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1;
    88         if (periode != NULL) { // also look for PAS values
    89                 if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
    90                 if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
    91         }
    92 
    93         // ---------- Parse the TE Factors into an array -----------------
    94         if (!Energy.ParseIndices()) return 1;
    95         if (Hcorrected) Hcorrection.ParseIndices();
    96 
    97         // ---------- Parse the Force indices into an array ---------------
    98         if (!Force.ParseIndices(argv[1])) return 1;
    99         if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
    100         if (!ForceFragments.ParseIndices(argv[1])) return 1;
    101 
    102         // ---------- Parse the shielding indices into an array ---------------
    103         if (periode != NULL) { // also look for PAS values
    104                 if(!Shielding.ParseIndices(argv[1])) return 1;
    105                 if(!ShieldingPAS.ParseIndices(argv[1])) return 1;
    106                 if (!ShieldingFragments.AllocateMatrix(Shielding.Header, Shielding.MatrixCounter, Shielding.RowCounter, Shielding.ColumnCounter)) return 1;
    107                 if (!ShieldingPASFragments.AllocateMatrix(ShieldingPAS.Header, ShieldingPAS.MatrixCounter, ShieldingPAS.RowCounter, ShieldingPAS.ColumnCounter)) return 1;
    108                 if(!ShieldingFragments.ParseIndices(argv[1])) return 1;
    109                 if(!ShieldingPASFragments.ParseIndices(argv[1])) return 1;
    110         }
    111 
    112         // ---------- Parse the KeySets into an array ---------------
    113         if (!KeySet.ParseKeySets(argv[1], Force.RowCounter, Force.MatrixCounter)) return 1;
    114         if (!KeySet.ParseManyBodyTerms()) return 1;
    115 
    116         // ---------- Parse fragment files created by 'joiner' into an array -------------
    117         if (!EnergyFragments.ParseFragmentMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
    118         if (Hcorrected) HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
    119         if (!ForceFragments.ParseFragmentMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
    120         if (periode != NULL) { // also look for PAS values
    121                 if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingFragmentSuffix, 1, 0)) return 1;
    122                 if (!ShieldingPASFragments.ParseFragmentMatrix(argv[1], dir, ShieldingPASFragmentSuffix, 1, 0)) return 1;
    123         }
    124 
    125         // +++++++++++++++ TESTING ++++++++++++++++++++++++++++++
    126 
    127         // print energy and forces to file
    128         filename.str("");
    129         filename << argv[3] << "/" << "energy-forces.all";
    130         output.open(filename.str().c_str(), ios::out);
    131         output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header << endl;
    132         for(int j=0;j<Energy.RowCounter[Energy.MatrixCounter];j++) {
    133                 for(int k=0;k<Energy.ColumnCounter;k++)
    134                         output << scientific << Energy.Matrix[ Energy.MatrixCounter ][j][k] << "\t";
    135                 output << endl;
    136         }
    137         output << endl;
    138 
    139         output << endl << "Total Forces" << endl << "===============" << endl << Force.Header << endl;
    140         for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
    141                 for(int k=0;k<Force.ColumnCounter;k++)
    142                         output << scientific << Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
    143                 output << endl;
    144         }
    145         output << endl;
    146 
    147         if (periode != NULL) { // also look for PAS values
    148                 output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header << endl;
    149                 for(int j=0;j<Shielding.RowCounter[Shielding.MatrixCounter];j++) {
    150                         for(int k=0;k<Shielding.ColumnCounter;k++)
    151                                 output << scientific << Shielding.Matrix[ Shielding.MatrixCounter ][j][k] << "\t";
    152                         output << endl;
    153                 }
    154                 output << endl;
    155 
    156                 output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header << endl;
    157                 for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
    158                         for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
    159                                 output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t";
    160                         output << endl;
    161                 }
    162                 output << endl;
    163         }
    164 
    165         output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
    166         for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
    167                 for(int k=0;k<Time.ColumnCounter;k++) {
    168                         output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
    169                 }
    170                 output << endl;
    171         }
    172         output << endl;
    173         output.close();
    174         for(int k=0;k<Time.ColumnCounter;k++)
    175                 Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];
    176 
    177         // +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++
    178 
    179         cout << "Analyzing ..." << endl;
    180 
    181         // ======================================= Creating the data files ==============================================================
    182 
    183         // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
    184         // +++++++++++++++++++++++++++++++++++++++ Plotting Delta Simtime vs Bond Order
    185         if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
    186         if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
    187         for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
    188                 for(int k=Time.ColumnCounter;k--;) {
    189                         Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
    190                 }
    191         counter = 0;
    192         output << "#Order\tFrag.No.\t" << Time.Header << endl;
    193         output2 << "#Order\tFrag.No.\t" << Time.Header << endl;
    194         for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    195                 for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
    196                         for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
    197                                 for(int k=Time.ColumnCounter;k--;) {
    198                                         Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
    199                                 }
    200                 counter += KeySet.FragmentsPerOrder[BondOrder];
    201                 output << BondOrder+1 << "\t" << counter;
    202                 output2 << BondOrder+1 << "\t" << counter;
    203                 for(int k=0;k<Time.ColumnCounter;k++) {
    204                         output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
    205                         if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
    206                                 output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
    207                         else
    208                                 output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
    209                 }
    210                 output << endl;
    211                 output2 << endl;
    212         }
    213         output.close();
    214         output2.close();
    215 
    216         // +++++++++++++++++++++++++++++++++++++++ Plotting shieldings
    217 
    218         if (periode != NULL) { // also look for PAS values
    219                 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;
    220                 if (!CreateDataForcesOrderPerAtom(ShieldingPASFragments, KeySet, argv[3], "ShieldingsPAS-Order", "Plot of approximated shieldings versus the Bond Order", datum)) return 1;
    221                 if (!AppendOutputFile(output, argv[3], "ShieldingsPAS-Order.dat" )) return false;
    222                 output << endl << "# Full" << endl;
    223                 for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
    224                         output << j << "\t";
    225                         for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
    226                                 output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
    227                         output << endl;
    228                 }
    229         }
    230         output.close();
    231 
    232 
    233         // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
    234         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;
    235 
    236         // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
    237         if (!CreateDataEnergyOrder(EnergyFragments, KeySet, argv[3], "Energies-Order", "Plot of approximated energies versus the Bond Order", datum)) return 1;
    238 
    239         // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
    240         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;
    241 
    242         // min force
    243         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;
    244 
    245         // mean force
    246         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;
    247 
    248         // max force
    249         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;
    250 
    251         // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
    252         if (!CreateDataForcesOrderPerAtom(ForceFragments, KeySet, argv[3], "Forces-Order", "Plot of approximated forces versus the Bond Order", datum)) return 1;
    253         if (!AppendOutputFile(output, argv[3], "Forces-Order.dat" )) return false;
    254         output << endl << "# Full" << endl;
    255         for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
    256                 output << j << "\t";
    257                 for(int k=0;k<Force.ColumnCounter;k++)
    258                         output << scientific << Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
    259                 output << endl;
    260         }
    261         output.close();
    262         // min force
    263         if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MinForces-Order", "Plot of min approximated forces versus the Bond Order", datum, CreateMinimumForce)) return 1;
    264 
    265         // mean force
    266         if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MeanForces-Order", "Plot of mean approximated forces versus the Bond Order", datum, CreateMeanForce)) return 1;
    267 
    268         // max force
    269         if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MaxForces-Order", "Plot of max approximated forces versus the Bond Order", datum, CreateMaximumForce)) return 1;
    270 
    271         // ++++++++++++++++++++++++++++++++++++++Plotting vector sum (should be 0) vs. bond order
    272         if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "VectorSum-Order", "Plot of vector sum of the approximated forces versus the Bond Order", datum, CreateVectorSumForce)) return 1;
    273 
    274         // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
    275         if (!CreateDataFragment(EnergyFragments, KeySet, argv[3], "Energies-Fragment", "Plot of fragment energy versus the Fragment No", datum, CreateEnergy)) return 1;
    276         if (!CreateDataFragment(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", "Plot of fragment energy of each Fragment No vs. Bond Order", datum, CreateEnergy)) return 1;
    277         if (!CreateDataFragmentOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", "Plot of maximum of fragment energy vs. Bond Order", datum, CreateMaxFragmentOrder)) return 1;
    278         if (!CreateDataFragmentOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", "Plot of minimum of fragment energy vs. Bond Order", datum, CreateMinFragmentOrder)) return 1;
    279 
    280         // +++++++++++++++++++++++++++++++Ploting min/mean/max forces for each fragment
    281         // min force
    282         if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MinForces-Fragment", "Plot of min approximated forces versus the Fragment No", datum, CreateMinimumForce)) return 1;
    283 
    284         // mean force
    285         if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MeanForces-Fragment", "Plot of mean approximated forces versus the Fragment No", datum, CreateMeanForce)) return 1;
    286 
    287         // max force
    288         if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MaxForces-Fragment", "Plot of max approximated forces versus the Fragment No", datum, CreateMaximumForce)) return 1;
    289 
    290         // +++++++++++++++++++++++++++++++Ploting min/mean/max forces for each fragment per order
    291         // min force
    292         if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MinForces-FragmentOrder", "Plot of min approximated forces of each Fragment No vs. Bond Order", datum, CreateMinimumForce)) return 1;
    293 
    294         // mean force
    295         if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MeanForces-FragmentOrder", "Plot of mean approximated forces of each Fragment No vs. Bond Order", datum, CreateMeanForce)) return 1;
    296 
    297         // max force
    298         if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MaxForces-FragmentOrder", "Plot of max approximated forces of each Fragment No vs. Bond Order", datum, CreateMaximumForce)) return 1;
    299 
    300         // ======================================= Creating the plot files ==============================================================
    301 
    302         Orderxrange << "[1:" << KeySet.Order << "]";
    303         Fragmentxrange << "[0:" << KeySet.FragmentCounter+1 << "]";
    304         yrange.str("[1e-8:1e+1]");
    305 
    306         // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
    307         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;
    308 
    309         // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
    310         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;
    311 
    312         // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
    313         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;
    314 
    315         // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
    316         yrange.str("[1e-8:1e+0]");
    317         // min force
    318         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;
    319 
    320         // mean force
    321         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;
    322 
    323         // max force
    324         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;
    325 
    326         // min/mean/max comparison for total force
    327         if(!OpenOutputFile(output, argv[3], "DeltaMinMeanMaxTotalForce-Order.pyx")) return 1;
    328         CreatePlotHeader(output, "DeltaMinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL,        1, 1, "bond order k", "absolute error in total forces [Ht/a.u.]");
    329         output << "plot " << Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;
    330         output << "'DeltaMinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
    331         output << "'DeltaMeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
    332         output << "'DeltaMaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
    333         output.close();
    334 
    335         // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
    336         // min force
    337         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;
    338 
    339         // mean force
    340         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;
    341 
    342         // max force
    343         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;
    344 
    345         // min/mean/max comparison for total force
    346         if(!OpenOutputFile(output, argv[3],"MinMeanMaxTotalForce-Order.pyx")) return 1;
    347         CreatePlotHeader(output, "MinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL, 1, 1, "bond order k", "absolute total force [Ht/a.u.]");
    348         output << "plot "<< Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;
    349         output << "'MinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
    350         output << "'MeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
    351         output << "'MaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
    352         output.close();
    353 
    354         // ++++++++++++++++++++++++++++++++++++++Plotting vector sum vs. Order
    355 
    356         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;
    357 
    358         // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
    359         yrange.str("");
    360         yrange << "[" << EnergyFragments.FindMinValue() << ":" << EnergyFragments.FindMaxValue() << "]";
    361         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;
    362         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;
    363         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;
    364         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;
    365 
    366         // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment
    367         yrange.str("");
    368         yrange << "[" << ForceFragments.FindMinValue() << ":" << ForceFragments.FindMaxValue()<< "]";
    369         // min
    370         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;
    371 
    372         // mean
    373         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;
    374 
    375         // max
    376         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;
    377 
    378         // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment per bond order
    379         // min
    380         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;
    381 
    382         // mean
    383         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;
    384 
    385         // max
    386         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;
    387 
    388         // +++++++++++++++++++++++++++++++=Ploting approximated and true shielding for each atom
    389         if (periode != NULL) { // also look for PAS values
    390                 if(!OpenOutputFile(output, argv[3], "ShieldingsPAS-Order.pyx")) return 1;
    391                 if(!OpenOutputFile(output2, argv[3], "DeltaShieldingsPAS-Order.pyx")) return 1;
    392                 CreatePlotHeader(output, "ShieldingsPAS-Order", 1, "top right", NULL, NULL,     1, 5, "nuclei index", "iso chemical shielding value [ppm]");
    393                 CreatePlotHeader(output2, "DeltaShieldingsPAS-Order", 1, "top right", NULL, NULL,       1, 5, "nuclei index", "iso chemical shielding value [ppm]");
    394                 double step=0.8/KeySet.Order;
    395                 output << "set boxwidth " << step << endl;
    396                 output << "plot [0:" << ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter]+10 << "]\\" << endl;
    397                 output2 << "plot [0:" << ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter]+10 << "]\\" << endl;
    398                 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    399                         output << "'ShieldingsPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
    400                         output2 << "'DeltaShieldingsPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
    401                         if (BondOrder-1 != KeySet.Order)
    402                                 output2 << ", \\" << endl;
    403                 }
    404                 output << "'ShieldingsPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
    405                 output.close();
    406                 output2.close();
    407         }
    408 
    409         // create Makefile
    410         if(!OpenOutputFile(output, argv[3], "Makefile")) return 1;
    411         output << "PYX = $(shell ls *.pyx)" << endl << endl;
    412         output << "EPS = $(PYX:.pyx=.eps)" << endl << endl;
    413         output << "%.eps: %.pyx" << endl;
    414         output << "\t~/build/pyxplot/pyxplot $<" << endl << endl;
    415         output << "all: $(EPS)" << endl << endl;
    416         output << ".PHONY: clean" << endl;
    417         output << "clean:" << endl;
    418         output << "\trm -rf $(EPS)" << endl;
    419         output.close();
    420 
    421         // ++++++++++++++++ exit ++++++++++++++++++++++++++++++++++
    422         delete(periode);
    423         Free((void **)&dir, "main: *dir");
    424         cout << "done." << endl;
    425         return 0;
     25  periodentafel *periode = NULL; // and a period table of all elements
     26  EnergyMatrix Energy;
     27  EnergyMatrix Hcorrection;
     28  ForceMatrix Force;
     29  ForceMatrix Shielding;
     30  ForceMatrix ShieldingPAS;
     31  EnergyMatrix Time;
     32  EnergyMatrix EnergyFragments;
     33  EnergyMatrix HcorrectionFragments;
     34  ForceMatrix ForceFragments;
     35  ForceMatrix ShieldingFragments;
     36  ForceMatrix ShieldingPASFragments;
     37  KeySetsContainer KeySet;
     38  ofstream output;
     39  ofstream output2;
     40  ofstream output3;
     41  ofstream output4;
     42  ifstream input;
     43  stringstream filename;
     44  time_t t = time(NULL);
     45  struct tm *ts = localtime(&t);
     46  char *datum = asctime(ts);
     47  stringstream Orderxrange;
     48  stringstream Fragmentxrange;
     49  stringstream yrange;
     50  char *dir = NULL;
     51  bool Hcorrected = true;
     52  int counter;
     53
     54  cout << "ANOVA Analyzer" << endl;
     55  cout << "==============" << endl;
     56
     57  // Get the command line options
     58  if (argc < 4) {
     59    cout << "Usage: " << argv[0] << " <inputdir> <prefix> <outputdir> [elementsdb]" << endl;
     60    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;
     61    cout << "<prefix>\tprefix of energy and forces file." << endl;
     62    cout << "<outputdir>\tcreated plotfiles and datafiles are placed into this directory " << endl;
     63    cout << "[elementsdb]\tpath to elements database, needed for shieldings." << endl;
     64    return 1;
     65  } else {
     66    dir = (char *) Malloc(sizeof(char)*(strlen(argv[2])+2), "main: *dir");
     67    strcpy(dir, "/");
     68    strcat(dir, argv[2]);
     69  }
     70
     71  if (argc > 4) {
     72    cout << "Loading periodentafel." << endl;
     73    periode = new periodentafel;
     74    periode->LoadPeriodentafel(argv[4]);
     75  }
     76
     77  // Test the given directory
     78  if (!TestParams(argc, argv))
     79    return 1;
     80
     81  // +++++++++++++++++ PARSING +++++++++++++++++++++++++++++++
     82
     83  // ------------- Parse through all Fragment subdirs --------
     84  if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
     85  Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0);
     86  if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
     87  if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1;
     88  if (periode != NULL) { // also look for PAS values
     89    if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
     90    if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
     91  }
     92
     93  // ---------- Parse the TE Factors into an array -----------------
     94  if (!Energy.ParseIndices()) return 1;
     95  if (Hcorrected) Hcorrection.ParseIndices();
     96
     97  // ---------- Parse the Force indices into an array ---------------
     98  if (!Force.ParseIndices(argv[1])) return 1;
     99  if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
     100  if (!ForceFragments.ParseIndices(argv[1])) return 1;
     101
     102  // ---------- Parse the shielding indices into an array ---------------
     103  if (periode != NULL) { // also look for PAS values
     104    if(!Shielding.ParseIndices(argv[1])) return 1;
     105    if(!ShieldingPAS.ParseIndices(argv[1])) return 1;
     106    if (!ShieldingFragments.AllocateMatrix(Shielding.Header, Shielding.MatrixCounter, Shielding.RowCounter, Shielding.ColumnCounter)) return 1;
     107    if (!ShieldingPASFragments.AllocateMatrix(ShieldingPAS.Header, ShieldingPAS.MatrixCounter, ShieldingPAS.RowCounter, ShieldingPAS.ColumnCounter)) return 1;
     108    if(!ShieldingFragments.ParseIndices(argv[1])) return 1;
     109    if(!ShieldingPASFragments.ParseIndices(argv[1])) return 1;
     110  }
     111
     112  // ---------- Parse the KeySets into an array ---------------
     113  if (!KeySet.ParseKeySets(argv[1], Force.RowCounter, Force.MatrixCounter)) return 1;
     114  if (!KeySet.ParseManyBodyTerms()) return 1;
     115
     116  // ---------- Parse fragment files created by 'joiner' into an array -------------
     117  if (!EnergyFragments.ParseFragmentMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
     118  if (Hcorrected) HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
     119  if (!ForceFragments.ParseFragmentMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
     120  if (periode != NULL) { // also look for PAS values
     121    if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingFragmentSuffix, 1, 0)) return 1;
     122    if (!ShieldingPASFragments.ParseFragmentMatrix(argv[1], dir, ShieldingPASFragmentSuffix, 1, 0)) return 1;
     123  }
     124
     125  // +++++++++++++++ TESTING ++++++++++++++++++++++++++++++
     126
     127  // print energy and forces to file
     128  filename.str("");
     129  filename << argv[3] << "/" << "energy-forces.all";
     130  output.open(filename.str().c_str(), ios::out);
     131  output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header << endl;
     132  for(int j=0;j<Energy.RowCounter[Energy.MatrixCounter];j++) {
     133    for(int k=0;k<Energy.ColumnCounter;k++)
     134      output << scientific << Energy.Matrix[ Energy.MatrixCounter ][j][k] << "\t";
     135    output << endl;
     136  }
     137  output << endl;
     138
     139  output << endl << "Total Forces" << endl << "===============" << endl << Force.Header << endl;
     140  for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
     141    for(int k=0;k<Force.ColumnCounter;k++)
     142      output << scientific << Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
     143    output << endl;
     144  }
     145  output << endl;
     146
     147  if (periode != NULL) { // also look for PAS values
     148    output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header << endl;
     149    for(int j=0;j<Shielding.RowCounter[Shielding.MatrixCounter];j++) {
     150      for(int k=0;k<Shielding.ColumnCounter;k++)
     151        output << scientific << Shielding.Matrix[ Shielding.MatrixCounter ][j][k] << "\t";
     152      output << endl;
     153    }
     154    output << endl;
     155
     156    output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header << endl;
     157    for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
     158      for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
     159        output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t";
     160      output << endl;
     161    }
     162    output << endl;
     163  }
     164
     165  output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
     166  for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
     167    for(int k=0;k<Time.ColumnCounter;k++) {
     168      output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
     169    }
     170    output << endl;
     171  }
     172  output << endl;
     173  output.close();
     174  for(int k=0;k<Time.ColumnCounter;k++)
     175    Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];
     176
     177  // +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++
     178
     179  cout << "Analyzing ..." << endl;
     180
     181  // ======================================= Creating the data files ==============================================================
     182
     183  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
     184  // +++++++++++++++++++++++++++++++++++++++ Plotting Delta Simtime vs Bond Order
     185  if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
     186  if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
     187  for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
     188    for(int k=Time.ColumnCounter;k--;) {
     189      Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
     190    }
     191  counter = 0;
     192  output << "#Order\tFrag.No.\t" << Time.Header << endl;
     193  output2 << "#Order\tFrag.No.\t" << Time.Header << endl;
     194  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     195    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
     196      for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
     197        for(int k=Time.ColumnCounter;k--;) {
     198          Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     199        }
     200    counter += KeySet.FragmentsPerOrder[BondOrder];
     201    output << BondOrder+1 << "\t" << counter;
     202    output2 << BondOrder+1 << "\t" << counter;
     203    for(int k=0;k<Time.ColumnCounter;k++) {
     204      output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
     205      if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
     206        output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
     207      else
     208        output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
     209    }
     210    output << endl;
     211    output2 << endl;
     212  }
     213  output.close();
     214  output2.close();
     215
     216  // +++++++++++++++++++++++++++++++++++++++ Plotting shieldings
     217
     218  if (periode != NULL) { // also look for PAS values
     219    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;
     220    if (!CreateDataForcesOrderPerAtom(ShieldingPASFragments, KeySet, argv[3], "ShieldingsPAS-Order", "Plot of approximated shieldings versus the Bond Order", datum)) return 1;
     221    if (!AppendOutputFile(output, argv[3], "ShieldingsPAS-Order.dat" )) return false;
     222    output << endl << "# Full" << endl;
     223    for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
     224      output << j << "\t";
     225      for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
     226        output << scientific <<  ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
     227      output << endl;
     228    }
     229  }
     230  output.close();
     231
     232
     233  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
     234  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;
     235
     236  // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
     237  if (!CreateDataEnergyOrder(EnergyFragments, KeySet, argv[3], "Energies-Order", "Plot of approximated energies versus the Bond Order", datum)) return 1;
     238
     239  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
     240  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;
     241
     242  // min force
     243  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;
     244
     245  // mean force
     246  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;
     247
     248  // max force
     249  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;
     250
     251  // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
     252  if (!CreateDataForcesOrderPerAtom(ForceFragments, KeySet, argv[3], "Forces-Order", "Plot of approximated forces versus the Bond Order", datum)) return 1;
     253  if (!AppendOutputFile(output, argv[3], "Forces-Order.dat" )) return false;
     254  output << endl << "# Full" << endl;
     255  for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
     256    output << j << "\t";
     257    for(int k=0;k<Force.ColumnCounter;k++)
     258      output << scientific <<  Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
     259    output << endl;
     260  }
     261  output.close();
     262  // min force
     263  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MinForces-Order", "Plot of min approximated forces versus the Bond Order", datum, CreateMinimumForce)) return 1;
     264
     265  // mean force
     266  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MeanForces-Order", "Plot of mean approximated forces versus the Bond Order", datum, CreateMeanForce)) return 1;
     267
     268  // max force
     269  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MaxForces-Order", "Plot of max approximated forces versus the Bond Order", datum, CreateMaximumForce)) return 1;
     270
     271  // ++++++++++++++++++++++++++++++++++++++Plotting vector sum (should be 0) vs. bond order
     272  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "VectorSum-Order", "Plot of vector sum of the approximated forces versus the Bond Order", datum, CreateVectorSumForce)) return 1;
     273
     274  // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
     275  if (!CreateDataFragment(EnergyFragments, KeySet, argv[3], "Energies-Fragment", "Plot of fragment energy versus the Fragment No", datum, CreateEnergy)) return 1;
     276  if (!CreateDataFragment(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", "Plot of fragment energy of each Fragment No vs. Bond Order", datum, CreateEnergy)) return 1;
     277  if (!CreateDataFragmentOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", "Plot of maximum of fragment energy vs. Bond Order", datum, CreateMaxFragmentOrder)) return 1;
     278  if (!CreateDataFragmentOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", "Plot of minimum of fragment energy vs. Bond Order", datum, CreateMinFragmentOrder)) return 1;
     279
     280  // +++++++++++++++++++++++++++++++Ploting min/mean/max forces for each fragment
     281  // min force
     282  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MinForces-Fragment", "Plot of min approximated forces versus the Fragment No", datum, CreateMinimumForce)) return 1;
     283
     284  // mean force
     285  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MeanForces-Fragment", "Plot of mean approximated forces versus the Fragment No", datum, CreateMeanForce)) return 1;
     286
     287  // max force
     288  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MaxForces-Fragment", "Plot of max approximated forces versus the Fragment No", datum, CreateMaximumForce)) return 1;
     289
     290  // +++++++++++++++++++++++++++++++Ploting min/mean/max forces for each fragment per order
     291  // min force
     292  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MinForces-FragmentOrder", "Plot of min approximated forces of each Fragment No vs. Bond Order", datum, CreateMinimumForce)) return 1;
     293
     294  // mean force
     295  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MeanForces-FragmentOrder", "Plot of mean approximated forces of each Fragment No vs. Bond Order", datum, CreateMeanForce)) return 1;
     296
     297  // max force
     298  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MaxForces-FragmentOrder", "Plot of max approximated forces of each Fragment No vs. Bond Order", datum, CreateMaximumForce)) return 1;
     299
     300  // ======================================= Creating the plot files ==============================================================
     301
     302  Orderxrange << "[1:" << KeySet.Order << "]";
     303  Fragmentxrange << "[0:" << KeySet.FragmentCounter+1 << "]";
     304  yrange.str("[1e-8:1e+1]");
     305
     306  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
     307  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;
     308
     309  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
     310  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;
     311
     312  // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
     313  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;
     314
     315  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
     316  yrange.str("[1e-8:1e+0]");
     317  // min force
     318  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;
     319
     320  // mean force
     321  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;
     322
     323  // max force
     324  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;
     325
     326  // min/mean/max comparison for total force
     327  if(!OpenOutputFile(output, argv[3], "DeltaMinMeanMaxTotalForce-Order.pyx")) return 1;
     328  CreatePlotHeader(output, "DeltaMinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL,  1, 1, "bond order k", "absolute error in total forces [Ht/a.u.]");
     329  output << "plot " << Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;
     330  output << "'DeltaMinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
     331  output << "'DeltaMeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
     332  output << "'DeltaMaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
     333  output.close();
     334
     335  // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
     336  // min force
     337  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;
     338
     339  // mean force
     340  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;
     341
     342  // max force
     343  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;
     344
     345  // min/mean/max comparison for total force
     346  if(!OpenOutputFile(output, argv[3],"MinMeanMaxTotalForce-Order.pyx")) return 1;
     347  CreatePlotHeader(output, "MinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL, 1, 1, "bond order k", "absolute total force [Ht/a.u.]");
     348  output << "plot "<< Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;
     349  output << "'MinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
     350  output << "'MeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
     351  output << "'MaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
     352  output.close();
     353
     354  // ++++++++++++++++++++++++++++++++++++++Plotting vector sum vs. Order
     355
     356  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;
     357
     358  // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
     359  yrange.str("");
     360  yrange << "[" << EnergyFragments.FindMinValue() << ":" << EnergyFragments.FindMaxValue() << "]";
     361  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;
     362  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;
     363  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;
     364  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;
     365
     366  // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment
     367  yrange.str("");
     368  yrange << "[" << ForceFragments.FindMinValue() << ":" << ForceFragments.FindMaxValue()<< "]";
     369  // min
     370  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;
     371
     372  // mean
     373  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;
     374
     375  // max
     376  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;
     377
     378  // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment per bond order
     379  // min
     380  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;
     381
     382  // mean
     383  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;
     384
     385  // max
     386  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;
     387
     388  // +++++++++++++++++++++++++++++++=Ploting approximated and true shielding for each atom
     389  if (periode != NULL) { // also look for PAS values
     390    if(!OpenOutputFile(output, argv[3], "ShieldingsPAS-Order.pyx")) return 1;
     391    if(!OpenOutputFile(output2, argv[3], "DeltaShieldingsPAS-Order.pyx")) return 1;
     392    CreatePlotHeader(output, "ShieldingsPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical shielding value [ppm]");
     393    CreatePlotHeader(output2, "DeltaShieldingsPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical shielding value [ppm]");
     394    double step=0.8/KeySet.Order;
     395    output << "set boxwidth " << step << endl;
     396    output << "plot [0:" << ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter]+10 << "]\\" << endl;
     397    output2 << "plot [0:" << ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter]+10 << "]\\" << endl;
     398    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     399      output << "'ShieldingsPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
     400      output2 << "'DeltaShieldingsPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
     401      if (BondOrder-1 != KeySet.Order)
     402        output2 << ", \\" << endl;
     403    }
     404    output << "'ShieldingsPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
     405    output.close();
     406    output2.close();
     407  }
     408
     409  // create Makefile
     410  if(!OpenOutputFile(output, argv[3], "Makefile")) return 1;
     411  output << "PYX = $(shell ls *.pyx)" << endl << endl;
     412  output << "EPS = $(PYX:.pyx=.eps)" << endl << endl;
     413  output << "%.eps: %.pyx" << endl;
     414  output << "\t~/build/pyxplot/pyxplot $<" << endl << endl;
     415  output << "all: $(EPS)" << endl << endl;
     416  output << ".PHONY: clean" << endl;
     417  output << "clean:" << endl;
     418  output << "\trm -rf $(EPS)" << endl;
     419  output.close();
     420
     421  // ++++++++++++++++ exit ++++++++++++++++++++++++++++++++++
     422  delete(periode);
     423  Free((void **)&dir, "main: *dir");
     424  cout << "done." << endl;
     425  return 0;
    426426};
    427427
Note: See TracChangeset for help on using the changeset viewer.