source: src/datacreator.cpp@ 958457

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since 958457 was 390248, checked in by Frederik Heber <heber@…>, 16 years ago

BIG change: Hcorrection included and bugfix of analyzer (wrt forces)

  • Hcorrection is supplied via molecuilder which gives correction values if hydrogen atoms are to close and thus already feel each other's influence. This correction energy matrix is also parsed ans subtracted from the approximated energy matrix generated from the fragments. This impacts both on joiner, analyzer and molecuilder
  • Forces were not working which had multiple causes:
    • first, we stepped back down to absolute errors which rather gives a feel about convergence (1e-6 is simply small and neglegible)
    • there may have been bugs where (either in Force or in ForceFragments) adding up took place, this is cleaned up
    • more routines have been abstracted from analyzer into datacreator functions
      • the order of ions had changed, re-run with current molecuilder, calculation and subsequent joining/analyzing brought correct results finally
  • plots were bad still due to wrong axes (two more functions finding max/min values), use of wrong columns in plot file
  • Property mode set to 100644
File size: 27.4 KB
Line 
1/** \file datacreator.cpp
2 *
3 * Declarations of assisting functions in creating data and plot files.
4 *
5 */
6
7//============================ INCLUDES ===========================
8
9#include "datacreator.hpp"
10
11//=========================== FUNCTIONS============================
12
13/** Opens a file with \a *filename in \a *dir.
14 * \param output file handle on return
15 * \param *dir directory
16 * \param *filename name of file
17 * \return true if file has been opened
18 */
19bool OpenOutputFile(ofstream &output, const char *dir, const char *filename)
20{
21 stringstream name;
22 name << dir << "/" << filename;
23 output.open(name.str().c_str(), ios::out);
24 if (output == NULL) {
25 cout << "Unable to open " << name.str() << " for writing, is directory correct?" << endl;
26 return false;
27 }
28 return true;
29};
30
31/** Plots an energy vs. order.
32 * \param &Fragments EnergyMatrix class containing matrix values
33 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
34 * \param *prefix prefix in filename (without ending)
35 * \param *msg message to be place in first line as a comment
36 * \return true if file was written successfully
37 */
38bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
39{
40 stringstream filename;
41 ofstream output;
42
43 filename << prefix << ".dat";
44 if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
45 cout << msg << endl;
46 output << "# " << msg << ", created on " << datum;
47 output << "#Order\tFrag.No.\t" << Fragments.Header << endl;
48 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
49 for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
50 for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
51 for(int k=Fragments.ColumnCounter;k--;)
52 Fragments.Matrix[Fragments.MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
53 }
54 output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
55 for (int l=0;l<Fragments.ColumnCounter;l++)
56 output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
57 output << endl;
58 }
59 output.close();
60 return true;
61};
62
63/** Plots an energy error vs. order.
64 * \param &Energy EnergyMatrix class containing reference values (in MatrixCounter matrix)
65 * \param &Fragments EnergyMatrix class containing matrix values
66 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
67 * \param *prefix prefix in filename (without ending)
68 * \param *msg message to be place in first line as a comment
69 * \return true if file was written successfully
70 */
71bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
72{
73 stringstream filename;
74 ofstream output;
75
76 filename << prefix << ".dat";
77 if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
78 cout << msg << endl;
79 output << "# " << msg << ", created on " << datum;
80 output << "#Order\tFrag.No.\t" << Fragments.Header << endl;
81 Fragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0);
82 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
83 for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
84 for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
85 for(int k=Fragments.ColumnCounter;k--;)
86 Fragments.Matrix[Fragments.MatrixCounter][j][k] -= Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
87 }
88 output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
89 for (int l=0;l<Fragments.ColumnCounter;l++)
90 if (fabs(Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]) < MYEPSILON)
91 output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
92 else
93 output << scientific << "\t" << (Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l] / Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]);
94 output << endl;
95 }
96 output.close();
97 return true;
98};
99
100/** Plot forces vs. order.
101 * \param &Fragments ForceMatrix class containing matrix values
102 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
103 * \param *prefix prefix in filename (without ending)
104 * \param *msg message to be place in first line as a comment
105 * \param *datum current date and time
106 * \return true if file was written successfully
107 */
108bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
109{
110 stringstream filename;
111 ofstream output;
112
113 filename << prefix << ".dat";
114 if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
115 cout << msg << endl;
116 output << "# " << msg << ", created on " << datum;
117 output << "# Order\tFrag.No.\t" << Fragments.Header << endl;
118 Fragments.SetLastMatrix(0.,0);
119 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
120 Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
121 output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
122 CreateForce(Fragments, Fragments.MatrixCounter);
123 for (int l=0;l<Fragments.ColumnCounter;l++)
124 output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
125 output << endl;
126 }
127 output.close();
128 return true;
129};
130
131/** Plot forces error vs. order.
132 * \param &Force ForceMatrix containing reference values (in MatrixCounter matrix)
133 * \param &Fragments ForceMatrix class containing matrix values
134 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
135 * \param *prefix prefix in filename (without ending)
136 * \param *msg message to be place in first line as a comment
137 * \param *datum current date and time
138 * \return true if file was written successfully
139 */
140bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
141{
142 stringstream filename;
143 ofstream output;
144
145 filename << prefix << ".dat";
146 if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
147 cout << msg << endl;
148 output << "# " << msg << ", created on " << datum;
149 output << "# Order\tFrag.No.\t" << Fragments.Header << endl;
150 Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0);
151 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
152 Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
153 output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
154 CreateForce(Fragments, Fragments.MatrixCounter);
155 for (int l=0;l<Fragments.ColumnCounter;l++)
156 output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
157 output << endl;
158 }
159 output.close();
160 return true;
161};
162
163/** Plot forces error vs. vs atom vs. order.
164 * \param &Force ForceMatrix containing reference values (in MatrixCounter matrix)
165 * \param &Fragments ForceMatrix class containing matrix values
166 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
167 * \param *prefix prefix in filename (without ending)
168 * \param *msg message to be place in first line as a comment
169 * \param *datum current date and time
170 * \return true if file was written successfully
171 */
172bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
173{
174 stringstream filename;
175 ofstream output;
176 double norm = 0.;
177
178 filename << prefix << ".dat";
179 if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
180 cout << msg << endl;
181 output << "# " << msg << ", created on " << datum;
182 output << "# AtomNo\t" << Fragments.Header << endl;
183 Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter], 0);
184 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
185 //cout << "Current order is " << BondOrder << "." << endl;
186 Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
187 // errors per atom
188 output << "#Order\t" << BondOrder+1 << endl;
189 for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
190 output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
191 for (int l=0;l<Fragments.ColumnCounter;l++) {
192 if (((l+1) % 3) == 0) {
193 norm = 0.;
194 for (int m=0;m<NDIM;m++)
195 norm += Force.Matrix[Force.MatrixCounter][ j ][l+m]*Force.Matrix[Force.MatrixCounter][ j ][l+m];
196 norm = sqrt(norm);
197 }
198// if (norm < MYEPSILON)
199 output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
200// else
201// output << scientific << (Fragments.Matrix[Fragments.MatrixCounter][ j ][l] / norm) << "\t";
202 }
203 output << endl;
204 }
205 output << endl;
206 }
207 output.close();
208 return true;
209};
210
211/** Plot forces error vs. vs atom vs. order.
212 * \param &Fragments ForceMatrix class containing matrix values
213 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
214 * \param *prefix prefix in filename (without ending)
215 * \param *msg message to be place in first line as a comment
216 * \param *datum current date and time
217 * \return true if file was written successfully
218 */
219bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
220{
221 stringstream filename;
222 ofstream output;
223
224 filename << prefix << ".dat";
225 if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
226 cout << msg << endl;
227 output << "# " << msg << ", created on " << datum;
228 output << "# AtomNo\t" << Fragments.Header << endl;
229 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
230 //cout << "Current order is " << BondOrder << "." << endl;
231 Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
232 // errors per atom
233 output << "#Order\t" << BondOrder+1 << endl;
234 for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
235 output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
236 for (int l=0;l<Fragments.ColumnCounter;l++)
237 output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
238 output << endl;
239 }
240 output << endl;
241 }
242 output.close();
243 return true;
244};
245
246/** Plot matrix vs. fragment.
247 */
248bool CreateDataFragment(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateFragment)(class MatrixContainer &, int))
249{
250 stringstream filename;
251 ofstream output;
252
253 filename << prefix << ".dat";
254 if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
255 cout << msg << endl;
256 output << "# " << msg << ", created on " << datum << endl;
257 output << "#Order\tFrag.No.\t" << Fragment.Header << endl;
258 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
259 for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
260 output << BondOrder+1 << "\t" << KeySet.OrderSet[BondOrder][i]+1;
261 CreateFragment(Fragment, KeySet.OrderSet[BondOrder][i]);
262 for (int l=0;l<Fragment.ColumnCounter;l++)
263 output << scientific << "\t" << Fragment.Matrix[ KeySet.OrderSet[BondOrder][i] ][ Fragment.RowCounter[ KeySet.OrderSet[BondOrder][i] ] ][l];
264 output << endl;
265 }
266 }
267 output.close();
268 return true;
269};
270
271/** Copies fragment energy values into last matrix of \a Matrix with greatest total energy.
272 * \param &Matrix MatrixContainer with all fragment energy values
273 * \param &KeySet KeySetsContainer with associations of each fragment to a bond order
274 * \param BondOrder current bond order
275 */
276void CreateMaxFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySet, int BondOrder)
277{
278 for(int j=Fragments.RowCounter[ Fragments.MatrixCounter ];j--;) {
279 for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
280 if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
281 for (int k=Fragments.ColumnCounter;k--;)
282 Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
283 }
284 }
285 }
286};
287
288/** Copies fragment energy values into last matrix of \a Matrix with smallest total energy.
289 * \param &Matrix MatrixContainer with all fragment energy values
290 * \param &KeySet KeySetsContainer with associations of each fragment to a bond order
291 * \param BondOrder current bond order
292 */
293void CreateMinFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySet, int BondOrder)
294{
295 for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
296 int i=0;
297 do { // first get a minimum value unequal to 0
298 for (int k=Fragments.ColumnCounter;k--;)
299 Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
300 i++;
301 } while ((fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < MYEPSILON) && (i<KeySet.FragmentsPerOrder[BondOrder]));
302 for(;i<KeySet.FragmentsPerOrder[BondOrder];i++) { // then find lowest
303 if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) > fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
304 for (int k=Fragments.ColumnCounter;k--;)
305 Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
306 }
307 }
308 }
309};
310
311/** Plot matrix vs. fragment.
312 */
313bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int))
314{
315 stringstream filename;
316 ofstream output;
317
318 filename << prefix << ".dat";
319 if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
320 cout << msg << endl;
321 output << "# " << msg << ", created on " << datum;
322 output << "#Order\tFrag.No.\t" << Fragment.Header << endl;
323 // max
324 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
325 Fragment.SetLastMatrix(0.,0);
326 CreateFragmentOrder(Fragment, KeySet, BondOrder);
327 output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
328 for (int l=0;l<Fragment.ColumnCounter;l++)
329 output << scientific << "\t" << Fragment.Matrix[ Fragment.MatrixCounter ][ Fragment.RowCounter[ Fragment.MatrixCounter ]-1 ][l];
330 output << endl;
331 }
332 output.close();
333 return true;
334};
335
336/** Takes last but one row and copies into final row.
337 * \param Energy MatrixContainer with matrix values
338 * \param MatrixNumber the index for the ForceMatrix::matrix array
339 */
340void CreateEnergy(class MatrixContainer &Energy, int MatrixNumber)
341{
342 for(int k=0;k<Energy.ColumnCounter;k++)
343 Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber] ] [k] = Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber]-1 ] [k];
344};
345
346/** Scans forces for the minimum in magnitude.
347 * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
348 * \param Force ForceMatrix class containing matrix values
349 * \param MatrixNumber the index for the ForceMatrix::matrix array
350 */
351void CreateMinimumForce(class MatrixContainer &Force, int MatrixNumber)
352{
353 for (int l=Force.ColumnCounter;l--;)
354 Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
355 for (int l=5;l<Force.ColumnCounter;l+=3) {
356 double stored = 0;
357 int k=0;
358 do {
359 for (int m=NDIM;m--;) {
360 stored += Force.Matrix[MatrixNumber][ k ][l+m]
361 * Force.Matrix[MatrixNumber][ k ][l+m];
362 Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m] = Force.Matrix[MatrixNumber][ k ][l+m];
363 }
364 k++;
365 } while ((fabs(stored) < MYEPSILON) && (k<Force.RowCounter[MatrixNumber]));
366 for (;k<Force.RowCounter[MatrixNumber];k++) {
367 double tmp = 0;
368 for (int m=NDIM;m--;)
369 tmp += Force.Matrix[MatrixNumber][ k ][l+m]
370 * Force.Matrix[MatrixNumber][ k ][l+m];
371 if ((fabs(tmp) > MYEPSILON) && (tmp < stored)) { // current force is greater than stored
372 for (int m=NDIM;m--;)
373 Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m] = Force.Matrix[MatrixNumber][ k ][l+m];
374 stored = tmp;
375 }
376 }
377 }
378};
379
380/** Scans forces for the mean in magnitude.
381 * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
382 * \param Force ForceMatrix class containing matrix values
383 * \param MatrixNumber the index for the ForceMatrix::matrix array
384 */
385void CreateMeanForce(class MatrixContainer &Force, int MatrixNumber)
386{
387 int divisor = 0;
388 for (int l=Force.ColumnCounter;l--;)
389 Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
390 for (int l=5;l<Force.ColumnCounter;l+=3) {
391 double tmp = 0;
392 for (int k=Force.RowCounter[MatrixNumber];k--;) {
393 double norm = 0.;
394 for (int m=NDIM;m--;)
395 norm += Force.Matrix[MatrixNumber][ k ][l+m]
396 * Force.Matrix[MatrixNumber][ k ][l+m];
397 tmp += sqrt(norm);
398 if (fabs(norm) > MYEPSILON) divisor++;
399 }
400 tmp /= (double)divisor;
401 Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = tmp;
402 }
403};
404
405/** Scans forces for the maximum in magnitude.
406 * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
407 * \param Force ForceMatrix class containing matrix values
408 * \param MatrixNumber the index for the ForceMatrix::matrix array
409 */
410void CreateMaximumForce(class MatrixContainer &Force, int MatrixNumber)
411{
412 for (int l=5;l<Force.ColumnCounter;l+=3) {
413 double stored = 0;
414 for (int k=Force.RowCounter[MatrixNumber];k--;) {
415 double tmp = 0;
416 for (int m=NDIM;m--;)
417 tmp += Force.Matrix[MatrixNumber][ k ][l+m]
418 * Force.Matrix[MatrixNumber][ k ][l+m];
419 if (tmp > stored) { // current force is greater than stored
420 for (int m=NDIM;m--;)
421 Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m] = Force.Matrix[MatrixNumber][ k ][l+m];
422 stored = tmp;
423 }
424 }
425 }
426};
427
428/** Leaves the Force.Matrix as it is.
429 * \param Force ForceMatrix class containing matrix values
430 * \param MatrixNumber the index for the ForceMatrix::matrix array
431*/
432void CreateSameForce(class MatrixContainer &Force, int MatrixNumber)
433{
434 // does nothing
435};
436
437/** Adds vectorwise all forces.
438 * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
439 * \param Force ForceMatrix class containing matrix values
440 * \param MatrixNumber the index for the ForceMatrix::matrix array
441 */
442void CreateVectorSumForce(class MatrixContainer &Force, int MatrixNumber)
443{
444 for (int l=Force.ColumnCounter;l--;)
445 Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
446 for (int l=0;l<Force.ColumnCounter;l++) {
447 for (int k=Force.RowCounter[MatrixNumber];k--;)
448 Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] += Force.Matrix[MatrixNumber][k][l];
449 }
450};
451
452/** Writes the standard pyxplot header info.
453 * \param keycolumns number of columns of the key
454 * \param *key position of key
455 * \param *logscale axis for logscale
456 * \param *extraline extra set lines if desired
457 * \param mxtics small tics at ...
458 * \param xtics large tics at ...
459 * \param *xlabel label for x axis
460 * \param *ylabel label for y axis
461 */
462void CreatePlotHeader(ofstream &output, const char *prefix, const int keycolumns, const char *key, const char *logscale, const char *extraline, const int mxtics, const int xtics, const char *xlabel, const char *ylabel)
463{
464 //output << "#!/home/heber/build/pyxplot/pyxplot" << endl << endl;
465 output << "reset" << endl;
466 output << "set keycolumns "<< keycolumns << endl;
467 output << "set key " << key << endl;
468 output << "set mxtics "<< mxtics << endl;
469 output << "set xtics "<< xtics << endl;
470 if (logscale != NULL)
471 output << "set logscale " << logscale << endl;
472 if (extraline != NULL)
473 output << extraline << endl;
474 output << "set xlabel '" << xlabel << "'" << endl;
475 output << "set ylabel '" << ylabel << "'" << endl;
476 output << "set terminal eps color" << endl;
477 output << "set output '"<< prefix << ".eps'" << endl;
478};
479
480/** Creates the pyxplotfile for energy data.
481 * \param Matrix MatrixContainer with matrix values
482 * \param KeySet contains bond order
483 * \param *dir directory
484 * \param *prefix prefix for all filenames (without ending)
485 * \param keycolumns number of columns of the key
486 * \param *key position of key
487 * \param logscale axis for logscale
488 * \param mxtics small tics at ...
489 * \param xtics large tics at ...
490 * \param xlabel label for x axis
491 * \param xlabel label for x axis
492 * \param *xrange xrange
493 * \param *yrange yrange
494 * \param *xargument number of column to plot against (x axis)
495 * \param uses using string for plot command
496 * \param (*CreatePlotLines) function reference that writes a single plot line
497 * \return true if file was written successfully
498 */
499bool CreatePlotOrder(class MatrixContainer &Matrix, const class KeySetsContainer &KeySet, const char *dir, const char *prefix, const int keycolumns, const char *key, const char *logscale, const char *extraline, const int mxtics, const int xtics, const char *xlabel, const char *ylabel, const char *xrange, const char *yrange, const char *xargument, const char *uses, void (*CreatePlotLines)(ofstream &, class MatrixContainer &, const char *, const char *, const char *))
500{
501 stringstream filename;
502 ofstream output;
503
504 filename << prefix << ".pyx";
505 if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
506 CreatePlotHeader(output, prefix, keycolumns, key, logscale, extraline, mxtics, xtics, xlabel, ylabel);
507 output << "plot " << xrange << " " << yrange << " \\" << endl;
508 CreatePlotLines(output, Matrix, prefix, xargument, uses);
509 output.close();
510 return true;
511};
512
513/** Writes plot lines for absolute energies.
514 * \param output file handler
515 * \param Energy MatrixContainer with matrix values
516 * \param *prefix prefix of data file
517 * \param *xargument number of column to plot against (x axis)
518 * \param *uses uses command
519 */
520void AbsEnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses)
521{
522 char item[1024];
523 stringstream line(Energy.Header);
524
525 line >> item;
526 for (int i=2; i<= Energy.ColumnCounter;i++) {
527 line >> item;
528 output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":(abs($" << i+2 << ")) " << uses;
529 if (i != (Energy.ColumnCounter))
530 output << ", \\";
531 output << endl;
532 }
533};
534
535/** Writes plot lines for energies.
536 * \param output file handler
537 * \param Energy MatrixContainer with matrix values
538 * \param *prefix prefix of data file
539 * \param *xargument number of column to plot against (x axis)
540 * \param *uses uses command
541 */
542void EnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses)
543{
544 char item[1024];
545 stringstream line(Energy.Header);
546
547 line >> item;
548 for (int i=2; i<= Energy.ColumnCounter;i++) {
549 line >> item;
550 output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":" << i+2 << " " << uses;
551 if (i != (Energy.ColumnCounter))
552 output << ", \\";
553 output << endl;
554 }
555};
556
557/** Writes plot lines for absolute force magnitudes.
558 * \param output file handler
559 * \param Force MatrixContainer with matrix values
560 * \param *prefix prefix of data file
561 * \param *xargument number of column to plot against (x axis)
562 * \param *uses uses command
563 */
564void ForceMagnitudePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
565{
566 char item[1024];
567 stringstream line(Force.Header);
568
569 line >> item;
570 line >> item;
571 line >> item;
572 line >> item;
573 line >> item;
574 for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
575 line >> item;
576 item[strlen(item)-1] = '\0'; // kill residual index char (the '0')
577 output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":(sqrt($" << i+1 << "*$" << i+1 << "+$" << i+2 << "*$" << i+2 << "+$" << i+3 << "*$" << i+3 << ")) " << uses;
578 if (i != (Force.ColumnCounter-1))
579 output << ", \\";
580 output << endl;
581 line >> item;
582 line >> item;
583 }
584};
585
586/** Writes plot lines for first component of force vector.
587 * \param output file handler
588 * \param Force MatrixContainer with matrix values
589 * \param *prefix prefix of data file
590 * \param *xargument number of column to plot against (x axis)
591 * \param *uses uses command
592 */
593void AbsFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
594{
595 char item[1024];
596 stringstream line(Force.Header);
597
598 line >> item;
599 line >> item;
600 line >> item;
601 line >> item;
602 line >> item;
603 for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
604 line >> item;
605 item[strlen(item)-1] = '\0'; // kill residual index char (the '0')
606 output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":(abs($" << i+1 << ")) " << uses;
607 if (i != (Force.ColumnCounter-1))
608 output << ", \\";
609 output << endl;
610 line >> item;
611 line >> item;
612 }
613};
614
615/** Writes plot lines for force vector as boxes with a fillcolor.
616 * \param output file handler
617 * \param Force MatrixContainer with matrix values
618 * \param *prefix prefix of data file
619 * \param *xargument number of column to plot against (x axis)
620 * \param *uses uses command
621 */
622void BoxesForcePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
623{
624 char item[1024];
625 stringstream line(Force.Header);
626 char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
627
628 line >> item;
629 line >> item;
630 line >> item;
631 line >> item;
632 line >> item;
633 for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
634 line >> item;
635 item[strlen(item)-1] = '\0'; // kill residual index char (the '0')
636 output << "'" << prefix << ".dat' title '" << item << "' using ($" << xargument << "+" << fixed << setprecision(1) << (double)((i-7)/3)*0.2 << "):(sqrt($" << i+1 << "*$" << i+1 << "+$" << i+2 << "*$" << i+2 << "+$" << i+3 << "*$" << i+3 << ")) " << uses << " " << fillcolor[(i-7)/3];
637 if (i != (Force.ColumnCounter-1))
638 output << ", \\";
639 output << endl;
640 line >> item;
641 line >> item;
642 }
643};
644
645/** Writes plot lines for first force vector component as boxes with a fillcolor.
646 * \param output file handler
647 * \param Force MatrixContainer with matrix values
648 * \param *prefix prefix of data file
649 * \param *xargument number of column to plot against (x axis)
650 * \param *uses uses command
651 */
652void BoxesFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
653{
654 char item[1024];
655 stringstream line(Force.Header);
656 char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
657
658 line >> item;
659 line >> item;
660 line >> item;
661 line >> item;
662 line >> item;
663 for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
664 line >> item;
665 item[strlen(item)-1] = '\0'; // kill residual index char (the '0')
666 output << "'" << prefix << ".dat' title '" << item << "' using ($" << xargument << "+" << fixed << setprecision(1) << (double)((i-7)/3)*0.2 << "):" << i+1 << " " << uses << " " << fillcolor[(i-7)/3];
667 if (i != (Force.ColumnCounter-1))
668 output << ", \\";
669 output << endl;
670 line >> item;
671 line >> item;
672 }
673};
Note: See TracBrowser for help on using the repository browser.