source: src/datacreator.cpp@ 02da9e

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 02da9e was 6ac7ee, checked in by Frederik Heber <heber@…>, 16 years ago

Merge branch 'ConcaveHull' of ../espack2 into ConcaveHull

Conflicts:

molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/linkedcell.cpp
molecuilder/src/linkedcell.hpp
molecuilder/src/vector.cpp
molecuilder/src/vector.hpp
util/src/NanoCreator.c

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