source: src/builder.cpp@ d8b94a

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 Candidate_v1.7.0 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 d8b94a was d8b94a, checked in by Frederik Heber <heber@…>, 17 years ago

Gaussian basis for MPQC input files can now be specified with -B switch

builder.cpp: -B in ParseCommandLineOptions
molecules.hpp: configuration::basis added
config.cpp: SaveMPQC() takes configuration::basis as name
  • Property mode set to 100755
File size: 61.3 KB
RevLine 
[14de469]1/** \file builder.cpp
[a8bcea6]2 *
[14de469]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
[a8bcea6]8 *
[14de469]9 */
10
11/*! \mainpage Molecuilder - a molecular set builder
[a8bcea6]12 *
[14de469]13 * This introductory shall briefly make aquainted with the program, helping in installing and a first run.
[a8bcea6]14 *
[14de469]15 * \section about About the Program
[a8bcea6]16 *
[6ac7ee]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.
[a8bcea6]20 *
[6ac7ee]21 * A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello
22 * molecular dynamics implementation.
[a8bcea6]23 *
[14de469]24 * \section install Installation
[a8bcea6]25 *
[6ac7ee]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
[a8bcea6]30 *
[6ac7ee]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
[a8bcea6]34 *
[14de469]35 * \section run Running
[a8bcea6]36 *
[6ac7ee]37 * The program can be executed by running: ./molecuilder
[a8bcea6]38 *
[6ac7ee]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.
[a8bcea6]42 *
[14de469]43 * \section ref References
[a8bcea6]44 *
[6ac7ee]45 * For the special configuration file format, see the documentation of pcp.
[a8bcea6]46 *
[14de469]47 */
48
49
50using namespace std;
51
[6ac7ee]52#include "boundary.hpp"
53#include "ellipsoid.hpp"
[14de469]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 */
[7f3b9d]63static void AddAtoms(periodentafel *periode, molecule *mol)
[14de469]64{
[6ac7ee]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;
[14de469]146/*
[6ac7ee]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 };
[14de469]255};
256
257/** Submenu for centering the atoms in the molecule.
258 * \param *mol the molecule with all the atoms
259 */
[7f3b9d]260static void CenterAtoms(molecule *mol)
[14de469]261{
[6ac7ee]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 }
[14de469]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 */
[7f3b9d]318static void AlignAtoms(periodentafel *periode, molecule *mol)
[14de469]319{
[6ac7ee]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);
[14de469]381};
382
383/** Submenu for mirroring the atoms in the molecule.
384 * \param *mol the molecule with all the atoms
385 */
[7f3b9d]386static void MirrorAtoms(molecule *mol)
[14de469]387{
[6ac7ee]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);
[14de469]428};
429
430/** Submenu for removing the atoms from the molecule.
431 * \param *mol the molecule with all the atoms
432 */
[7f3b9d]433static void RemoveAtoms(molecule *mol)
[14de469]434{
[6ac7ee]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';
[14de469]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 */
[d52ea1b]491static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
[14de469]492{
[6ac7ee]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 }
[14de469]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 */
[7f3b9d]602static void FragmentAtoms(molecule *mol, config *configuration)
[14de469]603{
[6ac7ee]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;
[14de469]617};
618
619/********************************************** Test routine **************************************/
620
621/** Is called always as option 'T' in the menu.
622 */
[7f3b9d]623static void testroutine(molecule *mol)
[14de469]624{
[6ac7ee]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 }
[14de469]688};
689
[dbe929]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 */
[7f3b9d]696static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, molecule *mol)
[dbe929]697{
[6ac7ee]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 }
[dbe929]763};
764
[ca2b83]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
[d7d29c]772 * \param *PathToDatabases pointer to db's path in **argv
[ca2b83]773 * \return exit code (0 - successful, all else - something's wrong)
774 */
[d7d29c]775static int ParseCommandLineOptions(int argc, char **argv, molecule *&mol, periodentafel *&periode, config& configuration, char *&ConfigFileName, char *&PathToDatabases)
[14de469]776{
[6ac7ee]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;
[d8b94a]809 cout << "\t-B <basis>\tSetting basis to store to MPQC config files." << endl;
[6ac7ee]810 cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
811 cout << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
812 cout << "\t-O\tCenter atoms in origin." << endl;
813 cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
814 cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
815 cout << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
816 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;
817 cout << "\t-h/-H/-?\tGive this help screen." << endl;
818 cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
819 cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
820 cout << "\t-N\tGet non-convex-envelope." << endl;
821 cout << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;
822 cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
823 cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
824 cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
825 cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
826 cout << "\t-T <file> Store temperatures from the config file in <file>." << endl;
827 cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
828 cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl;
829 cout << "\t-v/-V\t\tGives version information." << endl;
830 cout << "Note: config files must not begin with '-' !" << endl;
831 delete(mol);
832 delete(periode);
833 return (1);
834 break;
835 case 'v':
836 case 'V':
837 cout << argv[0] << " " << VERSIONSTRING << endl;
838 cout << "Build your own molecule position set." << endl;
839 delete(mol);
840 delete(periode);
841 return (1);
842 break;
843 case 'e':
844 if ((argptr >= argc) || (argv[argptr][0] == '-')) {
845 cerr << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl;
846 } else {
847 cout << "Using " << argv[argptr] << " as elements database." << endl;
848 PathToDatabases = argv[argptr];
849 argptr+=1;
850 }
851 break;
852 case 'n':
853 cout << "I won't parse trajectories." << endl;
854 configuration.FastParsing = true;
855 break;
856 default: // no match? Step on
857 argptr++;
858 break;
859 }
860 } else
861 argptr++;
862 } while (argptr < argc);
863
864 // 2. Parse the element database
865 if (periode->LoadPeriodentafel(PathToDatabases)) {
866 cout << Verbose(0) << "Element list loaded successfully." << endl;
867 //periode->Output((ofstream *)&cout);
868 } else {
869 cout << Verbose(0) << "Element list loading failed." << endl;
870 return 1;
871 }
872 // 3. Find config file name and parse if possible
873 if (argv[1][0] != '-') {
874 cout << Verbose(0) << "Config file given." << endl;
875 test.open(argv[1], ios::in);
876 if (test == NULL) {
877 //return (1);
878 output.open(argv[1], ios::out);
879 if (output == NULL) {
880 cout << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl;
881 config_present = absent;
882 } else {
883 cout << "Empty configuration file." << endl;
884 ConfigFileName = argv[1];
885 config_present = empty;
886 output.close();
887 }
888 } else {
889 test.close();
890 ConfigFileName = argv[1];
891 cout << Verbose(1) << "Specified config file found, parsing ... ";
892 switch (configuration.TestSyntax(ConfigFileName, periode, mol)) {
893 case 1:
894 cout << "new syntax." << endl;
895 configuration.Load(ConfigFileName, periode, mol);
896 config_present = present;
897 break;
898 case 0:
899 cout << "old syntax." << endl;
900 configuration.LoadOld(ConfigFileName, periode, mol);
901 config_present = present;
902 break;
903 default:
904 cout << "Unknown syntax or empty, yet present file." << endl;
905 config_present = empty;
906 }
907 }
908 } else
909 config_present = absent;
910 // 4. parse again through options, now for those depending on elements db and config presence
911 argptr = 1;
912 do {
913 cout << "Current Command line argument: " << argv[argptr] << "." << endl;
914 if (argv[argptr][0] == '-') {
915 argptr++;
916 if ((config_present == present) || (config_present == empty)) {
917 switch(argv[argptr-1][1]) {
918 case 'p':
919 ExitFlag = 1;
920 if ((argptr >= argc) || (argv[argptr][0] == '-')) {
921 ExitFlag = 255;
922 cerr << "Not enough arguments for parsing: -p <xyz file>" << endl;
923 } else {
924 SaveFlag = true;
925 cout << Verbose(1) << "Parsing xyz file for new atoms." << endl;
926 if (!mol->AddXYZFile(argv[argptr]))
927 cout << Verbose(2) << "File not found." << endl;
928 else {
929 cout << Verbose(2) << "File found and parsed." << endl;
930 config_present = present;
931 }
932 }
933 break;
934 case 'a':
935 ExitFlag = 1;
936 if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1]))) {
937 ExitFlag = 255;
938 cerr << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl;
939 } else {
940 SaveFlag = true;
941 cout << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";
942 first = new atom;
943 first->type = periode->FindElement(atoi(argv[argptr]));
944 if (first->type != NULL)
945 cout << Verbose(2) << "found element " << first->type->name << endl;
946 for (int i=NDIM;i--;)
947 first->x.x[i] = atof(argv[argptr+1+i]);
948 if (first->type != NULL) {
949 mol->AddAtom(first); // add to molecule
950 if ((config_present == empty) && (mol->AtomCount != 0))
951 config_present = present;
952 } else
953 cerr << Verbose(1) << "Could not find the specified element." << endl;
954 argptr+=4;
955 }
956 break;
957 default: // no match? Don't step on (this is done in next switch's default)
958 break;
959 }
960 }
961 if (config_present == present) {
962 switch(argv[argptr-1][1]) {
[d8b94a]963 case 'B':
964 if ((argptr >= argc) || (argv[argptr][0] == '-')) {
965 ExitFlag = 255;
966 cerr << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl;
967 } else {
968 configuration.basis = argv[argptr];
969 cout << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl;
970 argptr+=1;
971 }
972 break;
[6ac7ee]973 case 'D':
974 ExitFlag = 1;
975 {
976 cout << Verbose(1) << "Depth-First-Search Analysis." << endl;
977 MoleculeLeafClass *Subgraphs = NULL; // list of subgraphs from DFS analysis
978 int *MinimumRingSize = new int[mol->AtomCount];
979 atom ***ListOfLocalAtoms = NULL;
980 int FragmentCounter = 0;
981 class StackClass<bond *> *BackEdgeStack = NULL;
982 class StackClass<bond *> *LocalBackEdgeStack = NULL;
983 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr]), configuration.GetIsAngstroem());
984 mol->CreateListOfBondsPerAtom((ofstream *)&cout);
985 Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
986 if (Subgraphs != NULL) {
987 Subgraphs->next->FillBondStructureFromReference((ofstream *)&cout, mol, (FragmentCounter = 0), ListOfLocalAtoms, false); // we want to keep the created ListOfLocalAtoms
988 while (Subgraphs->next != NULL) {
989 Subgraphs = Subgraphs->next;
990 LocalBackEdgeStack = new StackClass<bond *> (Subgraphs->Leaf->BondCount);
991 Subgraphs->Leaf->PickLocalBackEdges((ofstream *)&cout, ListOfLocalAtoms[FragmentCounter++], BackEdgeStack, LocalBackEdgeStack);
992 Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize);
993 delete(LocalBackEdgeStack);
994 delete(Subgraphs->previous);
995 }
996 delete(Subgraphs);
997 for (int i=0;i<FragmentCounter;i++)
998 Free((void **)&ListOfLocalAtoms[FragmentCounter], "ParseCommandLineOptions: **ListOfLocalAtoms[]");
999 Free((void **)&ListOfLocalAtoms, "ParseCommandLineOptions: ***ListOfLocalAtoms");
1000 }
1001 delete(BackEdgeStack);
1002 delete[](MinimumRingSize);
1003 }
1004 //argptr+=1;
1005 break;
1006 case 'E':
1007 ExitFlag = 1;
1008 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) {
1009 ExitFlag = 255;
1010 cerr << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl;
1011 } else {
1012 SaveFlag = true;
1013 cout << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl;
1014 first = mol->FindAtom(atoi(argv[argptr]));
1015 first->type = periode->FindElement(atoi(argv[argptr+1]));
1016 argptr+=2;
1017 }
1018 break;
1019 case 'A':
1020 ExitFlag = 1;
1021 if ((argptr >= argc) || (argv[argptr][0] == '-')) {
1022 ExitFlag =255;
1023 cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
1024 } else {
1025 cout << "Parsing bonds from " << argv[argptr] << "." << endl;
1026 ifstream *input = new ifstream(argv[argptr]);
1027 mol->CreateAdjacencyList2((ofstream *)&cout, input);
1028 input->close();
1029 argptr+=1;
1030 }
1031 break;
1032 case 'N':
1033 ExitFlag = 1;
1034 if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
1035 ExitFlag = 255;
1036 cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
1037 } else {
1038 class Tesselation T;
1039 int N = 15;
1040 int number = 100;
1041 string filename(argv[argptr+1]);
1042 filename.append(".csv");
1043 cout << Verbose(0) << "Evaluating non-convex envelope.";
1044 cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
1045 LinkedCell LCList(mol, atof(argv[argptr])); // \NOTE not twice the radius??
1046 Find_non_convex_border((ofstream *)&cout, mol, &T, &LCList, argv[argptr+1], atof(argv[argptr]));
1047 FindDistributionOfEllipsoids((ofstream *)&cout, &T, &LCList, N, number, filename.c_str());
1048 argptr+=2;
1049 }
1050 break;
1051 case 'T':
1052 ExitFlag = 1;
1053 if ((argptr >= argc) || (argv[argptr][0] == '-')) {
1054 ExitFlag = 255;
1055 cerr << "Not enough or invalid arguments given for storing tempature: -T <temperature file>" << endl;
1056 } else {
1057 cout << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;
1058 ofstream *output = new ofstream(argv[argptr], ios::trunc);
1059 if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output))
1060 cout << Verbose(2) << "File could not be written." << endl;
1061 else
1062 cout << Verbose(2) << "File stored." << endl;
1063 output->close();
1064 delete(output);
1065 argptr+=1;
1066 }
1067 break;
1068 case 'P':
1069 ExitFlag = 1;
1070 if ((argptr >= argc) || (argv[argptr][0] == '-')) {
1071 ExitFlag = 255;
1072 cerr << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl;
1073 } else {
1074 SaveFlag = true;
1075 cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
1076 if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem()))
1077 cout << Verbose(2) << "File not found." << endl;
1078 else
1079 cout << Verbose(2) << "File found and parsed." << endl;
1080 argptr+=1;
1081 }
1082 break;
1083 case 't':
1084 ExitFlag = 1;
1085 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
1086 ExitFlag = 255;
1087 cerr << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl;
1088 } else {
1089 ExitFlag = 1;
1090 SaveFlag = true;
1091 cout << Verbose(1) << "Translating all ions to new origin." << endl;
1092 for (int i=NDIM;i--;)
1093 x.x[i] = atof(argv[argptr+i]);
1094 mol->Translate((const Vector *)&x);
1095 argptr+=3;
1096 }
1097 break;
1098 case 's':
1099 ExitFlag = 1;
1100 if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
1101 ExitFlag = 255;
1102 cerr << "Not enough or invalid arguments given for scaling: -s <factor/[factor_x]> [factor_y] [factor_z]" << endl;
1103 } else {
1104 SaveFlag = true;
1105 j = -1;
1106 cout << Verbose(1) << "Scaling all ion positions by factor." << endl;
1107 factor = new double[NDIM];
1108 factor[0] = atof(argv[argptr]);
1109 if ((argptr < argc) && (IsValidNumber(argv[argptr])))
1110 argptr++;
1111 factor[1] = atof(argv[argptr]);
1112 if ((argptr < argc) && (IsValidNumber(argv[argptr])))
1113 argptr++;
1114 factor[2] = atof(argv[argptr]);
1115 mol->Scale(&factor);
1116 for (int i=0;i<NDIM;i++) {
1117 j += i+1;
1118 x.x[i] = atof(argv[NDIM+i]);
1119 mol->cell_size[j]*=factor[i];
1120 }
1121 delete[](factor);
1122 argptr+=1;
1123 }
1124 break;
1125 case 'b':
1126 ExitFlag = 1;
1127 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
1128 ExitFlag = 255;
1129 cerr << "Not enough or invalid arguments given for centering in box: -b <length_x> <length_y> <length_z>" << endl;
1130 } else {
1131 SaveFlag = true;
1132 j = -1;
1133 cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
1134 j=-1;
1135 for (int i=0;i<NDIM;i++) {
1136 j += i+1;
1137 x.x[i] = atof(argv[argptr++]);
1138 mol->cell_size[j] += x.x[i]*2.;
1139 }
1140 // center
1141 mol->CenterInBox((ofstream *)&cout, &x);
1142 // update Box of atoms by boundary
1143 mol->SetBoxDimension(&x);
1144 }
1145 break;
1146 case 'c':
1147 ExitFlag = 1;
1148 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
1149 ExitFlag = 255;
1150 cerr << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl;
1151 } else {
1152 SaveFlag = true;
1153 j = -1;
1154 cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
1155 // make every coordinate positive
1156 mol->CenterEdge((ofstream *)&cout, &x);
1157 // update Box of atoms by boundary
1158 mol->SetBoxDimension(&x);
1159 // translate each coordinate by boundary
1160 j=-1;
1161 for (int i=0;i<NDIM;i++) {
1162 j += i+1;
1163 x.x[i] = atof(argv[argptr++]);
1164 mol->cell_size[j] += x.x[i]*2.;
1165 }
1166 mol->Translate((const Vector *)&x);
1167 }
1168 break;
1169 case 'O':
1170 ExitFlag = 1;
1171 SaveFlag = true;
1172 cout << Verbose(1) << "Centering atoms in origin." << endl;
1173 mol->CenterOrigin((ofstream *)&cout, &x);
1174 mol->SetBoxDimension(&x);
1175 break;
1176 case 'r':
1177 ExitFlag = 1;
1178 SaveFlag = true;
1179 cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
1180 break;
1181 case 'F':
1182 case 'f':
1183 ExitFlag = 1;
1184 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {
1185 ExitFlag = 255;
1186 cerr << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl;
1187 } else {
1188 cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;
1189 cout << Verbose(0) << "Creating connection matrix..." << endl;
1190 start = clock();
1191 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr++]), configuration.GetIsAngstroem());
1192 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
1193 if (mol->first->next != mol->last) {
1194 ExitFlag = mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);
1195 }
1196 end = clock();
1197 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
1198 argptr+=2;
1199 }
1200 break;
1201 case 'm':
1202 ExitFlag = 1;
1203 j = atoi(argv[argptr++]);
1204 if ((j<0) || (j>1)) {
1205 cerr << Verbose(1) << "ERROR: Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;
1206 j = 0;
1207 }
1208 if (j) {
1209 SaveFlag = true;
1210 cout << Verbose(0) << "Converting to prinicipal axis system." << endl;
1211 } else
1212 cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
1213 mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j);
1214 break;
1215 case 'o':
1216 ExitFlag = 1;
1217 if ((argptr >= argc) || (argv[argptr][0] == '-')){
1218 ExitFlag = 255;
1219 cerr << "Not enough or invalid arguments given for convex envelope: -o <tecplot output file>" << endl;
1220 } else {
1221 cout << Verbose(0) << "Evaluating volume of the convex envelope.";
1222 cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
1223 VolumeOfConvexEnvelope((ofstream *)&cout, argv[argptr], &configuration, NULL, mol);
1224 argptr+=1;
1225 }
1226 break;
1227 case 'U':
1228 ExitFlag = 1;
1229 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
1230 ExitFlag = 255;
1231 cerr << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl;
1232 volume = -1; // for case 'u': don't print error again
1233 } else {
1234 volume = atof(argv[argptr++]);
1235 cout << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;
1236 }
1237 case 'u':
1238 ExitFlag = 1;
1239 if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
1240 if (volume != -1)
1241 ExitFlag = 255;
1242 cerr << "Not enough arguments given for suspension: -u <density>" << endl;
1243 } else {
1244 double density;
1245 SaveFlag = true;
1246 cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
1247 density = atof(argv[argptr++]);
1248 if (density < 1.0) {
1249 cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
1250 density = 1.3;
1251 }
1252// for(int i=0;i<NDIM;i++) {
1253// repetition[i] = atoi(argv[argptr++]);
1254// if (repetition[i] < 1)
1255// cerr << Verbose(0) << "ERROR: repetition value must be greater 1!" << endl;
1256// repetition[i] = 1;
1257// }
1258 PrepareClustersinWater((ofstream *)&cout, &configuration, mol, volume, density); // if volume == 0, will calculate from ConvexEnvelope
1259 }
1260 break;
1261 case 'd':
1262 ExitFlag = 1;
1263 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
1264 ExitFlag = 255;
1265 cerr << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl;
1266 } else {
1267 SaveFlag = true;
1268 for (int axis = 1; axis <= NDIM; axis++) {
1269 int faktor = atoi(argv[argptr++]);
1270 int count;
1271 element ** Elements;
1272 Vector ** vectors;
1273 if (faktor < 1) {
1274 cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;
1275 faktor = 1;
1276 }
1277 mol->CountAtoms((ofstream *)&cout); // recount atoms
1278 if (mol->AtomCount != 0) { // if there is more than none
1279 count = mol->AtomCount; // is changed becausing of adding, thus has to be stored away beforehand
1280 Elements = new element *[count];
1281 vectors = new Vector *[count];
1282 j = 0;
1283 first = mol->start;
1284 while (first->next != mol->end) { // make a list of all atoms with coordinates and element
1285 first = first->next;
1286 Elements[j] = first->type;
1287 vectors[j] = &first->x;
1288 j++;
1289 }
1290 if (count != j)
1291 cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
1292 x.Zero();
1293 y.Zero();
1294 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
1295 for (int i=1;i<faktor;i++) { // then add this list with respective translation factor times
1296 x.AddVector(&y); // per factor one cell width further
1297 for (int k=count;k--;) { // go through every atom of the original cell
1298 first = new atom(); // create a new body
1299 first->x.CopyVector(vectors[k]); // use coordinate of original atom
1300 first->x.AddVector(&x); // translate the coordinates
1301 first->type = Elements[k]; // insert original element
1302 mol->AddAtom(first); // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
1303 }
1304 }
1305 // free memory
1306 delete[](Elements);
1307 delete[](vectors);
1308 // correct cell size
1309 if (axis < 0) { // if sign was negative, we have to translate everything
1310 x.Zero();
1311 x.AddVector(&y);
1312 x.Scale(-(faktor-1));
1313 mol->Translate(&x);
1314 }
1315 mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
1316 }
1317 }
1318 }
1319 break;
1320 default: // no match? Step on
1321 if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step!
1322 argptr++;
1323 break;
1324 }
1325 }
1326 } else argptr++;
1327 } while (argptr < argc);
1328 if (SaveFlag)
1329 SaveConfig(ConfigFileName, &configuration, periode, mol);
1330 if ((ExitFlag >= 1)) {
1331 delete(mol);
1332 delete(periode);
1333 return (ExitFlag);
1334 }
1335 } else { // no arguments, hence scan the elements db
1336 if (periode->LoadPeriodentafel(PathToDatabases))
1337 cout << Verbose(0) << "Element list loaded successfully." << endl;
1338 else
1339 cout << Verbose(0) << "Element list loading failed." << endl;
1340 configuration.RetrieveConfigPathAndName("main_pcp_linux");
1341 }
1342 return(0);
[ca2b83]1343};
1344
1345/********************************************** Main routine **************************************/
[14de469]1346
[ca2b83]1347int main(int argc, char **argv)
1348{
[6ac7ee]1349 periodentafel *periode = new periodentafel; // and a period table of all elements
1350 molecule *mol = new molecule(periode); // first we need an empty molecule
1351 config configuration;
1352 double tmp1;
1353 double bond, min_bond;
1354 atom *first, *second;
1355 char choice; // menu choice char
1356 Vector x,y,z,n; // coordinates for absolute point in cell volume
[ca2b83]1357 double *factor; // unit factor if desired
[6ac7ee]1358 bool valid; // flag if input was valid or not
1359 ifstream test;
1360 ofstream output;
1361 string line;
1362 char filename[MAXSTRINGSIZE];
1363 char *ConfigFileName = NULL;
1364 char *ElementsFileName = NULL;
1365 int Z;
1366 int j, axis, count, faktor;
1367 clock_t start,end;
1368// int *MinimumRingSize = NULL;
1369 MoleculeLeafClass *Subgraphs = NULL;
1370// class StackClass<bond *> *BackEdgeStack = NULL;
1371 element **Elements;
1372 Vector **vectors;
1373
1374 // =========================== PARSE COMMAND LINE OPTIONS ====================================
1375 j = ParseCommandLineOptions(argc, argv, mol, periode, configuration, ConfigFileName, ElementsFileName);
1376 if (j == 1) return 0; // just for -v and -h options
1377 if (j) return j; // something went wrong
1378
1379 // General stuff
1380 if (mol->cell_size[0] == 0.) {
1381 cout << Verbose(0) << "enter lower triadiagonal form of basis matrix" << endl << endl;
1382 for (int i=0;i<6;i++) {
1383 cout << Verbose(1) << "Cell size" << i << ": ";
1384 cin >> mol->cell_size[i];
1385 }
1386 }
1387
1388 // =========================== START INTERACTIVE SESSION ====================================
1389
1390 // now the main construction loop
1391 cout << Verbose(0) << endl << "Now comes the real construction..." << endl;
1392 do {
1393 cout << Verbose(0) << endl << endl;
1394 cout << Verbose(0) << "============Element list=======================" << endl;
1395 mol->Checkout((ofstream *)&cout);
1396 cout << Verbose(0) << "============Atom list==========================" << endl;
1397 mol->Output((ofstream *)&cout);
1398 cout << Verbose(0) << "============Menu===============================" << endl;
1399 cout << Verbose(0) << "a - add an atom" << endl;
1400 cout << Verbose(0) << "r - remove an atom" << endl;
1401 cout << Verbose(0) << "b - scale a bond between atoms" << endl;
1402 cout << Verbose(0) << "u - change an atoms element" << endl;
1403 cout << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
1404 cout << Verbose(0) << "-----------------------------------------------" << endl;
1405 cout << Verbose(0) << "p - Parse xyz file" << endl;
1406 cout << Verbose(0) << "e - edit the current configuration" << endl;
1407 cout << Verbose(0) << "o - create connection matrix" << endl;
1408 cout << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;
1409 cout << Verbose(0) << "-----------------------------------------------" << endl;
1410 cout << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
1411 cout << Verbose(0) << "i - realign molecule" << endl;
1412 cout << Verbose(0) << "m - mirror all molecules" << endl;
1413 cout << Verbose(0) << "t - translate molecule by vector" << endl;
1414 cout << Verbose(0) << "c - scale by unit transformation" << endl;
1415 cout << Verbose(0) << "g - center atoms in box" << endl;
1416 cout << Verbose(0) << "-----------------------------------------------" << endl;
1417 cout << Verbose(0) << "s - save current setup to config file" << endl;
1418 cout << Verbose(0) << "T - call the current test routine" << endl;
1419 cout << Verbose(0) << "q - quit" << endl;
1420 cout << Verbose(0) << "===============================================" << endl;
1421 cout << Verbose(0) << "Input: ";
1422 cin >> choice;
1423
1424 switch (choice) {
1425 default:
1426 case 'a': // add atom
1427 AddAtoms(periode, mol);
1428 choice = 'a';
1429 break;
1430
1431 case 'b': // scale a bond
1432 cout << Verbose(0) << "Scaling bond length between two atoms." << endl;
1433 first = mol->AskAtom("Enter first (fixed) atom: ");
1434 second = mol->AskAtom("Enter second (shifting) atom: ");
1435 min_bond = 0.;
1436 for (int i=NDIM;i--;)
1437 min_bond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
1438 min_bond = sqrt(min_bond);
1439 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;
1440 cout << Verbose(0) << "Enter new bond length [a.u.]: ";
1441 cin >> bond;
1442 for (int i=NDIM;i--;) {
1443 second->x.x[i] -= (second->x.x[i]-first->x.x[i])/min_bond*(min_bond-bond);
1444 }
1445 //cout << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
1446 //second->Output(second->type->No, 1, (ofstream *)&cout);
1447 break;
1448
1449 case 'c': // unit scaling of the metric
1450 cout << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;
1451 cout << Verbose(0) << "Enter three factors: ";
1452 factor = new double[NDIM];
1453 cin >> factor[0];
1454 cin >> factor[1];
1455 cin >> factor[2];
1456 valid = true;
1457 mol->Scale(&factor);
1458 delete[](factor);
1459 break;
1460
1461 case 'd': // duplicate the periodic cell along a given axis, given times
1462 cout << Verbose(0) << "State the axis [(+-)123]: ";
1463 cin >> axis;
1464 cout << Verbose(0) << "State the factor: ";
1465 cin >> faktor;
1466
1467 mol->CountAtoms((ofstream *)&cout); // recount atoms
1468 if (mol->AtomCount != 0) { // if there is more than none
1469 count = mol->AtomCount; // is changed becausing of adding, thus has to be stored away beforehand
1470 Elements = new element *[count];
1471 vectors = new Vector *[count];
1472 j = 0;
1473 first = mol->start;
1474 while (first->next != mol->end) { // make a list of all atoms with coordinates and element
1475 first = first->next;
1476 Elements[j] = first->type;
1477 vectors[j] = &first->x;
1478 j++;
1479 }
1480 if (count != j)
1481 cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
1482 x.Zero();
1483 y.Zero();
1484 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
1485 for (int i=1;i<faktor;i++) { // then add this list with respective translation factor times
1486 x.AddVector(&y); // per factor one cell width further
1487 for (int k=count;k--;) { // go through every atom of the original cell
1488 first = new atom(); // create a new body
1489 first->x.CopyVector(vectors[k]); // use coordinate of original atom
1490 first->x.AddVector(&x); // translate the coordinates
1491 first->type = Elements[k]; // insert original element
1492 mol->AddAtom(first); // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
1493 }
1494 }
1495 if (mol->first->next != mol->last) // if connect matrix is present already, redo it
1496 mol->CreateAdjacencyList((ofstream *)&cout, mol->BondDistance, configuration.GetIsAngstroem());
1497 // free memory
1498 delete[](Elements);
1499 delete[](vectors);
1500 // correct cell size
1501 if (axis < 0) { // if sign was negative, we have to translate everything
1502 x.Zero();
1503 x.AddVector(&y);
1504 x.Scale(-(faktor-1));
1505 mol->Translate(&x);
1506 }
1507 mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
1508 }
1509 break;
1510
1511 case 'e': // edit each field of the configuration
1512 configuration.Edit(mol);
1513 break;
1514
1515 case 'f':
1516 FragmentAtoms(mol, &configuration);
1517 break;
1518
1519 case 'g': // center the atoms
1520 CenterAtoms(mol);
1521 break;
1522
1523 case 'i': // align all atoms
1524 AlignAtoms(periode, mol);
1525 break;
1526
1527 case 'l': // measure distances or angles
1528 MeasureAtoms(periode, mol, &configuration);
1529 break;
1530
1531 case 'm': // mirror atoms along a given axis
1532 MirrorAtoms(mol);
1533 break;
1534
1535 case 'o': // create the connection matrix
1536 {
1537 cout << Verbose(0) << "What's the maximum bond distance: ";
1538 cin >> tmp1;
1539 start = clock();
1540 mol->CreateAdjacencyList((ofstream *)&cout, tmp1, configuration.GetIsAngstroem());
1541 mol->CreateListOfBondsPerAtom((ofstream *)&cout);
1542// Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
1543// while (Subgraphs->next != NULL) {
1544// Subgraphs = Subgraphs->next;
1545// Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize);
1546// delete(Subgraphs->previous);
1547// }
1548// delete(Subgraphs); // we don't need the list here, so free everything
1549// delete[](MinimumRingSize);
1550// Subgraphs = NULL;
1551 end = clock();
1552 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
1553 }
1554 break;
1555
1556 case 'p': // parse and XYZ file
1557 cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
1558 do {
1559 cout << Verbose(0) << "Enter file name: ";
1560 cin >> filename;
1561 } while (!mol->AddXYZFile(filename));
1562 break;
1563
1564 case 'q': // quit
1565 break;
1566
1567 case 'r': // remove atom
1568 RemoveAtoms(mol);
1569 break;
1570
1571 case 's': // save to config file
1572 SaveConfig(ConfigFileName, &configuration, periode, mol);
1573 break;
1574
1575 case 't': // translate all atoms
1576 cout << Verbose(0) << "Enter translation vector." << endl;
1577 x.AskPosition(mol->cell_size,0);
1578 mol->Translate((const Vector *)&x);
1579 break;
1580
1581 case 'T':
1582 testroutine(mol);
1583 break;
1584
1585 case 'u': // change an atom's element
1586 first = NULL;
1587 do {
1588 cout << Verbose(0) << "Change the element of which atom: ";
1589 cin >> Z;
1590 } while ((first = mol->FindAtom(Z)) == NULL);
1591 cout << Verbose(0) << "New element by atomic number Z: ";
1592 cin >> Z;
1593 first->type = periode->FindElement(Z);
1594 cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
1595 break;
1596 };
1597 } while (choice != 'q');
1598
1599 // save element data base
1600 if (periode->StorePeriodentafel(ElementsFileName)) //ElementsFileName
1601 cout << Verbose(0) << "Saving of elements.db successful." << endl;
1602 else
1603 cout << Verbose(0) << "Saving of elements.db failed." << endl;
1604
1605 // Free all
1606 if (Subgraphs != NULL) { // free disconnected subgraph list of DFS analysis was performed
1607 while (Subgraphs->next != NULL) {
1608 Subgraphs = Subgraphs->next;
1609 delete(Subgraphs->previous);
1610 }
1611 delete(Subgraphs);
1612 }
1613 delete(mol);
1614 delete(periode);
1615 return (0);
[14de469]1616}
1617
1618/********************************************** E N D **************************************************/
Note: See TracBrowser for help on using the repository browser.