source: src/builder.cpp@ a1fe77

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

Basic implementation of Constrained MD is done, missing testing.

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