source: src/builder.cpp@ 6ac7ee

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 6ac7ee was 6ac7ee, checked in by Frederik Heber <heber@…>, 17 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: 60.9 KB
Line 
1/** \file builder.cpp
2 *
3 * By stating absolute positions or binding angles and distances atomic positions of a molecule can be constructed.
4 * The output is the complete configuration file for PCP for direct use.
5 * Features:
6 * -# Atomic data is retrieved from a file, if not found requested and stored there for later re-use
7 * -# step-by-step construction of the molecule beginning either at a centre of with a certain atom
8 *
9 */
10
11/*! \mainpage Molecuilder - a molecular set builder
12 *
13 * This introductory shall briefly make aquainted with the program, helping in installing and a first run.
14 *
15 * \section about About the Program
16 *
17 * Molecuilder is a short program, written in C++, that enables the construction of a coordinate set for the
18 * atoms making up an molecule by the successive statement of binding angles and distances and referencing to
19 * already constructed atoms.
20 *
21 * A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello
22 * molecular dynamics implementation.
23 *
24 * \section install Installation
25 *
26 * Installation should without problems succeed as follows:
27 * -# ./configure (or: mkdir build;mkdir run;cd build; ../configure --bindir=../run)
28 * -# make
29 * -# make install
30 *
31 * Further useful commands are
32 * -# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n
33 * -# make doxygen-doc: Creates these html pages out of the documented source
34 *
35 * \section run Running
36 *
37 * The program can be executed by running: ./molecuilder
38 *
39 * Note, that it uses a database, called "elements.db", in the executable's directory. If the file is not found,
40 * it is created and any given data on elements of the periodic table will be stored therein and re-used on
41 * later re-execution.
42 *
43 * \section ref References
44 *
45 * For the special configuration file format, see the documentation of pcp.
46 *
47 */
48
49
50using namespace std;
51
52#include "boundary.hpp"
53#include "ellipsoid.hpp"
54#include "helpers.hpp"
55#include "molecules.hpp"
56
57/********************************************** Submenu routine **************************************/
58
59/** Submenu for adding atoms to the molecule.
60 * \param *periode periodentafel
61 * \param *mol the molecule to add to
62 */
63static void AddAtoms(periodentafel *periode, molecule *mol)
64{
65 atom *first, *second, *third, *fourth;
66 Vector **atoms;
67 Vector x,y,z,n; // coordinates for absolute point in cell volume
68 double a,b,c;
69 char choice; // menu choice char
70 bool valid;
71
72 cout << Verbose(0) << "===========ADD ATOM============================" << endl;
73 cout << Verbose(0) << " a - state absolute coordinates of atom" << endl;
74 cout << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
75 cout << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
76 cout << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
77 cout << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
78 cout << Verbose(0) << "all else - go back" << endl;
79 cout << Verbose(0) << "===============================================" << endl;
80 cout << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
81 cout << Verbose(0) << "INPUT: ";
82 cin >> choice;
83
84 switch (choice) {
85 case 'a': // absolute coordinates of atom
86 cout << Verbose(0) << "Enter absolute coordinates." << endl;
87 first = new atom;
88 first->x.AskPosition(mol->cell_size, false);
89 first->type = periode->AskElement(); // give type
90 mol->AddAtom(first); // add to molecule
91 break;
92
93 case 'b': // relative coordinates of atom wrt to reference point
94 first = new atom;
95 valid = true;
96 do {
97 if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
98 cout << Verbose(0) << "Enter reference coordinates." << endl;
99 x.AskPosition(mol->cell_size, true);
100 cout << Verbose(0) << "Enter relative coordinates." << endl;
101 first->x.AskPosition(mol->cell_size, false);
102 first->x.AddVector((const Vector *)&x);
103 cout << Verbose(0) << "\n";
104 } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
105 first->type = periode->AskElement(); // give type
106 mol->AddAtom(first); // add to molecule
107 break;
108
109 case 'c': // relative coordinates of atom wrt to already placed atom
110 first = new atom;
111 valid = true;
112 do {
113 if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
114 second = mol->AskAtom("Enter atom number: ");
115 cout << Verbose(0) << "Enter relative coordinates." << endl;
116 first->x.AskPosition(mol->cell_size, false);
117 for (int i=NDIM;i--;) {
118 first->x.x[i] += second->x.x[i];
119 }
120 } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
121 first->type = periode->AskElement(); // give type
122 mol->AddAtom(first); // add to molecule
123 break;
124
125 case 'd': // two atoms, two angles and a distance
126 first = new atom;
127 valid = true;
128 do {
129 if (!valid) {
130 cout << Verbose(0) << "Resulting coordinates out of cell - ";
131 first->x.Output((ofstream *)&cout);
132 cout << Verbose(0) << endl;
133 }
134 cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
135 second = mol->AskAtom("Enter central atom: ");
136 third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
137 fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
138 a = ask_value("Enter distance between central (first) and new atom: ");
139 b = ask_value("Enter angle between new, first and second atom (degrees): ");
140 b *= M_PI/180.;
141 bound(&b, 0., 2.*M_PI);
142 c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
143 c *= M_PI/180.;
144 bound(&c, -M_PI, M_PI);
145 cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
146/*
147 second->Output(1,1,(ofstream *)&cout);
148 third->Output(1,2,(ofstream *)&cout);
149 fourth->Output(1,3,(ofstream *)&cout);
150 n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
151 x.Copyvector(&second->x);
152 x.SubtractVector(&third->x);
153 x.Copyvector(&fourth->x);
154 x.SubtractVector(&third->x);
155
156 if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
157 cout << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
158 continue;
159 }
160 cout << Verbose(0) << "resulting relative coordinates: ";
161 z.Output((ofstream *)&cout);
162 cout << Verbose(0) << endl;
163 */
164 // calc axis vector
165 x.CopyVector(&second->x);
166 x.SubtractVector(&third->x);
167 x.Normalize();
168 cout << "x: ",
169 x.Output((ofstream *)&cout);
170 cout << endl;
171 z.MakeNormalVector(&second->x,&third->x,&fourth->x);
172 cout << "z: ",
173 z.Output((ofstream *)&cout);
174 cout << endl;
175 y.MakeNormalVector(&x,&z);
176 cout << "y: ",
177 y.Output((ofstream *)&cout);
178 cout << endl;
179
180 // rotate vector around first angle
181 first->x.CopyVector(&x);
182 first->x.RotateVector(&z,b - M_PI);
183 cout << "Rotated vector: ",
184 first->x.Output((ofstream *)&cout);
185 cout << endl;
186 // remove the projection onto the rotation plane of the second angle
187 n.CopyVector(&y);
188 n.Scale(first->x.Projection(&y));
189 cout << "N1: ",
190 n.Output((ofstream *)&cout);
191 cout << endl;
192 first->x.SubtractVector(&n);
193 cout << "Subtracted vector: ",
194 first->x.Output((ofstream *)&cout);
195 cout << endl;
196 n.CopyVector(&z);
197 n.Scale(first->x.Projection(&z));
198 cout << "N2: ",
199 n.Output((ofstream *)&cout);
200 cout << endl;
201 first->x.SubtractVector(&n);
202 cout << "2nd subtracted vector: ",
203 first->x.Output((ofstream *)&cout);
204 cout << endl;
205
206 // rotate another vector around second angle
207 n.CopyVector(&y);
208 n.RotateVector(&x,c - M_PI);
209 cout << "2nd Rotated vector: ",
210 n.Output((ofstream *)&cout);
211 cout << endl;
212
213 // add the two linear independent vectors
214 first->x.AddVector(&n);
215 first->x.Normalize();
216 first->x.Scale(a);
217 first->x.AddVector(&second->x);
218
219 cout << Verbose(0) << "resulting coordinates: ";
220 first->x.Output((ofstream *)&cout);
221 cout << Verbose(0) << endl;
222 } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
223 first->type = periode->AskElement(); // give type
224 mol->AddAtom(first); // add to molecule
225 break;
226
227 case 'e': // least square distance position to a set of atoms
228 first = new atom;
229 atoms = new (Vector*[128]);
230 valid = true;
231 for(int i=128;i--;)
232 atoms[i] = NULL;
233 int i=0, j=0;
234 cout << Verbose(0) << "Now we need at least three molecules.\n";
235 do {
236 cout << Verbose(0) << "Enter " << i+1 << "th atom: ";
237 cin >> j;
238 if (j != -1) {
239 second = mol->FindAtom(j);
240 atoms[i++] = &(second->x);
241 }
242 } while ((j != -1) && (i<128));
243 if (i >= 2) {
244 first->x.LSQdistance(atoms, i);
245
246 first->x.Output((ofstream *)&cout);
247 first->type = periode->AskElement(); // give type
248 mol->AddAtom(first); // add to molecule
249 } else {
250 delete first;
251 cout << Verbose(0) << "Please enter at least two vectors!\n";
252 }
253 break;
254 };
255};
256
257/** Submenu for centering the atoms in the molecule.
258 * \param *mol the molecule with all the atoms
259 */
260static void CenterAtoms(molecule *mol)
261{
262 Vector x, y, helper;
263 char choice; // menu choice char
264
265 cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
266 cout << Verbose(0) << " a - on origin" << endl;
267 cout << Verbose(0) << " b - on center of gravity" << endl;
268 cout << Verbose(0) << " c - within box with additional boundary" << endl;
269 cout << Verbose(0) << " d - within given simulation box" << endl;
270 cout << Verbose(0) << "all else - go back" << endl;
271 cout << Verbose(0) << "===============================================" << endl;
272 cout << Verbose(0) << "INPUT: ";
273 cin >> choice;
274
275 switch (choice) {
276 default:
277 cout << Verbose(0) << "Not a valid choice." << endl;
278 break;
279 case 'a':
280 cout << Verbose(0) << "Centering atoms in config file on origin." << endl;
281 mol->CenterOrigin((ofstream *)&cout, &x);
282 break;
283 case 'b':
284 cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
285 mol->CenterGravity((ofstream *)&cout, &x);
286 break;
287 case 'c':
288 cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
289 for (int i=0;i<NDIM;i++) {
290 cout << Verbose(0) << "Enter axis " << i << " boundary: ";
291 cin >> y.x[i];
292 }
293 mol->CenterEdge((ofstream *)&cout, &x); // make every coordinate positive
294 mol->Translate(&y); // translate by boundary
295 helper.CopyVector(&y);
296 helper.Scale(2.);
297 helper.AddVector(&x);
298 mol->SetBoxDimension(&helper); // update Box of atoms by boundary
299 break;
300 case 'd':
301 cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
302 for (int i=0;i<NDIM;i++) {
303 cout << Verbose(0) << "Enter axis " << i << " boundary: ";
304 cin >> x.x[i];
305 }
306 // center
307 mol->CenterInBox((ofstream *)&cout, &x);
308 // update Box of atoms by boundary
309 mol->SetBoxDimension(&x);
310 break;
311 }
312};
313
314/** Submenu for aligning the atoms in the molecule.
315 * \param *periode periodentafel
316 * \param *mol the molecule with all the atoms
317 */
318static void AlignAtoms(periodentafel *periode, molecule *mol)
319{
320 atom *first, *second, *third;
321 Vector x,n;
322 char choice; // menu choice char
323
324 cout << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
325 cout << Verbose(0) << " a - state three atoms defining align plane" << endl;
326 cout << Verbose(0) << " b - state alignment vector" << endl;
327 cout << Verbose(0) << " c - state two atoms in alignment direction" << endl;
328 cout << Verbose(0) << " d - align automatically by least square fit" << endl;
329 cout << Verbose(0) << "all else - go back" << endl;
330 cout << Verbose(0) << "===============================================" << endl;
331 cout << Verbose(0) << "INPUT: ";
332 cin >> choice;
333
334 switch (choice) {
335 default:
336 case 'a': // three atoms defining mirror plane
337 first = mol->AskAtom("Enter first atom: ");
338 second = mol->AskAtom("Enter second atom: ");
339 third = mol->AskAtom("Enter third atom: ");
340
341 n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
342 break;
343 case 'b': // normal vector of mirror plane
344 cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
345 n.AskPosition(mol->cell_size,0);
346 n.Normalize();
347 break;
348 case 'c': // three atoms defining mirror plane
349 first = mol->AskAtom("Enter first atom: ");
350 second = mol->AskAtom("Enter second atom: ");
351
352 n.CopyVector((const Vector *)&first->x);
353 n.SubtractVector((const Vector *)&second->x);
354 n.Normalize();
355 break;
356 case 'd':
357 char shorthand[4];
358 Vector a;
359 struct lsq_params param;
360 do {
361 fprintf(stdout, "Enter the element of atoms to be chosen: ");
362 fscanf(stdin, "%3s", shorthand);
363 } while ((param.type = periode->FindElement(shorthand)) == NULL);
364 cout << Verbose(0) << "Element is " << param.type->name << endl;
365 mol->GetAlignvector(&param);
366 for (int i=NDIM;i--;) {
367 x.x[i] = gsl_vector_get(param.x,i);
368 n.x[i] = gsl_vector_get(param.x,i+NDIM);
369 }
370 gsl_vector_free(param.x);
371 cout << Verbose(0) << "Offset vector: ";
372 x.Output((ofstream *)&cout);
373 cout << Verbose(0) << endl;
374 n.Normalize();
375 break;
376 };
377 cout << Verbose(0) << "Alignment vector: ";
378 n.Output((ofstream *)&cout);
379 cout << Verbose(0) << endl;
380 mol->Align(&n);
381};
382
383/** Submenu for mirroring the atoms in the molecule.
384 * \param *mol the molecule with all the atoms
385 */
386static void MirrorAtoms(molecule *mol)
387{
388 atom *first, *second, *third;
389 Vector n;
390 char choice; // menu choice char
391
392 cout << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;
393 cout << Verbose(0) << " a - state three atoms defining mirror plane" << endl;
394 cout << Verbose(0) << " b - state normal vector of mirror plane" << endl;
395 cout << Verbose(0) << " c - state two atoms in normal direction" << endl;
396 cout << Verbose(0) << "all else - go back" << endl;
397 cout << Verbose(0) << "===============================================" << endl;
398 cout << Verbose(0) << "INPUT: ";
399 cin >> choice;
400
401 switch (choice) {
402 default:
403 case 'a': // three atoms defining mirror plane
404 first = mol->AskAtom("Enter first atom: ");
405 second = mol->AskAtom("Enter second atom: ");
406 third = mol->AskAtom("Enter third atom: ");
407
408 n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
409 break;
410 case 'b': // normal vector of mirror plane
411 cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
412 n.AskPosition(mol->cell_size,0);
413 n.Normalize();
414 break;
415 case 'c': // three atoms defining mirror plane
416 first = mol->AskAtom("Enter first atom: ");
417 second = mol->AskAtom("Enter second atom: ");
418
419 n.CopyVector((const Vector *)&first->x);
420 n.SubtractVector((const Vector *)&second->x);
421 n.Normalize();
422 break;
423 };
424 cout << Verbose(0) << "Normal vector: ";
425 n.Output((ofstream *)&cout);
426 cout << Verbose(0) << endl;
427 mol->Mirror((const Vector *)&n);
428};
429
430/** Submenu for removing the atoms from the molecule.
431 * \param *mol the molecule with all the atoms
432 */
433static void RemoveAtoms(molecule *mol)
434{
435 atom *first, *second;
436 int axis;
437 double tmp1, tmp2;
438 char choice; // menu choice char
439
440 cout << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;
441 cout << Verbose(0) << " a - state atom for removal by number" << endl;
442 cout << Verbose(0) << " b - keep only in radius around atom" << endl;
443 cout << Verbose(0) << " c - remove this with one axis greater value" << endl;
444 cout << Verbose(0) << "all else - go back" << endl;
445 cout << Verbose(0) << "===============================================" << endl;
446 cout << Verbose(0) << "INPUT: ";
447 cin >> choice;
448
449 switch (choice) {
450 default:
451 case 'a':
452 if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: ")))
453 cout << Verbose(1) << "Atom removed." << endl;
454 else
455 cout << Verbose(1) << "Atom not found." << endl;
456 break;
457 case 'b':
458 second = mol->AskAtom("Enter number of atom as reference point: ");
459 cout << Verbose(0) << "Enter radius: ";
460 cin >> tmp1;
461 first = mol->start;
462 while(first->next != mol->end) {
463 first = first->next;
464 if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
465 mol->RemoveAtom(first);
466 }
467 break;
468 case 'c':
469 cout << Verbose(0) << "Which axis is it: ";
470 cin >> axis;
471 cout << Verbose(0) << "Left inward boundary: ";
472 cin >> tmp1;
473 cout << Verbose(0) << "Right inward boundary: ";
474 cin >> tmp2;
475 first = mol->start;
476 while(first->next != mol->end) {
477 first = first->next;
478 if ((first->x.x[axis] > tmp2) || (first->x.x[axis] < tmp1)) // out of boundary ...
479 mol->RemoveAtom(first);
480 }
481 break;
482 };
483 //mol->Output((ofstream *)&cout);
484 choice = 'r';
485};
486
487/** Submenu for measuring out the atoms in the molecule.
488 * \param *periode periodentafel
489 * \param *mol the molecule with all the atoms
490 */
491static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
492{
493 atom *first, *second, *third;
494 Vector x,y;
495 double min[256], tmp1, tmp2, tmp3;
496 int Z;
497 char choice; // menu choice char
498
499 cout << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;
500 cout << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;
501 cout << Verbose(0) << " b - calculate bond length between two atoms" << endl;
502 cout << Verbose(0) << " c - calculate bond angle" << endl;
503 cout << Verbose(0) << " d - calculate principal axis of the system" << endl;
504 cout << Verbose(0) << " e - calculate volume of the convex envelope" << endl;
505 cout << Verbose(0) << " f - calculate temperature from current velocity" << endl;
506 cout << Verbose(0) << " g - output all temperatures per step from velocities" << endl;
507 cout << Verbose(0) << "all else - go back" << endl;
508 cout << Verbose(0) << "===============================================" << endl;
509 cout << Verbose(0) << "INPUT: ";
510 cin >> choice;
511
512 switch(choice) {
513 default:
514 cout << Verbose(1) << "Not a valid choice." << endl;
515 break;
516 case 'a':
517 first = mol->AskAtom("Enter first atom: ");
518 for (int i=MAX_ELEMENTS;i--;)
519 min[i] = 0.;
520
521 second = mol->start;
522 while ((second->next != mol->end)) {
523 second = second->next; // advance
524 Z = second->type->Z;
525 tmp1 = 0.;
526 if (first != second) {
527 x.CopyVector((const Vector *)&first->x);
528 x.SubtractVector((const Vector *)&second->x);
529 tmp1 = x.Norm();
530 }
531 if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
532 //cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
533 }
534 for (int i=MAX_ELEMENTS;i--;)
535 if (min[i] != 0.) cout << Verbose(0) << "Minimum Bond length between " << first->type->name << " Atom " << first->nr << " and next Ion of type " << (periode->FindElement(i))->name << ": " << min[i] << " a.u." << endl;
536 break;
537
538 case 'b':
539 first = mol->AskAtom("Enter first atom: ");
540 second = mol->AskAtom("Enter second atom: ");
541 for (int i=NDIM;i--;)
542 min[i] = 0.;
543 x.CopyVector((const Vector *)&first->x);
544 x.SubtractVector((const Vector *)&second->x);
545 tmp1 = x.Norm();
546 cout << Verbose(1) << "Distance vector is ";
547 x.Output((ofstream *)&cout);
548 cout << "." << endl << "Norm of distance is " << tmp1 << "." << endl;
549 break;
550
551 case 'c':
552 cout << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl;
553 first = mol->AskAtom("Enter first atom: ");
554 second = mol->AskAtom("Enter central atom: ");
555 third = mol->AskAtom("Enter last atom: ");
556 tmp1 = tmp2 = tmp3 = 0.;
557 x.CopyVector((const Vector *)&first->x);
558 x.SubtractVector((const Vector *)&second->x);
559 y.CopyVector((const Vector *)&third->x);
560 y.SubtractVector((const Vector *)&second->x);
561 cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
562 cout << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
563 break;
564 case 'd':
565 cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
566 cout << Verbose(0) << "Shall we rotate? [0/1]: ";
567 cin >> Z;
568 if ((Z >=0) && (Z <=1))
569 mol->PrincipalAxisSystem((ofstream *)&cout, (bool)Z);
570 else
571 mol->PrincipalAxisSystem((ofstream *)&cout, false);
572 break;
573 case 'e':
574 cout << Verbose(0) << "Evaluating volume of the convex envelope.";
575 VolumeOfConvexEnvelope((ofstream *)&cout, NULL, configuration, NULL, mol);
576 break;
577 case 'f':
578 mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps, (ofstream *)&cout);
579 break;
580 case 'g':
581 {
582 char filename[255];
583 cout << "Please enter filename: " << endl;
584 cin >> filename;
585 cout << Verbose(1) << "Storing temperatures in " << filename << "." << endl;
586 ofstream *output = new ofstream(filename, ios::trunc);
587 if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output))
588 cout << Verbose(2) << "File could not be written." << endl;
589 else
590 cout << Verbose(2) << "File stored." << endl;
591 output->close();
592 delete(output);
593 }
594 break;
595 }
596};
597
598/** Submenu for measuring out the atoms in the molecule.
599 * \param *mol the molecule with all the atoms
600 * \param *configuration configuration structure for the to be written config files of all fragments
601 */
602static void FragmentAtoms(molecule *mol, config *configuration)
603{
604 int Order1;
605 clock_t start, end;
606
607 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
608 cout << Verbose(0) << "What's the desired bond order: ";
609 cin >> Order1;
610 if (mol->first->next != mol->last) { // there are bonds
611 start = clock();
612 mol->FragmentMolecule((ofstream *)&cout, Order1, configuration);
613 end = clock();
614 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
615 } else
616 cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
617};
618
619/********************************************** Test routine **************************************/
620
621/** Is called always as option 'T' in the menu.
622 */
623static void testroutine(molecule *mol)
624{
625 // the current test routine checks the functionality of the KeySet&Graph concept:
626 // We want to have a multiindex (the KeySet) describing a unique subgraph
627 atom *Walker = mol->start;
628 int i, comp, counter=0;
629
630 // generate some KeySets
631 cout << "Generating KeySets." << endl;
632 KeySet TestSets[mol->AtomCount+1];
633 i=1;
634 while (Walker->next != mol->end) {
635 Walker = Walker->next;
636 for (int j=0;j<i;j++) {
637 TestSets[j].insert(Walker->nr);
638 }
639 i++;
640 }
641 cout << "Testing insertion of already present item in KeySets." << endl;
642 KeySetTestPair test;
643 test = TestSets[mol->AtomCount-1].insert(Walker->nr);
644 if (test.second) {
645 cout << Verbose(1) << "Insertion worked?!" << endl;
646 } else {
647 cout << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
648 }
649 TestSets[mol->AtomCount].insert(mol->end->previous->nr);
650 TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
651
652 // constructing Graph structure
653 cout << "Generating Subgraph class." << endl;
654 Graph Subgraphs;
655
656 // insert KeySets into Subgraphs
657 cout << "Inserting KeySets into Subgraph class." << endl;
658 for (int j=0;j<mol->AtomCount;j++) {
659 Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
660 }
661 cout << "Testing insertion of already present item in Subgraph." << endl;
662 GraphTestPair test2;
663 test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.)));
664 if (test2.second) {
665 cout << Verbose(1) << "Insertion worked?!" << endl;
666 } else {
667 cout << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;
668 }
669
670 // show graphs
671 cout << "Showing Subgraph's contents, checking that it's sorted." << endl;
672 Graph::iterator A = Subgraphs.begin();
673 while (A != Subgraphs.end()) {
674 cout << (*A).second.first << ": ";
675 KeySet::iterator key = (*A).first.begin();
676 comp = -1;
677 while (key != (*A).first.end()) {
678 if ((*key) > comp)
679 cout << (*key) << " ";
680 else
681 cout << (*key) << "! ";
682 comp = (*key);
683 key++;
684 }
685 cout << endl;
686 A++;
687 }
688};
689
690/** Tries given filename or standard on saving the config file.
691 * \param *ConfigFileName name of file
692 * \param *configuration pointer to configuration structure with all the values
693 * \param *periode pointer to periodentafel structure with all the elements
694 * \param *mol pointer to molecule structure with all the atoms and coordinates
695 */
696static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, molecule *mol)
697{
698 char filename[MAXSTRINGSIZE];
699 ofstream output;
700
701 cout << Verbose(0) << "Storing configuration ... " << endl;
702 // get correct valence orbitals
703 mol->CalculateOrbitals(*configuration);
704 configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble;
705 strcpy(filename, ConfigFileName);
706 if (ConfigFileName != NULL) { // test the file name
707 output.open(ConfigFileName, ios::trunc);
708 } else if (strlen(configuration->configname) != 0) {
709 strcpy(filename, configuration->configname);
710 output.open(configuration->configname, ios::trunc);
711 } else {
712 strcpy(filename, DEFAULTCONFIG);
713 output.open(DEFAULTCONFIG, ios::trunc);
714 }
715 output.close();
716 output.clear();
717 cout << Verbose(0) << "Saving of config file ";
718 if (configuration->Save(filename, periode, mol))
719 cout << "successful." << endl;
720 else
721 cout << "failed." << endl;
722
723 // and save to xyz file
724 if (ConfigFileName != NULL) {
725 strcpy(filename, ConfigFileName);
726 strcat(filename, ".xyz");
727 output.open(filename, ios::trunc);
728 }
729 if (output == NULL) {
730 strcpy(filename,"main_pcp_linux");
731 strcat(filename, ".xyz");
732 output.open(filename, ios::trunc);
733 }
734 cout << Verbose(0) << "Saving of XYZ file ";
735 if (mol->MDSteps <= 1) {
736 if (mol->OutputXYZ(&output))
737 cout << "successful." << endl;
738 else
739 cout << "failed." << endl;
740 } else {
741 if (mol->OutputTrajectoriesXYZ(&output))
742 cout << "successful." << endl;
743 else
744 cout << "failed." << endl;
745 }
746 output.close();
747 output.clear();
748
749 // and save as MPQC configuration
750 if (ConfigFileName != NULL)
751 strcpy(filename, ConfigFileName);
752 if (output == NULL)
753 strcpy(filename,"main_pcp_linux");
754 cout << Verbose(0) << "Saving as mpqc input ";
755 if (configuration->SaveMPQC(filename, mol))
756 cout << "done." << endl;
757 else
758 cout << "failed." << endl;
759
760 if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
761 cerr << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl;
762 }
763};
764
765/** Parses the command line options.
766 * \param argc argument count
767 * \param **argv arguments array
768 * \param *mol molecule structure
769 * \param *periode elements structure
770 * \param configuration config file structure
771 * \param *ConfigFileName pointer to config file name in **argv
772 * \param *PathToDatabases pointer to db's path in **argv
773 * \return exit code (0 - successful, all else - something's wrong)
774 */
775static int ParseCommandLineOptions(int argc, char **argv, molecule *&mol, periodentafel *&periode, config& configuration, char *&ConfigFileName, char *&PathToDatabases)
776{
777 Vector x,y,z,n; // coordinates for absolute point in cell volume
778 double *factor; // unit factor if desired
779 ifstream test;
780 ofstream output;
781 string line;
782 atom *first;
783 bool SaveFlag = false;
784 int ExitFlag = 0;
785 int j;
786 double volume = 0.;
787 enum ConfigStatus config_present = absent;
788 clock_t start,end;
789 int argptr;
790 PathToDatabases = LocalPath;
791
792 if (argc > 1) { // config file specified as option
793 // 1. : Parse options that just set variables or print help
794 argptr = 1;
795 do {
796 if (argv[argptr][0] == '-') {
797 cout << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n";
798 argptr++;
799 switch(argv[argptr-1][1]) {
800 case 'h':
801 case 'H':
802 case '?':
803 cout << "MoleCuilder suite" << endl << "==================" << endl << endl;
804 cout << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl;
805 cout << "or simply " << argv[0] << " without arguments for interactive session." << endl;
806 cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;
807 cout << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl;
808 cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
809 cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
810 cout << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
811 cout << "\t-O\tCenter atoms in origin." << endl;
812 cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
813 cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
814 cout << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
815 cout << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
816 cout << "\t-h/-H/-?\tGive this help screen." << endl;
817 cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
818 cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
819 cout << "\t-N\tGet non-convex-envelope." << endl;
820 cout << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;
821 cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
822 cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
823 cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
824 cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
825 cout << "\t-T <file> Store temperatures from the config file in <file>." << endl;
826 cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
827 cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl;
828 cout << "\t-v/-V\t\tGives version information." << endl;
829 cout << "Note: config files must not begin with '-' !" << endl;
830 delete(mol);
831 delete(periode);
832 return (1);
833 break;
834 case 'v':
835 case 'V':
836 cout << argv[0] << " " << VERSIONSTRING << endl;
837 cout << "Build your own molecule position set." << endl;
838 delete(mol);
839 delete(periode);
840 return (1);
841 break;
842 case 'e':
843 if ((argptr >= argc) || (argv[argptr][0] == '-')) {
844 cerr << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl;
845 } else {
846 cout << "Using " << argv[argptr] << " as elements database." << endl;
847 PathToDatabases = argv[argptr];
848 argptr+=1;
849 }
850 break;
851 case 'n':
852 cout << "I won't parse trajectories." << endl;
853 configuration.FastParsing = true;
854 break;
855 default: // no match? Step on
856 argptr++;
857 break;
858 }
859 } else
860 argptr++;
861 } while (argptr < argc);
862
863 // 2. Parse the element database
864 if (periode->LoadPeriodentafel(PathToDatabases)) {
865 cout << Verbose(0) << "Element list loaded successfully." << endl;
866 //periode->Output((ofstream *)&cout);
867 } else {
868 cout << Verbose(0) << "Element list loading failed." << endl;
869 return 1;
870 }
871 // 3. Find config file name and parse if possible
872 if (argv[1][0] != '-') {
873 cout << Verbose(0) << "Config file given." << endl;
874 test.open(argv[1], ios::in);
875 if (test == NULL) {
876 //return (1);
877 output.open(argv[1], ios::out);
878 if (output == NULL) {
879 cout << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl;
880 config_present = absent;
881 } else {
882 cout << "Empty configuration file." << endl;
883 ConfigFileName = argv[1];
884 config_present = empty;
885 output.close();
886 }
887 } else {
888 test.close();
889 ConfigFileName = argv[1];
890 cout << Verbose(1) << "Specified config file found, parsing ... ";
891 switch (configuration.TestSyntax(ConfigFileName, periode, mol)) {
892 case 1:
893 cout << "new syntax." << endl;
894 configuration.Load(ConfigFileName, periode, mol);
895 config_present = present;
896 break;
897 case 0:
898 cout << "old syntax." << endl;
899 configuration.LoadOld(ConfigFileName, periode, mol);
900 config_present = present;
901 break;
902 default:
903 cout << "Unknown syntax or empty, yet present file." << endl;
904 config_present = empty;
905 }
906 }
907 } else
908 config_present = absent;
909 // 4. parse again through options, now for those depending on elements db and config presence
910 argptr = 1;
911 do {
912 cout << "Current Command line argument: " << argv[argptr] << "." << endl;
913 if (argv[argptr][0] == '-') {
914 argptr++;
915 if ((config_present == present) || (config_present == empty)) {
916 switch(argv[argptr-1][1]) {
917 case 'p':
918 ExitFlag = 1;
919 if ((argptr >= argc) || (argv[argptr][0] == '-')) {
920 ExitFlag = 255;
921 cerr << "Not enough arguments for parsing: -p <xyz file>" << endl;
922 } else {
923 SaveFlag = true;
924 cout << Verbose(1) << "Parsing xyz file for new atoms." << endl;
925 if (!mol->AddXYZFile(argv[argptr]))
926 cout << Verbose(2) << "File not found." << endl;
927 else {
928 cout << Verbose(2) << "File found and parsed." << endl;
929 config_present = present;
930 }
931 }
932 break;
933 case 'a':
934 ExitFlag = 1;
935 if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1]))) {
936 ExitFlag = 255;
937 cerr << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl;
938 } else {
939 SaveFlag = true;
940 cout << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";
941 first = new atom;
942 first->type = periode->FindElement(atoi(argv[argptr]));
943 if (first->type != NULL)
944 cout << Verbose(2) << "found element " << first->type->name << endl;
945 for (int i=NDIM;i--;)
946 first->x.x[i] = atof(argv[argptr+1+i]);
947 if (first->type != NULL) {
948 mol->AddAtom(first); // add to molecule
949 if ((config_present == empty) && (mol->AtomCount != 0))
950 config_present = present;
951 } else
952 cerr << Verbose(1) << "Could not find the specified element." << endl;
953 argptr+=4;
954 }
955 break;
956 default: // no match? Don't step on (this is done in next switch's default)
957 break;
958 }
959 }
960 if (config_present == present) {
961 switch(argv[argptr-1][1]) {
962 case 'D':
963 ExitFlag = 1;
964 {
965 cout << Verbose(1) << "Depth-First-Search Analysis." << endl;
966 MoleculeLeafClass *Subgraphs = NULL; // list of subgraphs from DFS analysis
967 int *MinimumRingSize = new int[mol->AtomCount];
968 atom ***ListOfLocalAtoms = NULL;
969 int FragmentCounter = 0;
970 class StackClass<bond *> *BackEdgeStack = NULL;
971 class StackClass<bond *> *LocalBackEdgeStack = NULL;
972 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr]), configuration.GetIsAngstroem());
973 mol->CreateListOfBondsPerAtom((ofstream *)&cout);
974 Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
975 if (Subgraphs != NULL) {
976 Subgraphs->next->FillBondStructureFromReference((ofstream *)&cout, mol, (FragmentCounter = 0), ListOfLocalAtoms, false); // we want to keep the created ListOfLocalAtoms
977 while (Subgraphs->next != NULL) {
978 Subgraphs = Subgraphs->next;
979 LocalBackEdgeStack = new StackClass<bond *> (Subgraphs->Leaf->BondCount);
980 Subgraphs->Leaf->PickLocalBackEdges((ofstream *)&cout, ListOfLocalAtoms[FragmentCounter++], BackEdgeStack, LocalBackEdgeStack);
981 Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize);
982 delete(LocalBackEdgeStack);
983 delete(Subgraphs->previous);
984 }
985 delete(Subgraphs);
986 for (int i=0;i<FragmentCounter;i++)
987 Free((void **)&ListOfLocalAtoms[FragmentCounter], "ParseCommandLineOptions: **ListOfLocalAtoms[]");
988 Free((void **)&ListOfLocalAtoms, "ParseCommandLineOptions: ***ListOfLocalAtoms");
989 }
990 delete(BackEdgeStack);
991 delete[](MinimumRingSize);
992 }
993 //argptr+=1;
994 break;
995 case 'E':
996 ExitFlag = 1;
997 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) {
998 ExitFlag = 255;
999 cerr << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl;
1000 } else {
1001 SaveFlag = true;
1002 cout << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl;
1003 first = mol->FindAtom(atoi(argv[argptr]));
1004 first->type = periode->FindElement(atoi(argv[argptr+1]));
1005 argptr+=2;
1006 }
1007 break;
1008 case 'A':
1009 ExitFlag = 1;
1010 if ((argptr >= argc) || (argv[argptr][0] == '-')) {
1011 ExitFlag =255;
1012 cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
1013 } else {
1014 cout << "Parsing bonds from " << argv[argptr] << "." << endl;
1015 ifstream *input = new ifstream(argv[argptr]);
1016 mol->CreateAdjacencyList2((ofstream *)&cout, input);
1017 input->close();
1018 argptr+=1;
1019 }
1020 break;
1021 case 'N':
1022 ExitFlag = 1;
1023 if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
1024 ExitFlag = 255;
1025 cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
1026 } else {
1027 class Tesselation T;
1028 int N = 15;
1029 int number = 100;
1030 string filename(argv[argptr+1]);
1031 filename.append(".csv");
1032 cout << Verbose(0) << "Evaluating non-convex envelope.";
1033 cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
1034 LinkedCell LCList(mol, atof(argv[argptr])); // \NOTE not twice the radius??
1035 Find_non_convex_border((ofstream *)&cout, mol, &T, &LCList, argv[argptr+1], atof(argv[argptr]));
1036 FindDistributionOfEllipsoids((ofstream *)&cout, &T, &LCList, N, number, filename.c_str());
1037 argptr+=2;
1038 }
1039 break;
1040 case 'T':
1041 ExitFlag = 1;
1042 if ((argptr >= argc) || (argv[argptr][0] == '-')) {
1043 ExitFlag = 255;
1044 cerr << "Not enough or invalid arguments given for storing tempature: -T <temperature file>" << endl;
1045 } else {
1046 cout << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;
1047 ofstream *output = new ofstream(argv[argptr], ios::trunc);
1048 if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output))
1049 cout << Verbose(2) << "File could not be written." << endl;
1050 else
1051 cout << Verbose(2) << "File stored." << endl;
1052 output->close();
1053 delete(output);
1054 argptr+=1;
1055 }
1056 break;
1057 case 'P':
1058 ExitFlag = 1;
1059 if ((argptr >= argc) || (argv[argptr][0] == '-')) {
1060 ExitFlag = 255;
1061 cerr << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl;
1062 } else {
1063 SaveFlag = true;
1064 cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
1065 if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem()))
1066 cout << Verbose(2) << "File not found." << endl;
1067 else
1068 cout << Verbose(2) << "File found and parsed." << endl;
1069 argptr+=1;
1070 }
1071 break;
1072 case 't':
1073 ExitFlag = 1;
1074 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
1075 ExitFlag = 255;
1076 cerr << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl;
1077 } else {
1078 ExitFlag = 1;
1079 SaveFlag = true;
1080 cout << Verbose(1) << "Translating all ions to new origin." << endl;
1081 for (int i=NDIM;i--;)
1082 x.x[i] = atof(argv[argptr+i]);
1083 mol->Translate((const Vector *)&x);
1084 argptr+=3;
1085 }
1086 break;
1087 case 's':
1088 ExitFlag = 1;
1089 if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
1090 ExitFlag = 255;
1091 cerr << "Not enough or invalid arguments given for scaling: -s <factor/[factor_x]> [factor_y] [factor_z]" << endl;
1092 } else {
1093 SaveFlag = true;
1094 j = -1;
1095 cout << Verbose(1) << "Scaling all ion positions by factor." << endl;
1096 factor = new double[NDIM];
1097 factor[0] = atof(argv[argptr]);
1098 if ((argptr < argc) && (IsValidNumber(argv[argptr])))
1099 argptr++;
1100 factor[1] = atof(argv[argptr]);
1101 if ((argptr < argc) && (IsValidNumber(argv[argptr])))
1102 argptr++;
1103 factor[2] = atof(argv[argptr]);
1104 mol->Scale(&factor);
1105 for (int i=0;i<NDIM;i++) {
1106 j += i+1;
1107 x.x[i] = atof(argv[NDIM+i]);
1108 mol->cell_size[j]*=factor[i];
1109 }
1110 delete[](factor);
1111 argptr+=1;
1112 }
1113 break;
1114 case 'b':
1115 ExitFlag = 1;
1116 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
1117 ExitFlag = 255;
1118 cerr << "Not enough or invalid arguments given for centering in box: -b <length_x> <length_y> <length_z>" << endl;
1119 } else {
1120 SaveFlag = true;
1121 j = -1;
1122 cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
1123 j=-1;
1124 for (int i=0;i<NDIM;i++) {
1125 j += i+1;
1126 x.x[i] = atof(argv[argptr++]);
1127 mol->cell_size[j] += x.x[i]*2.;
1128 }
1129 // center
1130 mol->CenterInBox((ofstream *)&cout, &x);
1131 // update Box of atoms by boundary
1132 mol->SetBoxDimension(&x);
1133 }
1134 break;
1135 case 'c':
1136 ExitFlag = 1;
1137 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
1138 ExitFlag = 255;
1139 cerr << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl;
1140 } else {
1141 SaveFlag = true;
1142 j = -1;
1143 cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
1144 // make every coordinate positive
1145 mol->CenterEdge((ofstream *)&cout, &x);
1146 // update Box of atoms by boundary
1147 mol->SetBoxDimension(&x);
1148 // translate each coordinate by boundary
1149 j=-1;
1150 for (int i=0;i<NDIM;i++) {
1151 j += i+1;
1152 x.x[i] = atof(argv[argptr++]);
1153 mol->cell_size[j] += x.x[i]*2.;
1154 }
1155 mol->Translate((const Vector *)&x);
1156 }
1157 break;
1158 case 'O':
1159 ExitFlag = 1;
1160 SaveFlag = true;
1161 cout << Verbose(1) << "Centering atoms in origin." << endl;
1162 mol->CenterOrigin((ofstream *)&cout, &x);
1163 mol->SetBoxDimension(&x);
1164 break;
1165 case 'r':
1166 ExitFlag = 1;
1167 SaveFlag = true;
1168 cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
1169 break;
1170 case 'F':
1171 case 'f':
1172 ExitFlag = 1;
1173 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {
1174 ExitFlag = 255;
1175 cerr << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl;
1176 } else {
1177 cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;
1178 cout << Verbose(0) << "Creating connection matrix..." << endl;
1179 start = clock();
1180 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr++]), configuration.GetIsAngstroem());
1181 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
1182 if (mol->first->next != mol->last) {
1183 ExitFlag = mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);
1184 }
1185 end = clock();
1186 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
1187 argptr+=2;
1188 }
1189 break;
1190 case 'm':
1191 ExitFlag = 1;
1192 j = atoi(argv[argptr++]);
1193 if ((j<0) || (j>1)) {
1194 cerr << Verbose(1) << "ERROR: Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;
1195 j = 0;
1196 }
1197 if (j) {
1198 SaveFlag = true;
1199 cout << Verbose(0) << "Converting to prinicipal axis system." << endl;
1200 } else
1201 cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
1202 mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j);
1203 break;
1204 case 'o':
1205 ExitFlag = 1;
1206 if ((argptr >= argc) || (argv[argptr][0] == '-')){
1207 ExitFlag = 255;
1208 cerr << "Not enough or invalid arguments given for convex envelope: -o <tecplot output file>" << endl;
1209 } else {
1210 cout << Verbose(0) << "Evaluating volume of the convex envelope.";
1211 cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
1212 VolumeOfConvexEnvelope((ofstream *)&cout, argv[argptr], &configuration, NULL, mol);
1213 argptr+=1;
1214 }
1215 break;
1216 case 'U':
1217 ExitFlag = 1;
1218 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
1219 ExitFlag = 255;
1220 cerr << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl;
1221 volume = -1; // for case 'u': don't print error again
1222 } else {
1223 volume = atof(argv[argptr++]);
1224 cout << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;
1225 }
1226 case 'u':
1227 ExitFlag = 1;
1228 if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
1229 if (volume != -1)
1230 ExitFlag = 255;
1231 cerr << "Not enough arguments given for suspension: -u <density>" << endl;
1232 } else {
1233 double density;
1234 SaveFlag = true;
1235 cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
1236 density = atof(argv[argptr++]);
1237 if (density < 1.0) {
1238 cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
1239 density = 1.3;
1240 }
1241// for(int i=0;i<NDIM;i++) {
1242// repetition[i] = atoi(argv[argptr++]);
1243// if (repetition[i] < 1)
1244// cerr << Verbose(0) << "ERROR: repetition value must be greater 1!" << endl;
1245// repetition[i] = 1;
1246// }
1247 PrepareClustersinWater((ofstream *)&cout, &configuration, mol, volume, density); // if volume == 0, will calculate from ConvexEnvelope
1248 }
1249 break;
1250 case 'd':
1251 ExitFlag = 1;
1252 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
1253 ExitFlag = 255;
1254 cerr << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl;
1255 } else {
1256 SaveFlag = true;
1257 for (int axis = 1; axis <= NDIM; axis++) {
1258 int faktor = atoi(argv[argptr++]);
1259 int count;
1260 element ** Elements;
1261 Vector ** vectors;
1262 if (faktor < 1) {
1263 cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;
1264 faktor = 1;
1265 }
1266 mol->CountAtoms((ofstream *)&cout); // recount atoms
1267 if (mol->AtomCount != 0) { // if there is more than none
1268 count = mol->AtomCount; // is changed becausing of adding, thus has to be stored away beforehand
1269 Elements = new element *[count];
1270 vectors = new Vector *[count];
1271 j = 0;
1272 first = mol->start;
1273 while (first->next != mol->end) { // make a list of all atoms with coordinates and element
1274 first = first->next;
1275 Elements[j] = first->type;
1276 vectors[j] = &first->x;
1277 j++;
1278 }
1279 if (count != j)
1280 cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
1281 x.Zero();
1282 y.Zero();
1283 y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
1284 for (int i=1;i<faktor;i++) { // then add this list with respective translation factor times
1285 x.AddVector(&y); // per factor one cell width further
1286 for (int k=count;k--;) { // go through every atom of the original cell
1287 first = new atom(); // create a new body
1288 first->x.CopyVector(vectors[k]); // use coordinate of original atom
1289 first->x.AddVector(&x); // translate the coordinates
1290 first->type = Elements[k]; // insert original element
1291 mol->AddAtom(first); // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
1292 }
1293 }
1294 // free memory
1295 delete[](Elements);
1296 delete[](vectors);
1297 // correct cell size
1298 if (axis < 0) { // if sign was negative, we have to translate everything
1299 x.Zero();
1300 x.AddVector(&y);
1301 x.Scale(-(faktor-1));
1302 mol->Translate(&x);
1303 }
1304 mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
1305 }
1306 }
1307 }
1308 break;
1309 default: // no match? Step on
1310 if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step!
1311 argptr++;
1312 break;
1313 }
1314 }
1315 } else argptr++;
1316 } while (argptr < argc);
1317 if (SaveFlag)
1318 SaveConfig(ConfigFileName, &configuration, periode, mol);
1319 if ((ExitFlag >= 1)) {
1320 delete(mol);
1321 delete(periode);
1322 return (ExitFlag);
1323 }
1324 } else { // no arguments, hence scan the elements db
1325 if (periode->LoadPeriodentafel(PathToDatabases))
1326 cout << Verbose(0) << "Element list loaded successfully." << endl;
1327 else
1328 cout << Verbose(0) << "Element list loading failed." << endl;
1329 configuration.RetrieveConfigPathAndName("main_pcp_linux");
1330 }
1331 return(0);
1332};
1333
1334/********************************************** Main routine **************************************/
1335
1336int main(int argc, char **argv)
1337{
1338 periodentafel *periode = new periodentafel; // and a period table of all elements
1339 molecule *mol = new molecule(periode); // first we need an empty molecule
1340 config configuration;
1341 double tmp1;
1342 double bond, min_bond;
1343 atom *first, *second;
1344 char choice; // menu choice char
1345 Vector x,y,z,n; // coordinates for absolute point in cell volume
1346 double *factor; // unit factor if desired
1347 bool valid; // flag if input was valid or not
1348 ifstream test;
1349 ofstream output;
1350 string line;
1351 char filename[MAXSTRINGSIZE];
1352 char *ConfigFileName = NULL;
1353 char *ElementsFileName = NULL;
1354 int Z;
1355 int j, axis, count, faktor;
1356 clock_t start,end;
1357// int *MinimumRingSize = NULL;
1358 MoleculeLeafClass *Subgraphs = NULL;
1359// class StackClass<bond *> *BackEdgeStack = NULL;
1360 element **Elements;
1361 Vector **vectors;
1362
1363 // =========================== PARSE COMMAND LINE OPTIONS ====================================
1364 j = ParseCommandLineOptions(argc, argv, mol, periode, configuration, ConfigFileName, ElementsFileName);
1365 if (j == 1) return 0; // just for -v and -h options
1366 if (j) return j; // something went wrong
1367
1368 // General stuff
1369 if (mol->cell_size[0] == 0.) {
1370 cout << Verbose(0) << "enter lower triadiagonal form of basis matrix" << endl << endl;
1371 for (int i=0;i<6;i++) {
1372 cout << Verbose(1) << "Cell size" << i << ": ";
1373 cin >> mol->cell_size[i];
1374 }
1375 }
1376
1377 // =========================== START INTERACTIVE SESSION ====================================
1378
1379 // now the main construction loop
1380 cout << Verbose(0) << endl << "Now comes the real construction..." << endl;
1381 do {
1382 cout << Verbose(0) << endl << endl;
1383 cout << Verbose(0) << "============Element list=======================" << endl;
1384 mol->Checkout((ofstream *)&cout);
1385 cout << Verbose(0) << "============Atom list==========================" << endl;
1386 mol->Output((ofstream *)&cout);
1387 cout << Verbose(0) << "============Menu===============================" << endl;
1388 cout << Verbose(0) << "a - add an atom" << endl;
1389 cout << Verbose(0) << "r - remove an atom" << endl;
1390 cout << Verbose(0) << "b - scale a bond between atoms" << endl;
1391 cout << Verbose(0) << "u - change an atoms element" << endl;
1392 cout << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
1393 cout << Verbose(0) << "-----------------------------------------------" << endl;
1394 cout << Verbose(0) << "p - Parse xyz file" << endl;
1395 cout << Verbose(0) << "e - edit the current configuration" << endl;
1396 cout << Verbose(0) << "o - create connection matrix" << endl;
1397 cout << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;
1398 cout << Verbose(0) << "-----------------------------------------------" << endl;
1399 cout << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
1400 cout << Verbose(0) << "i - realign molecule" << endl;
1401 cout << Verbose(0) << "m - mirror all molecules" << endl;
1402 cout << Verbose(0) << "t - translate molecule by vector" << endl;
1403 cout << Verbose(0) << "c - scale by unit transformation" << endl;
1404 cout << Verbose(0) << "g - center atoms in box" << endl;
1405 cout << Verbose(0) << "-----------------------------------------------" << endl;
1406 cout << Verbose(0) << "s - save current setup to config file" << endl;
1407 cout << Verbose(0) << "T - call the current test routine" << endl;
1408 cout << Verbose(0) << "q - quit" << endl;
1409 cout << Verbose(0) << "===============================================" << endl;
1410 cout << Verbose(0) << "Input: ";
1411 cin >> choice;
1412
1413 switch (choice) {
1414 default:
1415 case 'a': // add atom
1416 AddAtoms(periode, mol);
1417 choice = 'a';
1418 break;
1419
1420 case 'b': // scale a bond
1421 cout << Verbose(0) << "Scaling bond length between two atoms." << endl;
1422 first = mol->AskAtom("Enter first (fixed) atom: ");
1423 second = mol->AskAtom("Enter second (shifting) atom: ");
1424 min_bond = 0.;
1425 for (int i=NDIM;i--;)
1426 min_bond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
1427 min_bond = sqrt(min_bond);
1428 cout << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << min_bond << " a.u." << endl;
1429 cout << Verbose(0) << "Enter new bond length [a.u.]: ";
1430 cin >> bond;
1431 for (int i=NDIM;i--;) {
1432 second->x.x[i] -= (second->x.x[i]-first->x.x[i])/min_bond*(min_bond-bond);
1433 }
1434 //cout << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
1435 //second->Output(second->type->No, 1, (ofstream *)&cout);
1436 break;
1437
1438 case 'c': // unit scaling of the metric
1439 cout << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;
1440 cout << Verbose(0) << "Enter three factors: ";
1441 factor = new double[NDIM];
1442 cin >> factor[0];
1443 cin >> factor[1];
1444 cin >> factor[2];
1445 valid = true;
1446 mol->Scale(&factor);
1447 delete[](factor);
1448 break;
1449
1450 case 'd': // duplicate the periodic cell along a given axis, given times
1451 cout << Verbose(0) << "State the axis [(+-)123]: ";
1452 cin >> axis;
1453 cout << Verbose(0) << "State the factor: ";
1454 cin >> faktor;
1455
1456 mol->CountAtoms((ofstream *)&cout); // recount atoms
1457 if (mol->AtomCount != 0) { // if there is more than none
1458 count = mol->AtomCount; // is changed becausing of adding, thus has to be stored away beforehand
1459 Elements = new element *[count];
1460 vectors = new Vector *[count];
1461 j = 0;
1462 first = mol->start;
1463 while (first->next != mol->end) { // make a list of all atoms with coordinates and element
1464 first = first->next;
1465 Elements[j] = first->type;
1466 vectors[j] = &first->x;
1467 j++;
1468 }
1469 if (count != j)
1470 cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
1471 x.Zero();
1472 y.Zero();
1473 y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
1474 for (int i=1;i<faktor;i++) { // then add this list with respective translation factor times
1475 x.AddVector(&y); // per factor one cell width further
1476 for (int k=count;k--;) { // go through every atom of the original cell
1477 first = new atom(); // create a new body
1478 first->x.CopyVector(vectors[k]); // use coordinate of original atom
1479 first->x.AddVector(&x); // translate the coordinates
1480 first->type = Elements[k]; // insert original element
1481 mol->AddAtom(first); // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
1482 }
1483 }
1484 if (mol->first->next != mol->last) // if connect matrix is present already, redo it
1485 mol->CreateAdjacencyList((ofstream *)&cout, mol->BondDistance, configuration.GetIsAngstroem());
1486 // free memory
1487 delete[](Elements);
1488 delete[](vectors);
1489 // correct cell size
1490 if (axis < 0) { // if sign was negative, we have to translate everything
1491 x.Zero();
1492 x.AddVector(&y);
1493 x.Scale(-(faktor-1));
1494 mol->Translate(&x);
1495 }
1496 mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
1497 }
1498 break;
1499
1500 case 'e': // edit each field of the configuration
1501 configuration.Edit(mol);
1502 break;
1503
1504 case 'f':
1505 FragmentAtoms(mol, &configuration);
1506 break;
1507
1508 case 'g': // center the atoms
1509 CenterAtoms(mol);
1510 break;
1511
1512 case 'i': // align all atoms
1513 AlignAtoms(periode, mol);
1514 break;
1515
1516 case 'l': // measure distances or angles
1517 MeasureAtoms(periode, mol, &configuration);
1518 break;
1519
1520 case 'm': // mirror atoms along a given axis
1521 MirrorAtoms(mol);
1522 break;
1523
1524 case 'o': // create the connection matrix
1525 {
1526 cout << Verbose(0) << "What's the maximum bond distance: ";
1527 cin >> tmp1;
1528 start = clock();
1529 mol->CreateAdjacencyList((ofstream *)&cout, tmp1, configuration.GetIsAngstroem());
1530 mol->CreateListOfBondsPerAtom((ofstream *)&cout);
1531// Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
1532// while (Subgraphs->next != NULL) {
1533// Subgraphs = Subgraphs->next;
1534// Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize);
1535// delete(Subgraphs->previous);
1536// }
1537// delete(Subgraphs); // we don't need the list here, so free everything
1538// delete[](MinimumRingSize);
1539// Subgraphs = NULL;
1540 end = clock();
1541 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
1542 }
1543 break;
1544
1545 case 'p': // parse and XYZ file
1546 cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
1547 do {
1548 cout << Verbose(0) << "Enter file name: ";
1549 cin >> filename;
1550 } while (!mol->AddXYZFile(filename));
1551 break;
1552
1553 case 'q': // quit
1554 break;
1555
1556 case 'r': // remove atom
1557 RemoveAtoms(mol);
1558 break;
1559
1560 case 's': // save to config file
1561 SaveConfig(ConfigFileName, &configuration, periode, mol);
1562 break;
1563
1564 case 't': // translate all atoms
1565 cout << Verbose(0) << "Enter translation vector." << endl;
1566 x.AskPosition(mol->cell_size,0);
1567 mol->Translate((const Vector *)&x);
1568 break;
1569
1570 case 'T':
1571 testroutine(mol);
1572 break;
1573
1574 case 'u': // change an atom's element
1575 first = NULL;
1576 do {
1577 cout << Verbose(0) << "Change the element of which atom: ";
1578 cin >> Z;
1579 } while ((first = mol->FindAtom(Z)) == NULL);
1580 cout << Verbose(0) << "New element by atomic number Z: ";
1581 cin >> Z;
1582 first->type = periode->FindElement(Z);
1583 cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
1584 break;
1585 };
1586 } while (choice != 'q');
1587
1588 // save element data base
1589 if (periode->StorePeriodentafel(ElementsFileName)) //ElementsFileName
1590 cout << Verbose(0) << "Saving of elements.db successful." << endl;
1591 else
1592 cout << Verbose(0) << "Saving of elements.db failed." << endl;
1593
1594 // Free all
1595 if (Subgraphs != NULL) { // free disconnected subgraph list of DFS analysis was performed
1596 while (Subgraphs->next != NULL) {
1597 Subgraphs = Subgraphs->next;
1598 delete(Subgraphs->previous);
1599 }
1600 delete(Subgraphs);
1601 }
1602 delete(mol);
1603 delete(periode);
1604 return (0);
1605}
1606
1607/********************************************** E N D **************************************************/
Note: See TracBrowser for help on using the repository browser.