source: src/builder.cpp@ 14de469

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 14de469 was 14de469, checked in by Frederik Heber <heber@…>, 17 years ago

-initial commit
-Minimum set of files needed from ESPACK SVN repository
-Switch to three tantamount package parts instead of all relating to pcp (as at some time Ralf's might find inclusion as well)

  • Property mode set to 100644
File size: 43.6 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 "helpers.hpp"
53#include "molecules.hpp"
54
55/********************************************** Submenu routine **************************************/
56
57/** Submenu for adding atoms to the molecule.
58 * \param *periode periodentafel
59 * \param *mol the molecule to add to
60 */
61void AddAtoms(periodentafel *periode, molecule *mol)
62{
63 atom *first, *second, *third, *fourth;
64 vector **atoms;
65 vector x,y,z,n; // coordinates for absolute point in cell volume
66 double a,b,c;
67 char choice; // menu choice char
68 bool valid;
69
70 cout << Verbose(0) << "===========ADD ATOM============================" << endl;
71 cout << Verbose(0) << " a - state absolute coordinates of atom" << endl;
72 cout << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
73 cout << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
74 cout << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
75 cout << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
76 cout << Verbose(0) << "all else - go back" << endl;
77 cout << Verbose(0) << "===============================================" << endl;
78 cout << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
79 cout << Verbose(0) << "INPUT: ";
80 cin >> choice;
81
82 switch (choice) {
83 case 'a': // absolute coordinates of atom
84 cout << Verbose(0) << "Enter absolute coordinates." << endl;
85 first = new atom;
86 first->x.AskPosition(mol->cell_size, false);
87 first->type = periode->AskElement(); // give type
88 mol->AddAtom(first); // add to molecule
89 break;
90
91 case 'b': // relative coordinates of atom wrt to reference point
92 first = new atom;
93 valid = true;
94 do {
95 if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
96 cout << Verbose(0) << "Enter reference coordinates." << endl;
97 x.AskPosition(mol->cell_size, true);
98 cout << Verbose(0) << "Enter relative coordinates." << endl;
99 first->x.AskPosition(mol->cell_size, false);
100 first->x.AddVector((const vector *)&x);
101 cout << Verbose(0) << "\n";
102 } while (!(valid = mol->CheckBounds((const vector *)&first->x)));
103 first->type = periode->AskElement(); // give type
104 mol->AddAtom(first); // add to molecule
105 break;
106
107 case 'c': // relative coordinates of atom wrt to already placed atom
108 first = new atom;
109 valid = true;
110 do {
111 if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
112 second = mol->AskAtom("Enter atom number: ");
113 cout << Verbose(0) << "Enter relative coordinates." << endl;
114 first->x.AskPosition(mol->cell_size, false);
115 for (int i=0;i<3;i++) {
116 first->x.x[i] += second->x.x[i];
117 }
118 } while (!(valid = mol->CheckBounds((const vector *)&first->x)));
119 first->type = periode->AskElement(); // give type
120 mol->AddAtom(first); // add to molecule
121 break;
122
123 case 'd': // two atoms, two angles and a distance
124 first = new atom;
125 valid = true;
126 do {
127 if (!valid) {
128 cout << Verbose(0) << "Resulting coordinates out of cell - ";
129 first->x.Output((ofstream *)&cout);
130 cout << Verbose(0) << endl;
131 }
132 cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
133 second = mol->AskAtom("Enter central atom: ");
134 third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
135 fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
136 a = ask_value("Enter distance between central (first) and new atom: ");
137 b = ask_value("Enter angle between new, first and second atom (degrees): ");
138 b *= M_PI/180.;
139 bound(&b, 0., 2.*M_PI);
140 c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
141 c *= M_PI/180.;
142 bound(&c, -M_PI, M_PI);
143 cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
144/*
145 second->Output(1,1,(ofstream *)&cout);
146 third->Output(1,2,(ofstream *)&cout);
147 fourth->Output(1,3,(ofstream *)&cout);
148 n.MakeNormalVector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
149 x.CopyVector(&second->x);
150 x.SubtractVector(&third->x);
151 x.CopyVector(&fourth->x);
152 x.SubtractVector(&third->x);
153
154 if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
155 cout << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
156 continue;
157 }
158 cout << Verbose(0) << "resulting relative coordinates: ";
159 z.Output((ofstream *)&cout);
160 cout << Verbose(0) << endl;
161 */
162 // calc axis vector
163 x.CopyVector(&second->x);
164 x.SubtractVector(&third->x);
165 x.Normalize();
166 cout << "x: ",
167 x.Output((ofstream *)&cout);
168 cout << endl;
169 z.MakeNormalVector(&second->x,&third->x,&fourth->x);
170 cout << "z: ",
171 z.Output((ofstream *)&cout);
172 cout << endl;
173 y.MakeNormalVector(&x,&z);
174 cout << "y: ",
175 y.Output((ofstream *)&cout);
176 cout << endl;
177
178 // rotate vector around first angle
179 first->x.CopyVector(&x);
180 first->x.RotateVector(&z,b - M_PI);
181 cout << "Rotated vector: ",
182 first->x.Output((ofstream *)&cout);
183 cout << endl;
184 // remove the projection onto the rotation plane of the second angle
185 n.CopyVector(&y);
186 n.Scale(first->x.Projection(&y));
187 cout << "N1: ",
188 n.Output((ofstream *)&cout);
189 cout << endl;
190 first->x.SubtractVector(&n);
191 cout << "Subtracted vector: ",
192 first->x.Output((ofstream *)&cout);
193 cout << endl;
194 n.CopyVector(&z);
195 n.Scale(first->x.Projection(&z));
196 cout << "N2: ",
197 n.Output((ofstream *)&cout);
198 cout << endl;
199 first->x.SubtractVector(&n);
200 cout << "2nd subtracted vector: ",
201 first->x.Output((ofstream *)&cout);
202 cout << endl;
203
204 // rotate another vector around second angle
205 n.CopyVector(&y);
206 n.RotateVector(&x,c - M_PI);
207 cout << "2nd Rotated vector: ",
208 n.Output((ofstream *)&cout);
209 cout << endl;
210
211 // add the two linear independent vectors
212 first->x.AddVector(&n);
213 first->x.Normalize();
214 first->x.Scale(a);
215 first->x.AddVector(&second->x);
216
217 cout << Verbose(0) << "resulting coordinates: ";
218 first->x.Output((ofstream *)&cout);
219 cout << Verbose(0) << endl;
220 } while (!(valid = mol->CheckBounds((const vector *)&first->x)));
221 first->type = periode->AskElement(); // give type
222 mol->AddAtom(first); // add to molecule
223 break;
224
225 case 'e': // least square distance position to a set of atoms
226 first = new atom;
227 atoms = new (vector*[128]);
228 valid = true;
229 for(int i=0;i<128;i++)
230 atoms[i] = NULL;
231 int i=0, j=0;
232 cout << Verbose(0) << "Now we need at least three molecules.\n";
233 do {
234 cout << Verbose(0) << "Enter " << i+1 << "th atom: ";
235 cin >> j;
236 if (j != -1) {
237 second = mol->FindAtom(j);
238 atoms[i++] = &(second->x);
239 }
240 } while ((j != -1) && (i<128));
241 if (i >= 2) {
242 first->x.LSQdistance(atoms, i);
243
244 first->x.Output((ofstream *)&cout);
245 first->type = periode->AskElement(); // give type
246 mol->AddAtom(first); // add to molecule
247 } else {
248 delete first;
249 cout << Verbose(0) << "Please enter at least two vectors!\n";
250 }
251 break;
252 };
253};
254
255/** Submenu for centering the atoms in the molecule.
256 * \param *mol the molecule with all the atoms
257 */
258void CenterAtoms(molecule *mol)
259{
260 vector x;
261 char choice; // menu choice char
262 double min[3];
263 int j;
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) << "all else - go back" << endl;
270 cout << Verbose(0) << "===============================================" << endl;
271 cout << Verbose(0) << "INPUT: ";
272 cin >> choice;
273
274 switch (choice) {
275 default:
276 cout << Verbose(0) << "Not a valid choice." << endl;
277 break;
278 case 'a':
279 cout << Verbose(0) << "Centering atoms in config file on origin." << endl;
280 mol->CenterOrigin((ofstream *)&cout, &x);
281 break;
282 case 'b':
283 cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
284 mol->CenterGravity((ofstream *)&cout, &x);
285 break;
286 case 'c':
287 cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
288 for (int i=0;i<3;i++) {
289 cout << Verbose(0) << "Enter axis " << i << " boundary: ";
290 cin >> min[i];
291 }
292 mol->CenterEdge((ofstream *)&cout, &x); // make every coordinate positive
293 mol->SetBoxDimension(&x); // update Box of atoms by boundary
294 // translate each coordinate by boundary
295 j=-1;
296 for (int i=0;i<NDIM;i++) {
297 j += i+1;
298 x.x[i] = min[i];
299 mol->cell_size[j] += x.x[i]*2.;
300 }
301 mol->Translate(&x);
302 break;
303 }
304};
305
306/** Submenu for aligning the atoms in the molecule.
307 * \param *periode periodentafel
308 * \param *mol the molecule with all the atoms
309 */
310void AlignAtoms(periodentafel *periode, molecule *mol)
311{
312 atom *first, *second, *third;
313 vector x,n;
314 char choice; // menu choice char
315
316 cout << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
317 cout << Verbose(0) << " a - state three atoms defining align plane" << endl;
318 cout << Verbose(0) << " b - state alignment vector" << endl;
319 cout << Verbose(0) << " c - state two atoms in alignment direction" << endl;
320 cout << Verbose(0) << " d - align automatically by least square fit" << endl;
321 cout << Verbose(0) << "all else - go back" << endl;
322 cout << Verbose(0) << "===============================================" << endl;
323 cout << Verbose(0) << "INPUT: ";
324 cin >> choice;
325
326 switch (choice) {
327 default:
328 case 'a': // three atoms defining mirror plane
329 first = mol->AskAtom("Enter first atom: ");
330 second = mol->AskAtom("Enter second atom: ");
331 third = mol->AskAtom("Enter third atom: ");
332
333 n.MakeNormalVector((const vector *)&first->x,(const vector *)&second->x,(const vector *)&third->x);
334 break;
335 case 'b': // normal vector of mirror plane
336 cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
337 n.AskPosition(mol->cell_size,0);
338 n.Normalize();
339 break;
340 case 'c': // three atoms defining mirror plane
341 first = mol->AskAtom("Enter first atom: ");
342 second = mol->AskAtom("Enter second atom: ");
343
344 n.CopyVector((const vector *)&first->x);
345 n.SubtractVector((const vector *)&second->x);
346 n.Normalize();
347 break;
348 case 'd':
349 char shorthand[4];
350 vector a;
351 struct lsq_params param;
352 do {
353 fprintf(stdout, "Enter the element of atoms to be chosen: ");
354 fscanf(stdin, "%3s", shorthand);
355 } while ((param.type = periode->FindElement(shorthand)) == NULL);
356 cout << Verbose(0) << "Element is " << param.type->name << endl;
357 mol->GetAlignVector(&param);
358 for (int i=0;i<3;i++) {
359 x.x[i] = gsl_vector_get(param.x,i);
360 n.x[i] = gsl_vector_get(param.x,i+3);
361 }
362 gsl_vector_free(param.x);
363 cout << Verbose(0) << "Offset vector: ";
364 x.Output((ofstream *)&cout);
365 cout << Verbose(0) << endl;
366 n.Normalize();
367 break;
368 };
369 cout << Verbose(0) << "Alignment vector: ";
370 n.Output((ofstream *)&cout);
371 cout << Verbose(0) << endl;
372 mol->Align(&n);
373};
374
375/** Submenu for mirroring the atoms in the molecule.
376 * \param *mol the molecule with all the atoms
377 */
378void MirrorAtoms(molecule *mol)
379{
380 atom *first, *second, *third;
381 vector n;
382 char choice; // menu choice char
383
384 cout << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;
385 cout << Verbose(0) << " a - state three atoms defining mirror plane" << endl;
386 cout << Verbose(0) << " b - state normal vector of mirror plane" << endl;
387 cout << Verbose(0) << " c - state two atoms in normal direction" << endl;
388 cout << Verbose(0) << "all else - go back" << endl;
389 cout << Verbose(0) << "===============================================" << endl;
390 cout << Verbose(0) << "INPUT: ";
391 cin >> choice;
392
393 switch (choice) {
394 default:
395 case 'a': // three atoms defining mirror plane
396 first = mol->AskAtom("Enter first atom: ");
397 second = mol->AskAtom("Enter second atom: ");
398 third = mol->AskAtom("Enter third atom: ");
399
400 n.MakeNormalVector((const vector *)&first->x,(const vector *)&second->x,(const vector *)&third->x);
401 break;
402 case 'b': // normal vector of mirror plane
403 cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
404 n.AskPosition(mol->cell_size,0);
405 n.Normalize();
406 break;
407 case 'c': // three atoms defining mirror plane
408 first = mol->AskAtom("Enter first atom: ");
409 second = mol->AskAtom("Enter second atom: ");
410
411 n.CopyVector((const vector *)&first->x);
412 n.SubtractVector((const vector *)&second->x);
413 n.Normalize();
414 break;
415 };
416 cout << Verbose(0) << "Normal vector: ";
417 n.Output((ofstream *)&cout);
418 cout << Verbose(0) << endl;
419 mol->Mirror((const vector *)&n);
420};
421
422/** Submenu for removing the atoms from the molecule.
423 * \param *mol the molecule with all the atoms
424 */
425void RemoveAtoms(molecule *mol)
426{
427 atom *first, *second;
428 int axis;
429 double tmp1, tmp2;
430 char choice; // menu choice char
431
432 cout << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;
433 cout << Verbose(0) << " a - state atom for removal by number" << endl;
434 cout << Verbose(0) << " b - keep only in radius around atom" << endl;
435 cout << Verbose(0) << " c - remove this with one axis greater value" << endl;
436 cout << Verbose(0) << "all else - go back" << endl;
437 cout << Verbose(0) << "===============================================" << endl;
438 cout << Verbose(0) << "INPUT: ";
439 cin >> choice;
440
441 switch (choice) {
442 default:
443 case 'a':
444 if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: ")))
445 cout << Verbose(1) << "Atom removed." << endl;
446 else
447 cout << Verbose(1) << "Atom not found." << endl;
448 break;
449 case 'b':
450 second = mol->AskAtom("Enter number of atom as reference point: ");
451 cout << Verbose(0) << "Enter radius: ";
452 cin >> tmp1;
453 first = mol->start;
454 while(first->next != mol->end) {
455 first = first->next;
456 if (first->x.Distance((const vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
457 mol->RemoveAtom(first);
458 }
459 break;
460 case 'c':
461 cout << Verbose(0) << "Which axis is it: ";
462 cin >> axis;
463 cout << Verbose(0) << "Left inward boundary: ";
464 cin >> tmp1;
465 cout << Verbose(0) << "Right inward boundary: ";
466 cin >> tmp2;
467 first = mol->start;
468 while(first->next != mol->end) {
469 first = first->next;
470 if ((first->x.x[axis] > tmp2) || (first->x.x[axis] < tmp1)) // out of boundary ...
471 mol->RemoveAtom(first);
472 }
473 break;
474 };
475 //mol->Output((ofstream *)&cout);
476 choice = 'r';
477};
478
479/** Submenu for measuring out the atoms in the molecule.
480 * \param *periode periodentafel
481 * \param *mol the molecule with all the atoms
482 */
483void MeasureAtoms(periodentafel *periode, molecule *mol)
484{
485 atom *first, *second, *third;
486 vector x,y;
487 double min[256], tmp1, tmp2, tmp3;
488 int Z;
489 char choice; // menu choice char
490
491 cout << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;
492 cout << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;
493 cout << Verbose(0) << " b - calculate bond length between two atoms" << endl;
494 cout << Verbose(0) << " c - calculate bond angle" << endl;
495 cout << Verbose(0) << "all else - go back" << endl;
496 cout << Verbose(0) << "===============================================" << endl;
497 cout << Verbose(0) << "INPUT: ";
498 cin >> choice;
499
500 switch(choice) {
501 default:
502 cout << Verbose(1) << "Not a valid choice." << endl;
503 break;
504 case 'a':
505 first = mol->AskAtom("Enter first atom: ");
506 for (int i=0;i<256;i++)
507 min[i] = 0.;
508
509 second = mol->start;
510 while ((second->next != mol->end)) {
511 second = second->next; // advance
512 Z = second->type->Z;
513 tmp1 = 0.;
514 if (first != second) {
515 x.CopyVector((const vector *)&first->x);
516 x.SubtractVector((const vector *)&second->x);
517 tmp1 = x.Norm();
518 }
519 if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
520 //cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
521 }
522 for (int i=0;i<256;i++)
523 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;
524 break;
525
526 case 'b':
527 first = mol->AskAtom("Enter first atom: ");
528 second = mol->AskAtom("Enter second atom: ");
529 for (int i=0;i<NDIM;i++)
530 min[i] = 0.;
531 x.CopyVector((const vector *)&first->x);
532 x.SubtractVector((const vector *)&second->x);
533 tmp1 = x.Norm();
534 cout << Verbose(1) << "Distance vector is ";
535 x.Output((ofstream *)&cout);
536 cout << "." << endl << "Norm of distance is " << tmp1 << "." << endl;
537 break;
538
539 case 'c':
540 cout << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl;
541 first = mol->AskAtom("Enter first atom: ");
542 second = mol->AskAtom("Enter central atom: ");
543 third = mol->AskAtom("Enter last atom: ");
544 tmp1 = tmp2 = tmp3 = 0.;
545 x.CopyVector((const vector *)&first->x);
546 x.SubtractVector((const vector *)&second->x);
547 y.CopyVector((const vector *)&third->x);
548 y.SubtractVector((const vector *)&second->x);
549 cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
550 cout << Verbose(0) << (acos(x.ScalarProduct((const vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
551 break;
552 }
553};
554
555/** Submenu for measuring out the atoms in the molecule.
556 * \param *mol the molecule with all the atoms
557 * \param *configuration configuration structure for the to be written config files of all fragments
558 */
559void FragmentAtoms(molecule *mol, config *configuration)
560{
561 enum BondOrderScheme Scheme = NoScheme;
562 enum CutCyclicBond CutCyclic;
563 char schema;
564 int Order1, Order2;
565 clock_t start, end;
566
567 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
568 cout << Verbose(0) << "What's the desired bond order: ";
569 cin >> Order1;
570 cout << Verbose(0) << "What's the desired bond scheme [(B)ottomUp/(T)opDown/(A)NOVA/(C)ombined]: ";
571 cin >> schema;
572 CutCyclic = SaturateBond;
573 switch (schema) {
574 default:
575 Scheme = NoScheme;
576 break;
577 case 'B':
578 CutCyclic = KeepBond;
579 case 'b':
580 Scheme = BottomUp;
581 Order2 = 0;
582 break;
583 case 'T':
584 CutCyclic = KeepBond;
585 case 't':
586 Scheme = TopDown;
587 Order2 = Order1;
588 Order1 = 0;
589 break;
590 case 'A':
591 CutCyclic = KeepBond;
592 case 'a':
593 Scheme = ANOVA;
594 Order2 = 0;
595 break;
596 case 'C':
597 CutCyclic = KeepBond;
598 case 'c':
599 cout << Verbose(0) << "Combined first selects subgraphs by BottomUp of bond order " << Order1 << " then TopDown to fragment these." << endl;
600 cout << Verbose(0) << "What's the desired bond order for TopDown: ";
601 cin >> Order2;
602 Scheme = Combined;
603 break;
604 };
605 if (mol->first->next != mol->last) { // there are bonds
606 start = clock();
607 mol->FragmentMolecule((ofstream *)&cout, Order1, Order2, Scheme, configuration, CutCyclic);
608 end = clock();
609 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
610 } else
611 cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
612};
613
614/********************************************** Test routine **************************************/
615
616/** Is called always as option 'T' in the menu.
617 */
618void testroutine(molecule *mol)
619{
620 // the current test routine checks the functionality of the KeySet&Graph concept:
621 // We want to have a multiindex (the KeySet) describing a unique subgraph
622 atom *Walker = mol->start;
623 int i, comp, counter=0;
624
625 // generate some KeySets
626 cout << "Generating KeySets." << endl;
627 KeySet TestSets[mol->AtomCount+1];
628 i=1;
629 while (Walker->next != mol->end) {
630 Walker = Walker->next;
631 for (int j=0;j<i;j++) {
632 TestSets[j].insert(Walker->nr);
633 }
634 i++;
635 }
636 cout << "Testing insertion of already present item in KeySets." << endl;
637 KeySetTestPair test;
638 test = TestSets[mol->AtomCount-1].insert(Walker->nr);
639 if (test.second) {
640 cout << Verbose(1) << "Insertion worked?!" << endl;
641 } else {
642 cout << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
643 }
644 TestSets[mol->AtomCount].insert(mol->end->previous->nr);
645 TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
646
647 // constructing Graph structure
648 cout << "Generating Subgraph class." << endl;
649 Graph Subgraphs;
650
651 // insert KeySets into Subgraphs
652 cout << "Inserting KeySets into Subgraph class." << endl;
653 for (int j=0;j<mol->AtomCount;j++) {
654 Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
655 }
656 cout << "Testing insertion of already present item in Subgraph." << endl;
657 GraphTestPair test2;
658 test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.)));
659 if (test2.second) {
660 cout << Verbose(1) << "Insertion worked?!" << endl;
661 } else {
662 cout << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;
663 }
664
665 // show graphs
666 cout << "Showing Subgraph's contents, checking that it's sorted." << endl;
667 Graph::iterator A = Subgraphs.begin();
668 while (A != Subgraphs.end()) {
669 cout << (*A).second.first << ": ";
670 KeySet::iterator key = (*A).first.begin();
671 comp = -1;
672 while (key != (*A).first.end()) {
673 if ((*key) > comp)
674 cout << (*key) << " ";
675 else
676 cout << (*key) << "! ";
677 comp = (*key);
678 key++;
679 }
680 cout << endl;
681 A++;
682 }
683};
684
685/********************************************** Main routine **************************************/
686
687int main(int argc, char **argv)
688{
689 periodentafel *periode = new periodentafel; // and a period table of all elements
690 molecule *mol = new molecule(periode); // first we need an empty molecule
691 double tmp1;
692 double bond, min_bond;
693 atom *first, *second;
694 element *finder;
695 char choice; // menu choice char
696 vector x,y,z,n; // coordinates for absolute point in cell volume
697 double *factor; // unit factor if desired
698 bool valid; // flag if input was valid or not
699 config configuration;
700 ifstream test;
701 ofstream output;
702 string line;
703 char filename[255];
704 int flag = 1;
705 int Z;
706 int j, axis, count, faktor;
707 int MinimumRingSize = -1;
708 enum BondOrderScheme Scheme = NoScheme;
709 enum CutCyclicBond CutCyclic;
710 MoleculeLeafClass *Subgraphs = NULL;
711 clock_t start,end;
712 element **Elements;
713 vector **Vectors;
714
715 // Greetings
716 cout << Verbose(0) << "Molecuilder v1" << endl << endl;
717 cout << Verbose(0) << "Build your own molecule position set." << endl;
718
719 // load element list
720 periode->LoadPeriodentafel();
721 cout << Verbose(1) << "Element list loaded." << endl;
722
723 // parse config file if given
724 if (argc >= 2) { // config file specified as option
725 cout << Verbose(0) << "Config file given." << endl;
726 test.open(argv[1], ios::in);
727 if (test == NULL) {
728 cout << Verbose(1) << "Specified config file not found." << endl;
729 //return (1);
730 } else {
731 cout << Verbose(1) << "Specified config file found, parsing ...";
732 switch (configuration.TestSyntax(&test, periode, mol)) {
733 case 1:
734 cout << "new syntax." << endl;
735 configuration.Load(&test, periode, mol);
736 break;
737 case 0:
738 cout << "old syntax." << endl;
739 configuration.LoadOld(&test, periode, mol);
740 break;
741 default:
742 cout << "unknown syntax or empty configuration file." << endl;
743 //return(1);
744 }
745 test.close();
746 }
747 }
748 // command line stuff
749 if (argc > 2) {
750 cout << Verbose(0) << "Recognized command line argument: " << argv[2][1] << ".\n";
751 switch(argv[2][1]) {
752 case 't':
753 cout << Verbose(1) << "Translating all ions to new origin." << endl;
754 for (int i=0;i<3;i++)
755 x.x[i] = atof(argv[3+i]);
756 mol->Translate((const vector *)&x);
757 break;
758 case 'a':
759 cout << Verbose(1) << "Adding new atom." << endl;
760 first = new atom;
761 for (int i=0;i<3;i++)
762 first->x.x[i] = atof(argv[4+i]);
763 finder = periode->start;
764 while (finder != periode->end) {
765 finder = finder->next;
766 if (strncmp(finder->symbol,argv[3],3) == 0) {
767 first->type = finder;
768 break;
769 }
770 }
771 mol->AddAtom(first); // add to molecule
772 break;
773 case 's':
774 j = -1;
775 cout << Verbose(1) << "Scaling all ion positions by factor." << endl;
776 factor = (double *) Malloc(sizeof(double)*NDIM, "main: *factor");
777 factor[0] = atof(argv[3]);
778 factor[1] = (argc >=4) ? atof(argv[4]) : atof(argv[3]);
779 factor[2] = (argc >=5) ? atof(argv[5]) : atof(argv[3]);
780 mol->Scale(&factor);
781 for (int i=0;i<3;i++) {
782 j += i+1;
783 x.x[i] = atof(argv[3+i]);
784 mol->cell_size[j]*=factor[i];
785 }
786 Free((void **)&factor, "main: *factor");
787 break;
788 case 'c':
789 j = -1;
790 cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
791 // make every coordinate positive
792 mol->CenterEdge((ofstream *)&cout, &x);
793 // update Box of atoms by boundary
794 mol->SetBoxDimension(&x);
795 // translate each coordinate by boundary
796 j=-1;
797 for (int i=0;i<3;i++) {
798 j += i+1;
799 x.x[i] = atof(argv[3+i]);
800 mol->cell_size[j] += x.x[i]*2.;
801 }
802 mol->Translate((const vector *)&x);
803 break;
804 case 'r':
805 cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
806 break;
807 case 'p':
808 cout << Verbose(1) << "Parsing xyz file for new atoms." << endl;
809 if (!mol->AddXYZFile(argv[3]))
810 cout << Verbose(2) << "File not found." << endl;
811 else
812 cout << Verbose(2) << "File found and parsed." << endl;
813 break;
814 case 'f':
815 int i,j;
816 flag = 0;
817 if (argc > 4) {
818 cout << Verbose(0) << "Creating connection matrix..." << endl;
819 start = clock();
820 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[3]));
821 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
822 cout << Verbose(0) << "What's the desired bond scheme [(B)ottomUp/(T)opDown/(A)NOVA/(C)ombined]: ";
823 CutCyclic = SaturateBond;
824 switch (argv[5][0]) {
825 default:
826 Scheme = NoScheme;
827 break;
828 case 'B':
829 CutCyclic = KeepBond;
830 case 'b':
831 Scheme = BottomUp;
832 i = atoi(argv[4]);
833 j = 0;
834 break;
835 case 'T':
836 CutCyclic = KeepBond;
837 case 't':
838 Scheme = TopDown;
839 i = 0;
840 j = atoi(argv[4]);
841 break;
842 case 'A':
843 CutCyclic = KeepBond;
844 case 'a':
845 Scheme = ANOVA;
846 i = atoi(argv[4]);
847 j = 0;
848 break;
849 case 'C':
850 CutCyclic = KeepBond;
851 case 'c':
852 if (argc > 5) {
853 i = atoi(argv[4]);
854 j = atoi(argv[6]);
855 Scheme = Combined;
856 }
857 else
858 cerr << "Missing second bond order for TopDown fragmentation in combined approach." << endl;
859 break;
860 };
861 if (mol->first->next != mol->last) {
862 mol->FragmentMolecule((ofstream *)&cout, i, j, Scheme, &configuration, CutCyclic);
863 }
864 end = clock();
865 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
866 }
867 break;
868 default:
869 break;
870 }
871 if (flag) {
872 mol->CountElements(); // this is a bugfix, atoms should should actually be added correctly to this fragment
873 mol->CalculateOrbitals(configuration);
874 configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.MaxPsiDouble;
875 if (argc >= 1)
876 output.open(argv[1], ios::trunc);
877 if (!output)
878 output.open("main_pcp_linux", ios::trunc);
879 if (configuration.Save(&output, periode, mol))
880 cout << Verbose(1) << "Saving of config file successful." << endl;
881 else
882 cout << Verbose(1) << "Saving of config file failed." << endl;
883 output.close();
884 output.clear();
885 // and save to xyz file
886 if (argc >= 1) {
887 strcpy(filename, argv[1]);
888 strcat(filename, ".xyz");
889 output.open(filename, ios::trunc);
890 }
891 if (output == NULL) {
892 strcpy(filename,"main_pcp_linux");
893 strcat(filename, ".xyz");
894 output.open(filename, ios::trunc);
895 }
896 if (mol->OutputXYZ(&output))
897 cout << Verbose(0) << "Saving of XYZ file successful." << endl;
898 else
899 cout << Verbose(0) << "Saving of XYZ file failed." << endl;
900 output.close();
901 output.clear();
902 }
903 delete(mol);
904 delete(periode);
905 return (0);
906 }
907
908
909 // General stuff
910 if (mol->cell_size[0] == 0.) {
911 cout << Verbose(0) << "enter lower triadiagonal form of basis matrix" << endl << endl;
912 for (int i=0;i<6;i++) {
913 cout << Verbose(1) << "Cell size" << i << ": ";
914 cin >> mol->cell_size[i];
915 }
916 }
917
918 // now the main construction loop
919 cout << Verbose(0) << endl << "Now comes the real construction..." << endl;
920 do {
921 cout << Verbose(0) << endl << endl;
922 cout << Verbose(0) << "============Element list=======================" << endl;
923 mol->Checkout((ofstream *)&cout);
924 cout << Verbose(0) << "============Atom list==========================" << endl;
925 mol->Output((ofstream *)&cout);
926 cout << Verbose(0) << "============Menu===============================" << endl;
927 cout << Verbose(0) << "a - add an atom" << endl;
928 cout << Verbose(0) << "r - remove an atom" << endl;
929 cout << Verbose(0) << "b - scale a bond between atoms" << endl;
930 cout << Verbose(0) << "u - change an atoms element" << endl;
931 cout << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
932 cout << Verbose(0) << "-----------------------------------------------" << endl;
933 cout << Verbose(0) << "p - Parse xyz file" << endl;
934 cout << Verbose(0) << "e - edit the current configuration" << endl;
935 cout << Verbose(0) << "o - create connection matrix" << endl;
936 cout << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;
937 cout << Verbose(0) << "-----------------------------------------------" << endl;
938 cout << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
939 cout << Verbose(0) << "i - realign molecule" << endl;
940 cout << Verbose(0) << "m - mirror all molecules" << endl;
941 cout << Verbose(0) << "t - translate molecule by vector" << endl;
942 cout << Verbose(0) << "c - scale by unit transformation" << endl;
943 cout << Verbose(0) << "g - center atoms in box" << endl;
944 cout << Verbose(0) << "-----------------------------------------------" << endl;
945 cout << Verbose(0) << "s - save current setup to config file" << endl;
946 cout << Verbose(0) << "T - call the current test routine" << endl;
947 cout << Verbose(0) << "q - quit" << endl;
948 cout << Verbose(0) << "===============================================" << endl;
949 cout << Verbose(0) << "Input: ";
950 cin >> choice;
951
952 switch (choice) {
953 default:
954 case 'q': // quit
955 break;
956
957 case 'a': // add atom
958 AddAtoms(periode, mol);
959 choice = 'a';
960 break;
961
962 case 'd': // duplicate the periodic cell along a given axis, given times
963 cout << Verbose(0) << "State the axis [(+-)123]: ";
964 cin >> axis;
965 cout << Verbose(0) << "State the factor: ";
966 cin >> faktor;
967
968 mol->CountAtoms((ofstream *)&cout); // recount atoms
969 if (mol->AtomCount != 0) { // if there is more than none
970 count = mol->AtomCount; // is changed becausing of adding, thus has to be stored away beforehand
971 Elements = (element **) Malloc(sizeof(element *)*count, "main: duplicateCell - **Elements");
972 Vectors = (vector **) Malloc(sizeof(vector *)*count, "main: duplicateCell - **Vectors");
973 j = 0;
974 first = mol->start;
975 while (first->next != mol->end) { // make a list of all atoms with coordinates and element
976 first = first->next;
977 Elements[j] = first->type;
978 Vectors[j] = &first->x;
979 j++;
980 }
981 if (count != j)
982 cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
983 x.Zero();
984 y.Zero();
985 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
986 for (int i=1;i<faktor;i++) { // then add this list with respective translation factor times
987 x.AddVector(&y); // per factor one cell width further
988 for (int k=0;k<count;k++) { // go through every atom of the original cell
989 first = new atom(); // create a new body
990 first->x.CopyVector(Vectors[k]); // use coordinate of original atom
991 first->x.AddVector(&x); // translate the coordinates
992 first->type = Elements[k]; // insert original element
993 mol->AddAtom(first); // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
994 }
995 }
996 if (mol->first->next != mol->last) // if connect matrix is present already, redo it
997 mol->CreateAdjacencyList((ofstream *)&cout, mol->BondDistance);
998 // free memory
999 Free((void **)&Elements, "main: duplicateCell - **Elements");
1000 Free((void **)&Vectors, "main: duplicateCell - **Vectors");
1001 // correct cell size
1002 if (axis < 0) { // if sign was negative, we have to translate everything
1003 x.Zero();
1004 x.AddVector(&y);
1005 x.Scale(-(faktor-1));
1006 mol->Translate(&x);
1007 }
1008 mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
1009 }
1010 break;
1011
1012 case 'g': // center the atoms
1013 CenterAtoms(mol);
1014 break;
1015
1016 case 'b': // scale a bond
1017 cout << Verbose(0) << "Scaling bond length between two atoms." << endl;
1018 first = mol->AskAtom("Enter first (fixed) atom: ");
1019 second = mol->AskAtom("Enter second (shifting) atom: ");
1020 min_bond = 0.;
1021 for (int i=0;i<3;i++)
1022 min_bond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
1023 min_bond = sqrt(min_bond);
1024 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;
1025 cout << Verbose(0) << "Enter new bond length [a.u.]: ";
1026 cin >> bond;
1027 for (int i=0;i<3;i++) {
1028 second->x.x[i] -= (second->x.x[i]-first->x.x[i])/min_bond*(min_bond-bond);
1029 }
1030 //cout << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
1031 //second->Output(second->type->No, 1, (ofstream *)&cout);
1032 break;
1033
1034 case 'i': // align all atoms
1035 AlignAtoms(periode, mol);
1036 break;
1037
1038 case 'm': // mirror atoms along a given axis
1039 MirrorAtoms(mol);
1040 break;
1041
1042 case 't': // translate all atoms
1043 cout << Verbose(0) << "Enter translation vector." << endl;
1044 x.AskPosition(mol->cell_size,0);
1045 mol->Translate((const vector *)&x);
1046 break;
1047
1048 case 'e': // edit each field of the configuration
1049 configuration.Edit(mol);
1050 break;
1051
1052 case 'c': // unit scaling of the metric
1053 cout << Verbose(0) << "Enter three factors: ";
1054 factor = (double *) Malloc(sizeof(double)*NDIM, "main: *factor");
1055 cin >> factor[0];
1056 cin >> factor[1];
1057 cin >> factor[2];
1058 valid = true;
1059 mol->Scale(&factor);
1060 Free((void **)&factor, "main: *factor");
1061 break;
1062
1063 case 'r': // remove atom
1064 RemoveAtoms(mol);
1065 break;
1066
1067 case 'l': // measure distances or angles
1068 MeasureAtoms(periode, mol);
1069 break;
1070
1071 case 'p': // parse and XYZ file
1072 cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
1073 do {
1074 cout << Verbose(0) << "Enter file name: ";
1075 cin >> filename;
1076 } while (!mol->AddXYZFile(filename));
1077 break;
1078
1079 case 'o': // create the connection matrix
1080 cout << Verbose(0) << "What's the maximum bond distance: ";
1081 cin >> tmp1;
1082 start = clock();
1083 mol->CreateAdjacencyList((ofstream *)&cout, tmp1);
1084 //mol->CreateListOfBondsPerAtom((ofstream *)&cout);
1085 Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, false, MinimumRingSize);
1086 while (Subgraphs->next != NULL) {
1087 Subgraphs = Subgraphs->next;
1088 delete(Subgraphs->previous);
1089 }
1090 delete(Subgraphs); // we don't need the list here, so free everything
1091 Subgraphs = NULL;
1092 end = clock();
1093 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
1094 break;
1095
1096 case 'f':
1097 FragmentAtoms(mol, &configuration);
1098 break;
1099
1100 case 'u': // change an atom's element
1101 first = NULL;
1102 do {
1103 cout << Verbose(0) << "Change the element of which atom: ";
1104 cin >> Z;
1105 } while ((first = mol->FindAtom(Z)) == NULL);
1106 cout << Verbose(0) << "New element by atomic number Z: ";
1107 cin >> Z;
1108 first->type = periode->FindElement(Z);
1109 cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
1110 break;
1111
1112 case 'T':
1113 testroutine(mol);
1114 break;
1115
1116 case 's': // save to config file
1117 // get correct valence orbitals
1118 mol->CalculateOrbitals(configuration);
1119 configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.MaxPsiDouble;
1120 if (argc >= 1)
1121 output.open(argv[1], ios::trunc);
1122 if (output == NULL)
1123 output.open("main_pcp_linux", ios::trunc);
1124 if (configuration.Save(&output, periode, mol))
1125 cout << Verbose(0) << "Saving of config file successful." << endl;
1126 else
1127 cout << Verbose(0) << "Saving of config file failed." << endl;
1128 output.close();
1129 output.clear();
1130 // and save to xyz file
1131 if (argc >= 1) {
1132 strcpy(filename, argv[1]);
1133 strcat(filename, ".xyz");
1134 output.open(filename, ios::trunc);
1135 }
1136 if (output == NULL) {
1137 strcpy(filename,"main_pcp_linux");
1138 strcat(filename, ".xyz");
1139 output.open(filename, ios::trunc);
1140 }
1141 if (mol->OutputXYZ(&output))
1142 cout << Verbose(0) << "Saving of XYZ file successful." << endl;
1143 else
1144 cout << Verbose(0) << "Saving of XYZ file failed." << endl;
1145 output.close();
1146 output.clear();
1147 break;
1148 };
1149 } while (choice != 'q');
1150
1151 // save element data base
1152 if (periode->StorePeriodentafel())
1153 cout << Verbose(0) << "Saving of elements.db successful." << endl;
1154 else
1155 cout << Verbose(0) << "Saving of elements.db failed." << endl;
1156
1157 // Free all
1158 if (Subgraphs != NULL) { // free disconnected subgraph list of DFS analysis was performed
1159 while (Subgraphs->next != NULL) {
1160 Subgraphs = Subgraphs->next;
1161 delete(Subgraphs->previous);
1162 }
1163 delete(Subgraphs);
1164 }
1165 delete(mol);
1166 delete(periode);
1167 return (0);
1168}
1169
1170/********************************************** E N D **************************************************/
Note: See TracBrowser for help on using the repository browser.