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 <cstring>
|
---|
11 |
|
---|
12 | #include "datacreator.hpp"
|
---|
13 | #include "helpers.hpp"
|
---|
14 | #include "memoryallocator.hpp"
|
---|
15 | #include "parser.hpp"
|
---|
16 | #include "periodentafel.hpp"
|
---|
17 |
|
---|
18 | // include config.h
|
---|
19 | #ifdef HAVE_CONFIG_H
|
---|
20 | #include <config.h>
|
---|
21 | #endif
|
---|
22 |
|
---|
23 |
|
---|
24 | //============================== MAIN =============================
|
---|
25 |
|
---|
26 | int main(int argc, char **argv)
|
---|
27 | {
|
---|
28 | periodentafel *periode = NULL; // and a period table of all elements
|
---|
29 | EnergyMatrix Energy;
|
---|
30 | EnergyMatrix EnergyFragments;
|
---|
31 | ForceMatrix Force;
|
---|
32 | ForceMatrix ForceFragments;
|
---|
33 | HessianMatrix Hessian;
|
---|
34 | HessianMatrix HessianFragments;
|
---|
35 | EnergyMatrix Hcorrection;
|
---|
36 | EnergyMatrix HcorrectionFragments;
|
---|
37 | ForceMatrix Shielding;
|
---|
38 | ForceMatrix ShieldingPAS;
|
---|
39 | ForceMatrix Chi;
|
---|
40 | ForceMatrix ChiPAS;
|
---|
41 | EnergyMatrix Time;
|
---|
42 | ForceMatrix ShieldingFragments;
|
---|
43 | ForceMatrix ShieldingPASFragments;
|
---|
44 | ForceMatrix ChiFragments;
|
---|
45 | ForceMatrix ChiPASFragments;
|
---|
46 | KeySetsContainer KeySet;
|
---|
47 | ofstream output;
|
---|
48 | ofstream output2;
|
---|
49 | ofstream output3;
|
---|
50 | ofstream output4;
|
---|
51 | ifstream input;
|
---|
52 | stringstream filename;
|
---|
53 | time_t t = time(NULL);
|
---|
54 | struct tm *ts = localtime(&t);
|
---|
55 | char *datum = asctime(ts);
|
---|
56 | stringstream Orderxrange;
|
---|
57 | stringstream Fragmentxrange;
|
---|
58 | stringstream yrange;
|
---|
59 | char *dir = NULL;
|
---|
60 | bool NoHessian = false;
|
---|
61 | bool NoTime = false;
|
---|
62 | bool NoHCorrection = true;
|
---|
63 | int counter = 0;
|
---|
64 |
|
---|
65 | DoLog(0) && (Log() << Verbose(0) << "ANOVA Analyzer" << endl);
|
---|
66 | DoLog(0) && (Log() << Verbose(0) << "==============" << endl);
|
---|
67 |
|
---|
68 | // Get the command line options
|
---|
69 | if (argc < 4) {
|
---|
70 | DoLog(0) && (Log() << Verbose(0) << "Usage: " << argv[0] << " <inputdir> <prefix> <outputdir> [elementsdb]" << endl);
|
---|
71 | DoLog(0) && (Log() << Verbose(0) << "<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);
|
---|
72 | DoLog(0) && (Log() << Verbose(0) << "<prefix>\tprefix of energy and forces file." << endl);
|
---|
73 | DoLog(0) && (Log() << Verbose(0) << "<outputdir>\tcreated plotfiles and datafiles are placed into this directory " << endl);
|
---|
74 | DoLog(0) && (Log() << Verbose(0) << "[elementsdb]\tpath to elements database, needed for shieldings." << endl);
|
---|
75 | return 1;
|
---|
76 | } else {
|
---|
77 | dir = new char[strlen(argv[2]) + 2];
|
---|
78 | strcpy(dir, "/");
|
---|
79 | strcat(dir, argv[2]);
|
---|
80 | }
|
---|
81 |
|
---|
82 | if (argc > 4) {
|
---|
83 | DoLog(0) && (Log() << Verbose(0) << "Loading periodentafel." << endl);
|
---|
84 | periode = new periodentafel;
|
---|
85 | periode->LoadPeriodentafel(argv[4]);
|
---|
86 | }
|
---|
87 |
|
---|
88 | // Test the given directory
|
---|
89 | if (!TestParams(argc, argv))
|
---|
90 | return 1;
|
---|
91 |
|
---|
92 | // +++++++++++++++++ PARSING +++++++++++++++++++++++++++++++
|
---|
93 |
|
---|
94 | // ------------- Parse through all Fragment subdirs --------
|
---|
95 | if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
|
---|
96 | if (!Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0)) {
|
---|
97 | NoHCorrection = true;
|
---|
98 | DoeLog(2) && (eLog()<< Verbose(2) << "No HCorrection file found, skipping these." << endl);
|
---|
99 | }
|
---|
100 |
|
---|
101 | if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
|
---|
102 | if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix,0,0)) {
|
---|
103 | NoHessian = true;
|
---|
104 | DoeLog(2) && (eLog()<< Verbose(2) << "No Hessian file found, skipping these." << endl);
|
---|
105 | }
|
---|
106 | if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) {
|
---|
107 | NoTime = true;
|
---|
108 | DoeLog(2) && (eLog()<< Verbose(2) << "No speed file found, skipping these." << endl);
|
---|
109 | }
|
---|
110 | if (periode != NULL) { // also look for PAS values
|
---|
111 | if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
|
---|
112 | if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
|
---|
113 | if (!Chi.ParseFragmentMatrix(argv[1], dir, ChiSuffix, 1, 0)) return 1;
|
---|
114 | if (!ChiPAS.ParseFragmentMatrix(argv[1], dir, ChiPASSuffix, 1, 0)) return 1;
|
---|
115 | }
|
---|
116 |
|
---|
117 | // ---------- Parse the TE Factors into an array -----------------
|
---|
118 | if (!Energy.ParseIndices()) return 1;
|
---|
119 | if (!NoHCorrection) Hcorrection.ParseIndices();
|
---|
120 |
|
---|
121 | // ---------- Parse the Force indices into an array ---------------
|
---|
122 | if (!Force.ParseIndices(argv[1])) return 1;
|
---|
123 | if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
|
---|
124 | if (!ForceFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
|
---|
125 |
|
---|
126 | // ---------- Parse hessian indices into an array -----------------
|
---|
127 | if (!NoHessian) {
|
---|
128 | if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
|
---|
129 | if (!HessianFragments.AllocateMatrix(Hessian.Header, Hessian.MatrixCounter, Hessian.RowCounter, Hessian.ColumnCounter)) return 1;
|
---|
130 | if (!HessianFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
|
---|
131 | }
|
---|
132 |
|
---|
133 | // ---------- Parse the shielding indices into an array ---------------
|
---|
134 | if (periode != NULL) { // also look for PAS values
|
---|
135 | if(!Shielding.ParseIndices(argv[1])) return 1;
|
---|
136 | if(!ShieldingPAS.ParseIndices(argv[1])) return 1;
|
---|
137 | if (!ShieldingFragments.AllocateMatrix(Shielding.Header, Shielding.MatrixCounter, Shielding.RowCounter, Shielding.ColumnCounter)) return 1;
|
---|
138 | if (!ShieldingPASFragments.AllocateMatrix(ShieldingPAS.Header, ShieldingPAS.MatrixCounter, ShieldingPAS.RowCounter, ShieldingPAS.ColumnCounter)) return 1;
|
---|
139 | if(!ShieldingFragments.ParseIndices(argv[1])) return 1;
|
---|
140 | if(!ShieldingPASFragments.ParseIndices(argv[1])) return 1;
|
---|
141 | if(!Chi.ParseIndices(argv[1])) return 1;
|
---|
142 | if(!ChiPAS.ParseIndices(argv[1])) return 1;
|
---|
143 | if (!ChiFragments.AllocateMatrix(Chi.Header, Chi.MatrixCounter, Chi.RowCounter, Chi.ColumnCounter)) return 1;
|
---|
144 | if (!ChiPASFragments.AllocateMatrix(ChiPAS.Header, ChiPAS.MatrixCounter, ChiPAS.RowCounter, ChiPAS.ColumnCounter)) return 1;
|
---|
145 | if(!ChiFragments.ParseIndices(argv[1])) return 1;
|
---|
146 | if(!ChiPASFragments.ParseIndices(argv[1])) return 1;
|
---|
147 | }
|
---|
148 |
|
---|
149 | // ---------- Parse the KeySets into an array ---------------
|
---|
150 | if (!KeySet.ParseKeySets(argv[1], Force.RowCounter, Force.MatrixCounter)) return 1;
|
---|
151 | if (!KeySet.ParseManyBodyTerms()) return 1;
|
---|
152 |
|
---|
153 | // ---------- Parse fragment files created by 'joiner' into an array -------------
|
---|
154 | if (!EnergyFragments.ParseFragmentMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
|
---|
155 | if (!NoHCorrection)
|
---|
156 | HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
|
---|
157 | if (!ForceFragments.ParseFragmentMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
|
---|
158 | if (!NoHessian)
|
---|
159 | if (!HessianFragments.ParseFragmentMatrix(argv[1], dir, HessianFragmentSuffix,0,0)) return 1;
|
---|
160 | if (periode != NULL) { // also look for PAS values
|
---|
161 | if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingFragmentSuffix, 1, 0)) return 1;
|
---|
162 | if (!ShieldingPASFragments.ParseFragmentMatrix(argv[1], dir, ShieldingPASFragmentSuffix, 1, 0)) return 1;
|
---|
163 | if (!ChiFragments.ParseFragmentMatrix(argv[1], dir, ChiFragmentSuffix, 1, 0)) return 1;
|
---|
164 | if (!ChiPASFragments.ParseFragmentMatrix(argv[1], dir, ChiPASFragmentSuffix, 1, 0)) return 1;
|
---|
165 | }
|
---|
166 |
|
---|
167 | // +++++++++++++++ TESTING ++++++++++++++++++++++++++++++
|
---|
168 |
|
---|
169 | // print energy and forces to file
|
---|
170 | filename.str("");
|
---|
171 | filename << argv[3] << "/" << "energy-forces.all";
|
---|
172 | output.open(filename.str().c_str(), ios::out);
|
---|
173 | output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header[Energy.MatrixCounter] << endl;
|
---|
174 | for(int j=0;j<Energy.RowCounter[Energy.MatrixCounter];j++) {
|
---|
175 | for(int k=0;k<Energy.ColumnCounter[Energy.MatrixCounter];k++)
|
---|
176 | output << scientific << Energy.Matrix[ Energy.MatrixCounter ][j][k] << "\t";
|
---|
177 | output << endl;
|
---|
178 | }
|
---|
179 | output << endl;
|
---|
180 |
|
---|
181 | output << endl << "Total Forces" << endl << "===============" << endl << Force.Header[Force.MatrixCounter] << endl;
|
---|
182 | for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
|
---|
183 | for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
|
---|
184 | output << scientific << Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
|
---|
185 | output << endl;
|
---|
186 | }
|
---|
187 | output << endl;
|
---|
188 |
|
---|
189 | if (!NoHessian) {
|
---|
190 | output << endl << "Total Hessian" << endl << "===============" << endl << Hessian.Header[Hessian.MatrixCounter] << endl;
|
---|
191 | for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
|
---|
192 | for(int k=0;k<Hessian.ColumnCounter[Hessian.MatrixCounter];k++)
|
---|
193 | output << scientific << Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
|
---|
194 | output << endl;
|
---|
195 | }
|
---|
196 | output << endl;
|
---|
197 | }
|
---|
198 |
|
---|
199 | if (periode != NULL) { // also look for PAS values
|
---|
200 | output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header[Hessian.MatrixCounter] << endl;
|
---|
201 | for(int j=0;j<Shielding.RowCounter[Shielding.MatrixCounter];j++) {
|
---|
202 | for(int k=0;k<Shielding.ColumnCounter[Shielding.MatrixCounter];k++)
|
---|
203 | output << scientific << Shielding.Matrix[ Shielding.MatrixCounter ][j][k] << "\t";
|
---|
204 | output << endl;
|
---|
205 | }
|
---|
206 | output << endl;
|
---|
207 |
|
---|
208 | output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header[ShieldingPAS.MatrixCounter] << endl;
|
---|
209 | for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
|
---|
210 | for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
|
---|
211 | output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t";
|
---|
212 | output << endl;
|
---|
213 | }
|
---|
214 | output << endl;
|
---|
215 |
|
---|
216 | output << endl << "Total Chis" << endl << "===============" << endl << Chi.Header[Chi.MatrixCounter] << endl;
|
---|
217 | for(int j=0;j<Chi.RowCounter[Chi.MatrixCounter];j++) {
|
---|
218 | for(int k=0;k<Chi.ColumnCounter[Chi.MatrixCounter];k++)
|
---|
219 | output << scientific << Chi.Matrix[ Chi.MatrixCounter ][j][k] << "\t";
|
---|
220 | output << endl;
|
---|
221 | }
|
---|
222 | output << endl;
|
---|
223 |
|
---|
224 | output << endl << "Total Chis PAS" << endl << "===============" << endl << ChiPAS.Header[ChiPAS.MatrixCounter] << endl;
|
---|
225 | for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
|
---|
226 | for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
|
---|
227 | output << scientific << ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t";
|
---|
228 | output << endl;
|
---|
229 | }
|
---|
230 | output << endl;
|
---|
231 | }
|
---|
232 |
|
---|
233 | if (!NoTime) {
|
---|
234 | output << endl << "Total Times" << endl << "===============" << endl << Time.Header[Time.MatrixCounter] << endl;
|
---|
235 | for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
|
---|
236 | for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
|
---|
237 | output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
|
---|
238 | }
|
---|
239 | output << endl;
|
---|
240 | }
|
---|
241 | output << endl;
|
---|
242 | }
|
---|
243 | output.close();
|
---|
244 | if (!NoTime)
|
---|
245 | for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++)
|
---|
246 | Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];
|
---|
247 |
|
---|
248 | // +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++
|
---|
249 |
|
---|
250 | DoLog(0) && (Log() << Verbose(0) << "Analyzing ..." << endl);
|
---|
251 |
|
---|
252 | // ======================================= Creating the data files ==============================================================
|
---|
253 |
|
---|
254 | // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
|
---|
255 | // +++++++++++++++++++++++++++++++++++++++ Plotting Delta Simtime vs Bond Order
|
---|
256 | if (!NoTime) {
|
---|
257 | if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
|
---|
258 | if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
|
---|
259 | for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
|
---|
260 | for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
|
---|
261 | Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
|
---|
262 | }
|
---|
263 | counter = 0;
|
---|
264 | output << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl;
|
---|
265 | output2 << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl;
|
---|
266 | for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
|
---|
267 | for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
|
---|
268 | for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
|
---|
269 | for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
|
---|
270 | Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
|
---|
271 | }
|
---|
272 | counter += KeySet.FragmentsPerOrder[BondOrder];
|
---|
273 | output << BondOrder+1 << "\t" << counter;
|
---|
274 | output2 << BondOrder+1 << "\t" << counter;
|
---|
275 | for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
|
---|
276 | output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
|
---|
277 | if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
|
---|
278 | output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
|
---|
279 | else
|
---|
280 | output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
|
---|
281 | }
|
---|
282 | output << endl;
|
---|
283 | output2 << endl;
|
---|
284 | }
|
---|
285 | output.close();
|
---|
286 | output2.close();
|
---|
287 | }
|
---|
288 |
|
---|
289 | if (!NoHessian) {
|
---|
290 | // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in hessian to full QM
|
---|
291 | if (!CreateDataDeltaHessianOrderPerAtom(Hessian, HessianFragments, KeySet, argv[3], "DeltaHessian_xx-Order", "Plot of error between approximated hessian and full hessian versus the Bond Order", datum)) return 1;
|
---|
292 |
|
---|
293 | if (!CreateDataDeltaFrobeniusOrderPerAtom(Hessian, HessianFragments, KeySet, argv[3], "DeltaFrobeniusHessian_xx-Order", "Plot of error between approximated hessian and full hessian in the frobenius norm versus the Bond Order", datum)) return 1;
|
---|
294 |
|
---|
295 | // ++++++++++++++++++++++++++++++++++++++Plotting Hessian vs. Order
|
---|
296 | if (!CreateDataHessianOrderPerAtom(HessianFragments, KeySet, argv[3], "Hessian_xx-Order", "Plot of approximated hessian versus the Bond Order", datum)) return 1;
|
---|
297 | if (!AppendOutputFile(output, argv[3], "Hessian_xx-Order.dat" )) return false;
|
---|
298 | output << endl << "# Full" << endl;
|
---|
299 | for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
|
---|
300 | output << j << "\t";
|
---|
301 | for(int k=0;k<Hessian.ColumnCounter[Force.MatrixCounter];k++)
|
---|
302 | output << scientific << Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
|
---|
303 | output << endl;
|
---|
304 | }
|
---|
305 | output.close();
|
---|
306 | }
|
---|
307 |
|
---|
308 | // +++++++++++++++++++++++++++++++++++++++ Plotting shieldings
|
---|
309 | if (periode != NULL) { // also look for PAS values
|
---|
310 | 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;
|
---|
311 | if (!CreateDataForcesOrderPerAtom(ShieldingPASFragments, KeySet, argv[3], "ShieldingsPAS-Order", "Plot of approximated shieldings versus the Bond Order", datum)) return 1;
|
---|
312 | if (!AppendOutputFile(output, argv[3], "ShieldingsPAS-Order.dat" )) return false;
|
---|
313 | output << endl << "# Full" << endl;
|
---|
314 | for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
|
---|
315 | output << j << "\t";
|
---|
316 | for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
|
---|
317 | output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
|
---|
318 | output << endl;
|
---|
319 | }
|
---|
320 | output.close();
|
---|
321 | if (!CreateDataDeltaForcesOrderPerAtom(ChiPAS, ChiPASFragments, KeySet, argv[3], "DeltaChisPAS-Order", "Plot of error between approximated Chis and full Chis versus the Bond Order", datum)) return 1;
|
---|
322 | if (!CreateDataForcesOrderPerAtom(ChiPASFragments, KeySet, argv[3], "ChisPAS-Order", "Plot of approximated Chis versus the Bond Order", datum)) return 1;
|
---|
323 | if (!AppendOutputFile(output, argv[3], "ChisPAS-Order.dat" )) return false;
|
---|
324 | output << endl << "# Full" << endl;
|
---|
325 | for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
|
---|
326 | output << j << "\t";
|
---|
327 | for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
|
---|
328 | output << scientific << ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
|
---|
329 | output << endl;
|
---|
330 | }
|
---|
331 | output.close();
|
---|
332 | }
|
---|
333 |
|
---|
334 |
|
---|
335 | // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
|
---|
336 | 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;
|
---|
337 |
|
---|
338 | // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
|
---|
339 | if (!CreateDataEnergyOrder(EnergyFragments, KeySet, argv[3], "Energies-Order", "Plot of approximated energies versus the Bond Order", datum)) return 1;
|
---|
340 |
|
---|
341 | // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
|
---|
342 | 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;
|
---|
343 |
|
---|
344 | // min force
|
---|
345 | 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;
|
---|
346 |
|
---|
347 | // mean force
|
---|
348 | 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;
|
---|
349 |
|
---|
350 | // max force
|
---|
351 | 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;
|
---|
352 |
|
---|
353 | // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
|
---|
354 | if (!CreateDataForcesOrderPerAtom(ForceFragments, KeySet, argv[3], "Forces-Order", "Plot of approximated forces versus the Bond Order", datum)) return 1;
|
---|
355 | if (!AppendOutputFile(output, argv[3], "Forces-Order.dat" )) return false;
|
---|
356 | output << endl << "# Full" << endl;
|
---|
357 | for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
|
---|
358 | output << j << "\t";
|
---|
359 | for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
|
---|
360 | output << scientific << Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
|
---|
361 | output << endl;
|
---|
362 | }
|
---|
363 | output.close();
|
---|
364 | // min force
|
---|
365 | if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MinForces-Order", "Plot of min approximated forces versus the Bond Order", datum, CreateMinimumForce)) return 1;
|
---|
366 |
|
---|
367 | // mean force
|
---|
368 | if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MeanForces-Order", "Plot of mean approximated forces versus the Bond Order", datum, CreateMeanForce)) return 1;
|
---|
369 |
|
---|
370 | // max force
|
---|
371 | if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MaxForces-Order", "Plot of max approximated forces versus the Bond Order", datum, CreateMaximumForce)) return 1;
|
---|
372 |
|
---|
373 | // ++++++++++++++++++++++++++++++++++++++Plotting vector sum (should be 0) vs. bond order
|
---|
374 | if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "VectorSum-Order", "Plot of vector sum of the approximated forces versus the Bond Order", datum, CreateVectorSumForce)) return 1;
|
---|
375 |
|
---|
376 | // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
|
---|
377 | if (!CreateDataFragment(EnergyFragments, KeySet, argv[3], "Energies-Fragment", "Plot of fragment energy versus the Fragment No", datum, CreateEnergy)) return 1;
|
---|
378 | if (!CreateDataFragment(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", "Plot of fragment energy of each Fragment No vs. Bond Order", datum, CreateEnergy)) return 1;
|
---|
379 | if (!CreateDataFragmentOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", "Plot of maximum of fragment energy vs. Bond Order", datum, CreateMaxFragmentOrder)) return 1;
|
---|
380 | if (!CreateDataFragmentOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", "Plot of minimum of fragment energy vs. Bond Order", datum, CreateMinFragmentOrder)) return 1;
|
---|
381 |
|
---|
382 | // +++++++++++++++++++++++++++++++Ploting min/mean/max forces for each fragment
|
---|
383 | // min force
|
---|
384 | if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MinForces-Fragment", "Plot of min approximated forces versus the Fragment No", datum, CreateMinimumForce)) return 1;
|
---|
385 |
|
---|
386 | // mean force
|
---|
387 | if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MeanForces-Fragment", "Plot of mean approximated forces versus the Fragment No", datum, CreateMeanForce)) return 1;
|
---|
388 |
|
---|
389 | // max force
|
---|
390 | if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MaxForces-Fragment", "Plot of max approximated forces versus the Fragment No", datum, CreateMaximumForce)) return 1;
|
---|
391 |
|
---|
392 | // +++++++++++++++++++++++++++++++Ploting min/mean/max forces for each fragment per order
|
---|
393 | // min force
|
---|
394 | if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MinForces-FragmentOrder", "Plot of min approximated forces of each Fragment No vs. Bond Order", datum, CreateMinimumForce)) return 1;
|
---|
395 |
|
---|
396 | // mean force
|
---|
397 | if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MeanForces-FragmentOrder", "Plot of mean approximated forces of each Fragment No vs. Bond Order", datum, CreateMeanForce)) return 1;
|
---|
398 |
|
---|
399 | // max force
|
---|
400 | if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MaxForces-FragmentOrder", "Plot of max approximated forces of each Fragment No vs. Bond Order", datum, CreateMaximumForce)) return 1;
|
---|
401 |
|
---|
402 | // ======================================= Creating the plot files ==============================================================
|
---|
403 |
|
---|
404 | Orderxrange << "[1:" << KeySet.Order << "]";
|
---|
405 | Fragmentxrange << "[0:" << KeySet.FragmentCounter+1 << "]";
|
---|
406 | yrange.str("[1e-8:1e+1]");
|
---|
407 |
|
---|
408 | if (!NoTime) {
|
---|
409 | // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
|
---|
410 | 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;
|
---|
411 | }
|
---|
412 |
|
---|
413 | // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
|
---|
414 | 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;
|
---|
415 |
|
---|
416 | // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
|
---|
417 | 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;
|
---|
418 |
|
---|
419 | // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
|
---|
420 | yrange.str("[1e-8:1e+0]");
|
---|
421 | // min force
|
---|
422 | 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;
|
---|
423 |
|
---|
424 | // mean force
|
---|
425 | 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;
|
---|
426 |
|
---|
427 | // max force
|
---|
428 | 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;
|
---|
429 |
|
---|
430 | // min/mean/max comparison for total force
|
---|
431 | if(!OpenOutputFile(output, argv[3], "DeltaMinMeanMaxTotalForce-Order.pyx")) return 1;
|
---|
432 | CreatePlotHeader(output, "DeltaMinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL, 1, 1, "bond order k", "absolute error in total forces [Ht/a.u.]");
|
---|
433 | output << "plot " << Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;
|
---|
434 | output << "'DeltaMinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
|
---|
435 | output << "'DeltaMeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
|
---|
436 | output << "'DeltaMaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
|
---|
437 | output.close();
|
---|
438 |
|
---|
439 | // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
|
---|
440 | // min force
|
---|
441 | 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;
|
---|
442 |
|
---|
443 | // mean force
|
---|
444 | 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;
|
---|
445 |
|
---|
446 | // max force
|
---|
447 | 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;
|
---|
448 |
|
---|
449 | // min/mean/max comparison for total force
|
---|
450 | if(!OpenOutputFile(output, argv[3],"MinMeanMaxTotalForce-Order.pyx")) return 1;
|
---|
451 | CreatePlotHeader(output, "MinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL, 1, 1, "bond order k", "absolute total force [Ht/a.u.]");
|
---|
452 | output << "plot "<< Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;
|
---|
453 | output << "'MinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
|
---|
454 | output << "'MeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
|
---|
455 | output << "'MaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
|
---|
456 | output.close();
|
---|
457 |
|
---|
458 | // ++++++++++++++++++++++++++++++++++++++Plotting vector sum vs. Order
|
---|
459 |
|
---|
460 | 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;
|
---|
461 |
|
---|
462 | // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
|
---|
463 | yrange.str("");
|
---|
464 | yrange << "[" << EnergyFragments.FindMinValue() << ":" << EnergyFragments.FindMaxValue() << "]";
|
---|
465 | 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;
|
---|
466 | 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;
|
---|
467 | 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;
|
---|
468 | 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;
|
---|
469 |
|
---|
470 | // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment
|
---|
471 | yrange.str("");
|
---|
472 | yrange << "[" << ForceFragments.FindMinValue() << ":" << ForceFragments.FindMaxValue()<< "]";
|
---|
473 | // min
|
---|
474 | 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;
|
---|
475 |
|
---|
476 | // mean
|
---|
477 | 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;
|
---|
478 |
|
---|
479 | // max
|
---|
480 | 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;
|
---|
481 |
|
---|
482 | // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment per bond order
|
---|
483 | // min
|
---|
484 | 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;
|
---|
485 |
|
---|
486 | // mean
|
---|
487 | 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;
|
---|
488 |
|
---|
489 | // max
|
---|
490 | 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;
|
---|
491 |
|
---|
492 | // +++++++++++++++++++++++++++++++=Ploting approximated and true shielding for each atom
|
---|
493 | if (periode != NULL) { // also look for PAS values
|
---|
494 | if(!OpenOutputFile(output, argv[3], "ShieldingsPAS-Order.pyx")) return 1;
|
---|
495 | if(!OpenOutputFile(output2, argv[3], "DeltaShieldingsPAS-Order.pyx")) return 1;
|
---|
496 | CreatePlotHeader(output, "ShieldingsPAS-Order", 1, "top right", NULL, NULL, 1, 5, "nuclei index", "iso chemical shielding value [ppm]");
|
---|
497 | CreatePlotHeader(output2, "DeltaShieldingsPAS-Order", 1, "top right", NULL, NULL, 1, 5, "nuclei index", "iso chemical shielding value [ppm]");
|
---|
498 | double step=0.8/KeySet.Order;
|
---|
499 | output << "set boxwidth " << step << endl;
|
---|
500 | output << "plot [0:" << ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter]+10 << "]\\" << endl;
|
---|
501 | output2 << "plot [0:" << ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter]+10 << "]\\" << endl;
|
---|
502 | for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
|
---|
503 | output << "'ShieldingsPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
|
---|
504 | output2 << "'DeltaShieldingsPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
|
---|
505 | if (BondOrder-1 != KeySet.Order)
|
---|
506 | output2 << ", \\" << endl;
|
---|
507 | }
|
---|
508 | output << "'ShieldingsPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
|
---|
509 | output2.close();
|
---|
510 |
|
---|
511 | if(!OpenOutputFile(output, argv[3], "ChisPAS-Order.pyx")) return 1;
|
---|
512 | if(!OpenOutputFile(output2, argv[3], "DeltaChisPAS-Order.pyx")) return 1;
|
---|
513 | CreatePlotHeader(output, "ChisPAS-Order", 1, "top right", NULL, NULL, 1, 5, "nuclei index", "iso chemical Chi value [ppm]");
|
---|
514 | CreatePlotHeader(output2, "DeltaChisPAS-Order", 1, "top right", NULL, NULL, 1, 5, "nuclei index", "iso chemical Chi value [ppm]");
|
---|
515 | output << "set boxwidth " << step << endl;
|
---|
516 | output << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
|
---|
517 | output2 << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
|
---|
518 | for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
|
---|
519 | output << "'ChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
|
---|
520 | output2 << "'DeltaChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
|
---|
521 | if (BondOrder-1 != KeySet.Order)
|
---|
522 | output2 << ", \\" << endl;
|
---|
523 | }
|
---|
524 | output << "'ChisPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
|
---|
525 | output.close();
|
---|
526 | output2.close();
|
---|
527 |
|
---|
528 | if(!OpenOutputFile(output, argv[3], "ChisPAS-Order.pyx")) return 1;
|
---|
529 | if(!OpenOutputFile(output2, argv[3], "DeltaChisPAS-Order.pyx")) return 1;
|
---|
530 | CreatePlotHeader(output, "ChisPAS-Order", 1, "top right", NULL, NULL, 1, 5, "nuclei index", "iso chemical Chi value [ppm]");
|
---|
531 | CreatePlotHeader(output2, "DeltaChisPAS-Order", 1, "top right", NULL, NULL, 1, 5, "nuclei index", "iso chemical Chi value [ppm]");
|
---|
532 | output << "set boxwidth " << step << endl;
|
---|
533 | output << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
|
---|
534 | output2 << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
|
---|
535 | for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
|
---|
536 | output << "'ChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
|
---|
537 | output2 << "'DeltaChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
|
---|
538 | if (BondOrder-1 != KeySet.Order)
|
---|
539 | output2 << ", \\" << endl;
|
---|
540 | }
|
---|
541 | output << "'ChisPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
|
---|
542 | output.close();
|
---|
543 | output2.close();
|
---|
544 | }
|
---|
545 |
|
---|
546 | // create Makefile
|
---|
547 | if(!OpenOutputFile(output, argv[3], "Makefile")) return 1;
|
---|
548 | output << "PYX = $(shell ls *.pyx)" << endl << endl;
|
---|
549 | output << "EPS = $(PYX:.pyx=.eps)" << endl << endl;
|
---|
550 | output << "%.eps: %.pyx" << endl;
|
---|
551 | output << "\t~/build/pyxplot/pyxplot $<" << endl << endl;
|
---|
552 | output << "all: $(EPS)" << endl << endl;
|
---|
553 | output << ".PHONY: clean" << endl;
|
---|
554 | output << "clean:" << endl;
|
---|
555 | output << "\trm -rf $(EPS)" << endl;
|
---|
556 | output.close();
|
---|
557 |
|
---|
558 | // ++++++++++++++++ exit ++++++++++++++++++++++++++++++++++
|
---|
559 | delete(periode);
|
---|
560 | delete[](dir);
|
---|
561 | DoLog(0) && (Log() << Verbose(0) << "done." << endl);
|
---|
562 | return 0;
|
---|
563 | };
|
---|
564 |
|
---|
565 | //============================ END ===========================
|
---|