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