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