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