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 | #include <boost/bind.hpp>
|
---|
51 |
|
---|
52 | using namespace std;
|
---|
53 |
|
---|
54 | #include <cstring>
|
---|
55 | #include <cstdlib>
|
---|
56 |
|
---|
57 | #include "analysis_bonds.hpp"
|
---|
58 | #include "analysis_correlation.hpp"
|
---|
59 | #include "atom.hpp"
|
---|
60 | #include "bond.hpp"
|
---|
61 | #include "bondgraph.hpp"
|
---|
62 | #include "boundary.hpp"
|
---|
63 | #include "CommandLineParser.hpp"
|
---|
64 | #include "config.hpp"
|
---|
65 | #include "element.hpp"
|
---|
66 | #include "ellipsoid.hpp"
|
---|
67 | #include "helpers.hpp"
|
---|
68 | #include "leastsquaremin.hpp"
|
---|
69 | #include "linkedcell.hpp"
|
---|
70 | #include "log.hpp"
|
---|
71 | #include "memoryusageobserver.hpp"
|
---|
72 | #include "molecule.hpp"
|
---|
73 | #include "periodentafel.hpp"
|
---|
74 | #include "UIElements/UIFactory.hpp"
|
---|
75 | #include "UIElements/MainWindow.hpp"
|
---|
76 | #include "UIElements/Dialog.hpp"
|
---|
77 | #include "Menu/ActionMenuItem.hpp"
|
---|
78 | #include "Actions/ActionRegistry.hpp"
|
---|
79 | #include "Actions/ActionHistory.hpp"
|
---|
80 | #include "Actions/MapOfActions.hpp"
|
---|
81 | #include "Actions/MethodAction.hpp"
|
---|
82 | #include "Actions/MoleculeAction/ChangeNameAction.hpp"
|
---|
83 | #include "World.hpp"
|
---|
84 | #include "version.h"
|
---|
85 | #include "World.hpp"
|
---|
86 | #include "Helpers/MemDebug.hpp"
|
---|
87 |
|
---|
88 | /********************************************* Subsubmenu routine ************************************/
|
---|
89 | #if 0
|
---|
90 | /** Submenu for adding atoms to the molecule.
|
---|
91 | * \param *periode periodentafel
|
---|
92 | * \param *molecule molecules with atoms
|
---|
93 | */
|
---|
94 | static void AddAtoms(periodentafel *periode, molecule *mol)
|
---|
95 | {
|
---|
96 | atom *first, *second, *third, *fourth;
|
---|
97 | Vector **atoms;
|
---|
98 | Vector x,y,z,n; // coordinates for absolute point in cell volume
|
---|
99 | double a,b,c;
|
---|
100 | char choice; // menu choice char
|
---|
101 | bool valid;
|
---|
102 |
|
---|
103 | cout << Verbose(0) << "===========ADD ATOM============================" << endl;
|
---|
104 | cout << Verbose(0) << " a - state absolute coordinates of atom" << endl;
|
---|
105 | cout << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
|
---|
106 | cout << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
|
---|
107 | cout << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
|
---|
108 | cout << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
|
---|
109 | cout << Verbose(0) << "all else - go back" << endl;
|
---|
110 | cout << Verbose(0) << "===============================================" << endl;
|
---|
111 | cout << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
|
---|
112 | cout << Verbose(0) << "INPUT: ";
|
---|
113 | cin >> choice;
|
---|
114 |
|
---|
115 | switch (choice) {
|
---|
116 | default:
|
---|
117 | DoeLog(2) && (eLog()<< Verbose(2) << "Not a valid choice." << endl);
|
---|
118 | break;
|
---|
119 | case 'a': // absolute coordinates of atom
|
---|
120 | cout << Verbose(0) << "Enter absolute coordinates." << endl;
|
---|
121 | first = new atom;
|
---|
122 | first->x.AskPosition(World::getInstance().getDomain(), false);
|
---|
123 | first->type = periode->AskElement(); // give type
|
---|
124 | mol->AddAtom(first); // add to molecule
|
---|
125 | break;
|
---|
126 |
|
---|
127 | case 'b': // relative coordinates of atom wrt to reference point
|
---|
128 | first = new atom;
|
---|
129 | valid = true;
|
---|
130 | do {
|
---|
131 | if (!valid) DoeLog(2) && (eLog()<< Verbose(2) << "Resulting position out of cell." << endl);
|
---|
132 | cout << Verbose(0) << "Enter reference coordinates." << endl;
|
---|
133 | x.AskPosition(World::getInstance().getDomain(), true);
|
---|
134 | cout << Verbose(0) << "Enter relative coordinates." << endl;
|
---|
135 | first->x.AskPosition(World::getInstance().getDomain(), false);
|
---|
136 | first->x.AddVector((const Vector *)&x);
|
---|
137 | cout << Verbose(0) << "\n";
|
---|
138 | } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
|
---|
139 | first->type = periode->AskElement(); // give type
|
---|
140 | mol->AddAtom(first); // add to molecule
|
---|
141 | break;
|
---|
142 |
|
---|
143 | case 'c': // relative coordinates of atom wrt to already placed atom
|
---|
144 | first = new atom;
|
---|
145 | valid = true;
|
---|
146 | do {
|
---|
147 | if (!valid) DoeLog(2) && (eLog()<< Verbose(2) << "Resulting position out of cell." << endl);
|
---|
148 | second = mol->AskAtom("Enter atom number: ");
|
---|
149 | DoLog(0) && (Log() << Verbose(0) << "Enter relative coordinates." << endl);
|
---|
150 | first->x.AskPosition(World::getInstance().getDomain(), false);
|
---|
151 | for (int i=NDIM;i--;) {
|
---|
152 | first->x.x[i] += second->x.x[i];
|
---|
153 | }
|
---|
154 | } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
|
---|
155 | first->type = periode->AskElement(); // give type
|
---|
156 | mol->AddAtom(first); // add to molecule
|
---|
157 | break;
|
---|
158 |
|
---|
159 | case 'd': // two atoms, two angles and a distance
|
---|
160 | first = new atom;
|
---|
161 | valid = true;
|
---|
162 | do {
|
---|
163 | if (!valid) {
|
---|
164 | DoeLog(2) && (eLog()<< Verbose(2) << "Resulting coordinates out of cell - " << first->x << endl);
|
---|
165 | }
|
---|
166 | cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
|
---|
167 | second = mol->AskAtom("Enter central atom: ");
|
---|
168 | third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
|
---|
169 | fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
|
---|
170 | a = ask_value("Enter distance between central (first) and new atom: ");
|
---|
171 | b = ask_value("Enter angle between new, first and second atom (degrees): ");
|
---|
172 | b *= M_PI/180.;
|
---|
173 | bound(&b, 0., 2.*M_PI);
|
---|
174 | c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
|
---|
175 | c *= M_PI/180.;
|
---|
176 | bound(&c, -M_PI, M_PI);
|
---|
177 | cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
|
---|
178 | /*
|
---|
179 | second->Output(1,1,(ofstream *)&cout);
|
---|
180 | third->Output(1,2,(ofstream *)&cout);
|
---|
181 | fourth->Output(1,3,(ofstream *)&cout);
|
---|
182 | n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
|
---|
183 | x.Copyvector(&second->x);
|
---|
184 | x.SubtractVector(&third->x);
|
---|
185 | x.Copyvector(&fourth->x);
|
---|
186 | x.SubtractVector(&third->x);
|
---|
187 |
|
---|
188 | if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
|
---|
189 | coutg() << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
|
---|
190 | continue;
|
---|
191 | }
|
---|
192 | DoLog(0) && (Log() << Verbose(0) << "resulting relative coordinates: ");
|
---|
193 | z.Output();
|
---|
194 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
195 | */
|
---|
196 | // calc axis vector
|
---|
197 | x.CopyVector(&second->x);
|
---|
198 | x.SubtractVector(&third->x);
|
---|
199 | x.Normalize();
|
---|
200 | Log() << Verbose(0) << "x: ",
|
---|
201 | x.Output();
|
---|
202 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
203 | z.MakeNormalVector(&second->x,&third->x,&fourth->x);
|
---|
204 | Log() << Verbose(0) << "z: ",
|
---|
205 | z.Output();
|
---|
206 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
207 | y.MakeNormalVector(&x,&z);
|
---|
208 | Log() << Verbose(0) << "y: ",
|
---|
209 | y.Output();
|
---|
210 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
211 |
|
---|
212 | // rotate vector around first angle
|
---|
213 | first->x.CopyVector(&x);
|
---|
214 | first->x.RotateVector(&z,b - M_PI);
|
---|
215 | Log() << Verbose(0) << "Rotated vector: ",
|
---|
216 | first->x.Output();
|
---|
217 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
218 | // remove the projection onto the rotation plane of the second angle
|
---|
219 | n.CopyVector(&y);
|
---|
220 | n.Scale(first->x.ScalarProduct(&y));
|
---|
221 | Log() << Verbose(0) << "N1: ",
|
---|
222 | n.Output();
|
---|
223 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
224 | first->x.SubtractVector(&n);
|
---|
225 | Log() << Verbose(0) << "Subtracted vector: ",
|
---|
226 | first->x.Output();
|
---|
227 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
228 | n.CopyVector(&z);
|
---|
229 | n.Scale(first->x.ScalarProduct(&z));
|
---|
230 | Log() << Verbose(0) << "N2: ",
|
---|
231 | n.Output();
|
---|
232 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
233 | first->x.SubtractVector(&n);
|
---|
234 | Log() << Verbose(0) << "2nd subtracted vector: ",
|
---|
235 | first->x.Output();
|
---|
236 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
237 |
|
---|
238 | // rotate another vector around second angle
|
---|
239 | n.CopyVector(&y);
|
---|
240 | n.RotateVector(&x,c - M_PI);
|
---|
241 | Log() << Verbose(0) << "2nd Rotated vector: ",
|
---|
242 | n.Output();
|
---|
243 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
244 |
|
---|
245 | // add the two linear independent vectors
|
---|
246 | first->x.AddVector(&n);
|
---|
247 | first->x.Normalize();
|
---|
248 | first->x.Scale(a);
|
---|
249 | first->x.AddVector(&second->x);
|
---|
250 |
|
---|
251 | DoLog(0) && (Log() << Verbose(0) << "resulting coordinates: ");
|
---|
252 | first->x.Output();
|
---|
253 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
254 | } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
|
---|
255 | first->type = periode->AskElement(); // give type
|
---|
256 | mol->AddAtom(first); // add to molecule
|
---|
257 | break;
|
---|
258 |
|
---|
259 | case 'e': // least square distance position to a set of atoms
|
---|
260 | first = new atom;
|
---|
261 | atoms = new (Vector*[128]);
|
---|
262 | valid = true;
|
---|
263 | for(int i=128;i--;)
|
---|
264 | atoms[i] = NULL;
|
---|
265 | int i=0, j=0;
|
---|
266 | cout << Verbose(0) << "Now we need at least three molecules.\n";
|
---|
267 | do {
|
---|
268 | cout << Verbose(0) << "Enter " << i+1 << "th atom: ";
|
---|
269 | cin >> j;
|
---|
270 | if (j != -1) {
|
---|
271 | second = mol->FindAtom(j);
|
---|
272 | atoms[i++] = &(second->x);
|
---|
273 | }
|
---|
274 | } while ((j != -1) && (i<128));
|
---|
275 | if (i >= 2) {
|
---|
276 | first->x.LSQdistance((const Vector **)atoms, i);
|
---|
277 | first->x.Output();
|
---|
278 | first->type = periode->AskElement(); // give type
|
---|
279 | mol->AddAtom(first); // add to molecule
|
---|
280 | } else {
|
---|
281 | delete first;
|
---|
282 | cout << Verbose(0) << "Please enter at least two vectors!\n";
|
---|
283 | }
|
---|
284 | break;
|
---|
285 | };
|
---|
286 | };
|
---|
287 |
|
---|
288 | /** Submenu for centering the atoms in the molecule.
|
---|
289 | * \param *mol molecule with all the atoms
|
---|
290 | */
|
---|
291 | static void CenterAtoms(molecule *mol)
|
---|
292 | {
|
---|
293 | Vector x, y, helper;
|
---|
294 | char choice; // menu choice char
|
---|
295 |
|
---|
296 | cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
|
---|
297 | cout << Verbose(0) << " a - on origin" << endl;
|
---|
298 | cout << Verbose(0) << " b - on center of gravity" << endl;
|
---|
299 | cout << Verbose(0) << " c - within box with additional boundary" << endl;
|
---|
300 | cout << Verbose(0) << " d - within given simulation box" << endl;
|
---|
301 | cout << Verbose(0) << "all else - go back" << endl;
|
---|
302 | cout << Verbose(0) << "===============================================" << endl;
|
---|
303 | cout << Verbose(0) << "INPUT: ";
|
---|
304 | cin >> choice;
|
---|
305 |
|
---|
306 | switch (choice) {
|
---|
307 | default:
|
---|
308 | cout << Verbose(0) << "Not a valid choice." << endl;
|
---|
309 | break;
|
---|
310 | case 'a':
|
---|
311 | cout << Verbose(0) << "Centering atoms in config file on origin." << endl;
|
---|
312 | mol->CenterOrigin();
|
---|
313 | break;
|
---|
314 | case 'b':
|
---|
315 | cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
|
---|
316 | mol->CenterPeriodic();
|
---|
317 | break;
|
---|
318 | case 'c':
|
---|
319 | cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
|
---|
320 | for (int i=0;i<NDIM;i++) {
|
---|
321 | cout << Verbose(0) << "Enter axis " << i << " boundary: ";
|
---|
322 | cin >> y.x[i];
|
---|
323 | }
|
---|
324 | mol->CenterEdge(&x); // make every coordinate positive
|
---|
325 | mol->Center.AddVector(&y); // translate by boundary
|
---|
326 | helper.CopyVector(&y);
|
---|
327 | helper.Scale(2.);
|
---|
328 | helper.AddVector(&x);
|
---|
329 | mol->SetBoxDimension(&helper); // update Box of atoms by boundary
|
---|
330 | break;
|
---|
331 | case 'd':
|
---|
332 | cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
|
---|
333 | for (int i=0;i<NDIM;i++) {
|
---|
334 | cout << Verbose(0) << "Enter axis " << i << " boundary: ";
|
---|
335 | cin >> x.x[i];
|
---|
336 | }
|
---|
337 | // update Box of atoms by boundary
|
---|
338 | mol->SetBoxDimension(&x);
|
---|
339 | // center
|
---|
340 | mol->CenterInBox();
|
---|
341 | break;
|
---|
342 | }
|
---|
343 | };
|
---|
344 |
|
---|
345 | /** Submenu for aligning the atoms in the molecule.
|
---|
346 | * \param *periode periodentafel
|
---|
347 | * \param *mol molecule with all the atoms
|
---|
348 | */
|
---|
349 | static void AlignAtoms(periodentafel *periode, molecule *mol)
|
---|
350 | {
|
---|
351 | atom *first, *second, *third;
|
---|
352 | Vector x,n;
|
---|
353 | char choice; // menu choice char
|
---|
354 |
|
---|
355 | cout << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
|
---|
356 | cout << Verbose(0) << " a - state three atoms defining align plane" << endl;
|
---|
357 | cout << Verbose(0) << " b - state alignment vector" << endl;
|
---|
358 | cout << Verbose(0) << " c - state two atoms in alignment direction" << endl;
|
---|
359 | cout << Verbose(0) << " d - align automatically by least square fit" << endl;
|
---|
360 | cout << Verbose(0) << "all else - go back" << endl;
|
---|
361 | cout << Verbose(0) << "===============================================" << endl;
|
---|
362 | cout << Verbose(0) << "INPUT: ";
|
---|
363 | cin >> choice;
|
---|
364 |
|
---|
365 | switch (choice) {
|
---|
366 | default:
|
---|
367 | case 'a': // three atoms defining mirror plane
|
---|
368 | first = mol->AskAtom("Enter first atom: ");
|
---|
369 | second = mol->AskAtom("Enter second atom: ");
|
---|
370 | third = mol->AskAtom("Enter third atom: ");
|
---|
371 |
|
---|
372 | n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
|
---|
373 | break;
|
---|
374 | case 'b': // normal vector of mirror plane
|
---|
375 | cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
|
---|
376 | n.AskPosition(World::getInstance().getDomain(),0);
|
---|
377 | n.Normalize();
|
---|
378 | break;
|
---|
379 | case 'c': // three atoms defining mirror plane
|
---|
380 | first = mol->AskAtom("Enter first atom: ");
|
---|
381 | second = mol->AskAtom("Enter second atom: ");
|
---|
382 |
|
---|
383 | n.CopyVector((const Vector *)&first->x);
|
---|
384 | n.SubtractVector((const Vector *)&second->x);
|
---|
385 | n.Normalize();
|
---|
386 | break;
|
---|
387 | case 'd':
|
---|
388 | char shorthand[4];
|
---|
389 | Vector a;
|
---|
390 | struct lsq_params param;
|
---|
391 | do {
|
---|
392 | fprintf(stdout, "Enter the element of atoms to be chosen: ");
|
---|
393 | fscanf(stdin, "%3s", shorthand);
|
---|
394 | } while ((param.type = periode->FindElement(shorthand)) == NULL);
|
---|
395 | cout << Verbose(0) << "Element is " << param.type->name << endl;
|
---|
396 | mol->GetAlignvector(¶m);
|
---|
397 | for (int i=NDIM;i--;) {
|
---|
398 | x.x[i] = gsl_vector_get(param.x,i);
|
---|
399 | n.x[i] = gsl_vector_get(param.x,i+NDIM);
|
---|
400 | }
|
---|
401 | gsl_vector_free(param.x);
|
---|
402 | cout << Verbose(0) << "Offset vector: ";
|
---|
403 | x.Output();
|
---|
404 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
405 | n.Normalize();
|
---|
406 | break;
|
---|
407 | };
|
---|
408 | DoLog(0) && (Log() << Verbose(0) << "Alignment vector: ");
|
---|
409 | n.Output();
|
---|
410 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
411 | mol->Align(&n);
|
---|
412 | };
|
---|
413 |
|
---|
414 | /** Submenu for mirroring the atoms in the molecule.
|
---|
415 | * \param *mol molecule with all the atoms
|
---|
416 | */
|
---|
417 | static void MirrorAtoms(molecule *mol)
|
---|
418 | {
|
---|
419 | atom *first, *second, *third;
|
---|
420 | Vector n;
|
---|
421 | char choice; // menu choice char
|
---|
422 |
|
---|
423 | DoLog(0) && (Log() << Verbose(0) << "===========MIRROR ATOMS=========================" << endl);
|
---|
424 | DoLog(0) && (Log() << Verbose(0) << " a - state three atoms defining mirror plane" << endl);
|
---|
425 | DoLog(0) && (Log() << Verbose(0) << " b - state normal vector of mirror plane" << endl);
|
---|
426 | DoLog(0) && (Log() << Verbose(0) << " c - state two atoms in normal direction" << endl);
|
---|
427 | DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
|
---|
428 | DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
|
---|
429 | DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
|
---|
430 | cin >> choice;
|
---|
431 |
|
---|
432 | switch (choice) {
|
---|
433 | default:
|
---|
434 | case 'a': // three atoms defining mirror plane
|
---|
435 | first = mol->AskAtom("Enter first atom: ");
|
---|
436 | second = mol->AskAtom("Enter second atom: ");
|
---|
437 | third = mol->AskAtom("Enter third atom: ");
|
---|
438 |
|
---|
439 | n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
|
---|
440 | break;
|
---|
441 | case 'b': // normal vector of mirror plane
|
---|
442 | DoLog(0) && (Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl);
|
---|
443 | n.AskPosition(World::getInstance().getDomain(),0);
|
---|
444 | n.Normalize();
|
---|
445 | break;
|
---|
446 | case 'c': // three atoms defining mirror plane
|
---|
447 | first = mol->AskAtom("Enter first atom: ");
|
---|
448 | second = mol->AskAtom("Enter second atom: ");
|
---|
449 |
|
---|
450 | n.CopyVector((const Vector *)&first->x);
|
---|
451 | n.SubtractVector((const Vector *)&second->x);
|
---|
452 | n.Normalize();
|
---|
453 | break;
|
---|
454 | };
|
---|
455 | DoLog(0) && (Log() << Verbose(0) << "Normal vector: ");
|
---|
456 | n.Output();
|
---|
457 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
458 | mol->Mirror((const Vector *)&n);
|
---|
459 | };
|
---|
460 |
|
---|
461 | /** Submenu for removing the atoms from the molecule.
|
---|
462 | * \param *mol molecule with all the atoms
|
---|
463 | */
|
---|
464 | static void RemoveAtoms(molecule *mol)
|
---|
465 | {
|
---|
466 | atom *first, *second;
|
---|
467 | int axis;
|
---|
468 | double tmp1, tmp2;
|
---|
469 | char choice; // menu choice char
|
---|
470 |
|
---|
471 | DoLog(0) && (Log() << Verbose(0) << "===========REMOVE ATOMS=========================" << endl);
|
---|
472 | DoLog(0) && (Log() << Verbose(0) << " a - state atom for removal by number" << endl);
|
---|
473 | DoLog(0) && (Log() << Verbose(0) << " b - keep only in radius around atom" << endl);
|
---|
474 | DoLog(0) && (Log() << Verbose(0) << " c - remove this with one axis greater value" << endl);
|
---|
475 | DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
|
---|
476 | DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
|
---|
477 | DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
|
---|
478 | cin >> choice;
|
---|
479 |
|
---|
480 | switch (choice) {
|
---|
481 | default:
|
---|
482 | case 'a':
|
---|
483 | if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: ")))
|
---|
484 | DoLog(1) && (Log() << Verbose(1) << "Atom removed." << endl);
|
---|
485 | else
|
---|
486 | DoLog(1) && (Log() << Verbose(1) << "Atom not found." << endl);
|
---|
487 | break;
|
---|
488 | case 'b':
|
---|
489 | second = mol->AskAtom("Enter number of atom as reference point: ");
|
---|
490 | DoLog(0) && (Log() << Verbose(0) << "Enter radius: ");
|
---|
491 | cin >> tmp1;
|
---|
492 | first = mol->start;
|
---|
493 | second = first->next;
|
---|
494 | while(second != mol->end) {
|
---|
495 | first = second;
|
---|
496 | second = first->next;
|
---|
497 | if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
|
---|
498 | mol->RemoveAtom(first);
|
---|
499 | }
|
---|
500 | break;
|
---|
501 | case 'c':
|
---|
502 | DoLog(0) && (Log() << Verbose(0) << "Which axis is it: ");
|
---|
503 | cin >> axis;
|
---|
504 | DoLog(0) && (Log() << Verbose(0) << "Lower boundary: ");
|
---|
505 | cin >> tmp1;
|
---|
506 | DoLog(0) && (Log() << Verbose(0) << "Upper boundary: ");
|
---|
507 | cin >> tmp2;
|
---|
508 | first = mol->start;
|
---|
509 | second = first->next;
|
---|
510 | while(second != mol->end) {
|
---|
511 | first = second;
|
---|
512 | second = first->next;
|
---|
513 | if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ...
|
---|
514 | //Log() << Verbose(0) << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
|
---|
515 | mol->RemoveAtom(first);
|
---|
516 | }
|
---|
517 | }
|
---|
518 | break;
|
---|
519 | };
|
---|
520 | //mol->Output();
|
---|
521 | choice = 'r';
|
---|
522 | };
|
---|
523 |
|
---|
524 | /** Submenu for measuring out the atoms in the molecule.
|
---|
525 | * \param *periode periodentafel
|
---|
526 | * \param *mol molecule with all the atoms
|
---|
527 | */
|
---|
528 | static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
|
---|
529 | {
|
---|
530 | atom *first, *second, *third;
|
---|
531 | Vector x,y;
|
---|
532 | double min[256], tmp1, tmp2, tmp3;
|
---|
533 | int Z;
|
---|
534 | char choice; // menu choice char
|
---|
535 |
|
---|
536 | DoLog(0) && (Log() << Verbose(0) << "===========MEASURE ATOMS=========================" << endl);
|
---|
537 | DoLog(0) && (Log() << Verbose(0) << " a - calculate bond length between one atom and all others" << endl);
|
---|
538 | DoLog(0) && (Log() << Verbose(0) << " b - calculate bond length between two atoms" << endl);
|
---|
539 | DoLog(0) && (Log() << Verbose(0) << " c - calculate bond angle" << endl);
|
---|
540 | DoLog(0) && (Log() << Verbose(0) << " d - calculate principal axis of the system" << endl);
|
---|
541 | DoLog(0) && (Log() << Verbose(0) << " e - calculate volume of the convex envelope" << endl);
|
---|
542 | DoLog(0) && (Log() << Verbose(0) << " f - calculate temperature from current velocity" << endl);
|
---|
543 | DoLog(0) && (Log() << Verbose(0) << " g - output all temperatures per step from velocities" << endl);
|
---|
544 | DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
|
---|
545 | DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
|
---|
546 | DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
|
---|
547 | cin >> choice;
|
---|
548 |
|
---|
549 | switch(choice) {
|
---|
550 | default:
|
---|
551 | DoLog(1) && (Log() << Verbose(1) << "Not a valid choice." << endl);
|
---|
552 | break;
|
---|
553 | case 'a':
|
---|
554 | first = mol->AskAtom("Enter first atom: ");
|
---|
555 | for (int i=MAX_ELEMENTS;i--;)
|
---|
556 | min[i] = 0.;
|
---|
557 |
|
---|
558 | second = mol->start;
|
---|
559 | while ((second->next != mol->end)) {
|
---|
560 | second = second->next; // advance
|
---|
561 | Z = second->type->Z;
|
---|
562 | tmp1 = 0.;
|
---|
563 | if (first != second) {
|
---|
564 | x.CopyVector((const Vector *)&first->x);
|
---|
565 | x.SubtractVector((const Vector *)&second->x);
|
---|
566 | tmp1 = x.Norm();
|
---|
567 | }
|
---|
568 | if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
|
---|
569 | //Log() << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
|
---|
570 | }
|
---|
571 | for (int i=MAX_ELEMENTS;i--;)
|
---|
572 | if (min[i] != 0.) Log() << 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;
|
---|
573 | break;
|
---|
574 |
|
---|
575 | case 'b':
|
---|
576 | first = mol->AskAtom("Enter first atom: ");
|
---|
577 | second = mol->AskAtom("Enter second atom: ");
|
---|
578 | for (int i=NDIM;i--;)
|
---|
579 | min[i] = 0.;
|
---|
580 | x.CopyVector((const Vector *)&first->x);
|
---|
581 | x.SubtractVector((const Vector *)&second->x);
|
---|
582 | tmp1 = x.Norm();
|
---|
583 | DoLog(1) && (Log() << Verbose(1) << "Distance vector is ");
|
---|
584 | x.Output();
|
---|
585 | DoLog(0) && (Log() << Verbose(0) << "." << endl << "Norm of distance is " << tmp1 << "." << endl);
|
---|
586 | break;
|
---|
587 |
|
---|
588 | case 'c':
|
---|
589 | DoLog(0) && (Log() << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl);
|
---|
590 | first = mol->AskAtom("Enter first atom: ");
|
---|
591 | second = mol->AskAtom("Enter central atom: ");
|
---|
592 | third = mol->AskAtom("Enter last atom: ");
|
---|
593 | tmp1 = tmp2 = tmp3 = 0.;
|
---|
594 | x.CopyVector((const Vector *)&first->x);
|
---|
595 | x.SubtractVector((const Vector *)&second->x);
|
---|
596 | y.CopyVector((const Vector *)&third->x);
|
---|
597 | y.SubtractVector((const Vector *)&second->x);
|
---|
598 | DoLog(0) && (Log() << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ");
|
---|
599 | DoLog(0) && (Log() << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl);
|
---|
600 | break;
|
---|
601 | case 'd':
|
---|
602 | DoLog(0) && (Log() << Verbose(0) << "Evaluating prinicipal axis." << endl);
|
---|
603 | DoLog(0) && (Log() << Verbose(0) << "Shall we rotate? [0/1]: ");
|
---|
604 | cin >> Z;
|
---|
605 | if ((Z >=0) && (Z <=1))
|
---|
606 | mol->PrincipalAxisSystem((bool)Z);
|
---|
607 | else
|
---|
608 | mol->PrincipalAxisSystem(false);
|
---|
609 | break;
|
---|
610 | case 'e':
|
---|
611 | {
|
---|
612 | DoLog(0) && (Log() << Verbose(0) << "Evaluating volume of the convex envelope.");
|
---|
613 | class Tesselation *TesselStruct = NULL;
|
---|
614 | const LinkedCell *LCList = NULL;
|
---|
615 | LCList = new LinkedCell(mol, 10.);
|
---|
616 | FindConvexBorder(mol, TesselStruct, LCList, NULL);
|
---|
617 | double clustervolume = VolumeOfConvexEnvelope(TesselStruct, configuration);
|
---|
618 | DoLog(0) && (Log() << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl);\
|
---|
619 | delete(LCList);
|
---|
620 | delete(TesselStruct);
|
---|
621 | }
|
---|
622 | break;
|
---|
623 | case 'f':
|
---|
624 | mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps);
|
---|
625 | break;
|
---|
626 | case 'g':
|
---|
627 | {
|
---|
628 | char filename[255];
|
---|
629 | DoLog(0) && (Log() << Verbose(0) << "Please enter filename: " << endl);
|
---|
630 | cin >> filename;
|
---|
631 | DoLog(1) && (Log() << Verbose(1) << "Storing temperatures in " << filename << "." << endl);
|
---|
632 | ofstream *output = new ofstream(filename, ios::trunc);
|
---|
633 | if (!mol->OutputTemperatureFromTrajectories(output, 0, mol->MDSteps))
|
---|
634 | DoLog(2) && (Log() << Verbose(2) << "File could not be written." << endl);
|
---|
635 | else
|
---|
636 | DoLog(2) && (Log() << Verbose(2) << "File stored." << endl);
|
---|
637 | output->close();
|
---|
638 | delete(output);
|
---|
639 | }
|
---|
640 | break;
|
---|
641 | }
|
---|
642 | };
|
---|
643 |
|
---|
644 | /** Submenu for measuring out the atoms in the molecule.
|
---|
645 | * \param *mol molecule with all the atoms
|
---|
646 | * \param *configuration configuration structure for the to be written config files of all fragments
|
---|
647 | */
|
---|
648 | static void FragmentAtoms(molecule *mol, config *configuration)
|
---|
649 | {
|
---|
650 | int Order1;
|
---|
651 | clock_t start, end;
|
---|
652 |
|
---|
653 | DoLog(0) && (Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl);
|
---|
654 | DoLog(0) && (Log() << Verbose(0) << "What's the desired bond order: ");
|
---|
655 | cin >> Order1;
|
---|
656 | if (mol->first->next != mol->last) { // there are bonds
|
---|
657 | start = clock();
|
---|
658 | mol->FragmentMolecule(Order1, configuration);
|
---|
659 | end = clock();
|
---|
660 | DoLog(0) && (Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl);
|
---|
661 | } else
|
---|
662 | DoLog(0) && (Log() << Verbose(0) << "Connection matrix has not yet been generated!" << endl);
|
---|
663 | };
|
---|
664 |
|
---|
665 | /********************************************** Submenu routine **************************************/
|
---|
666 |
|
---|
667 | /** Submenu for manipulating atoms.
|
---|
668 | * \param *periode periodentafel
|
---|
669 | * \param *molecules list of molecules whose atoms are to be manipulated
|
---|
670 | */
|
---|
671 | static void ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
|
---|
672 | {
|
---|
673 | atom *first, *second, *third;
|
---|
674 | molecule *mol = NULL;
|
---|
675 | Vector x,y,z,n; // coordinates for absolute point in cell volume
|
---|
676 | double *factor; // unit factor if desired
|
---|
677 | double bond, minBond;
|
---|
678 | char choice; // menu choice char
|
---|
679 | bool valid;
|
---|
680 |
|
---|
681 | DoLog(0) && (Log() << Verbose(0) << "=========MANIPULATE ATOMS======================" << endl);
|
---|
682 | DoLog(0) && (Log() << Verbose(0) << "a - add an atom" << endl);
|
---|
683 | DoLog(0) && (Log() << Verbose(0) << "r - remove an atom" << endl);
|
---|
684 | DoLog(0) && (Log() << Verbose(0) << "b - scale a bond between atoms" << endl);
|
---|
685 | DoLog(0) && (Log() << Verbose(0) << "t - turn an atom round another bond" << endl);
|
---|
686 | DoLog(0) && (Log() << Verbose(0) << "u - change an atoms element" << endl);
|
---|
687 | DoLog(0) && (Log() << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl);
|
---|
688 | DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
|
---|
689 | DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
|
---|
690 | if (molecules->NumberOfActiveMolecules() > 1)
|
---|
691 | DoeLog(2) && (eLog()<< Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl);
|
---|
692 | DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
|
---|
693 | cin >> choice;
|
---|
694 |
|
---|
695 | switch (choice) {
|
---|
696 | default:
|
---|
697 | DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl);
|
---|
698 | break;
|
---|
699 |
|
---|
700 | case 'a': // add atom
|
---|
701 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
702 | if ((*ListRunner)->ActiveFlag) {
|
---|
703 | mol = *ListRunner;
|
---|
704 | DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
|
---|
705 | AddAtoms(periode, mol);
|
---|
706 | }
|
---|
707 | break;
|
---|
708 |
|
---|
709 | case 'b': // scale a bond
|
---|
710 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
711 | if ((*ListRunner)->ActiveFlag) {
|
---|
712 | mol = *ListRunner;
|
---|
713 | DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
|
---|
714 | DoLog(0) && (Log() << Verbose(0) << "Scaling bond length between two atoms." << endl);
|
---|
715 | first = mol->AskAtom("Enter first (fixed) atom: ");
|
---|
716 | second = mol->AskAtom("Enter second (shifting) atom: ");
|
---|
717 | minBond = 0.;
|
---|
718 | for (int i=NDIM;i--;)
|
---|
719 | minBond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
|
---|
720 | minBond = sqrt(minBond);
|
---|
721 | DoLog(0) && (Log() << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << minBond << " a.u." << endl);
|
---|
722 | DoLog(0) && (Log() << Verbose(0) << "Enter new bond length [a.u.]: ");
|
---|
723 | cin >> bond;
|
---|
724 | for (int i=NDIM;i--;) {
|
---|
725 | second->x.x[i] -= (second->x.x[i]-first->x.x[i])/minBond*(minBond-bond);
|
---|
726 | }
|
---|
727 | //Log() << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
|
---|
728 | //second->Output(second->type->No, 1);
|
---|
729 | }
|
---|
730 | break;
|
---|
731 |
|
---|
732 | case 'c': // unit scaling of the metric
|
---|
733 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
734 | if ((*ListRunner)->ActiveFlag) {
|
---|
735 | mol = *ListRunner;
|
---|
736 | DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
|
---|
737 | DoLog(0) && (Log() << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl);
|
---|
738 | DoLog(0) && (Log() << Verbose(0) << "Enter three factors: ");
|
---|
739 | factor = new double[NDIM];
|
---|
740 | cin >> factor[0];
|
---|
741 | cin >> factor[1];
|
---|
742 | cin >> factor[2];
|
---|
743 | valid = true;
|
---|
744 | mol->Scale((const double ** const)&factor);
|
---|
745 | delete[](factor);
|
---|
746 | }
|
---|
747 | break;
|
---|
748 |
|
---|
749 | case 'l': // measure distances or angles
|
---|
750 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
751 | if ((*ListRunner)->ActiveFlag) {
|
---|
752 | mol = *ListRunner;
|
---|
753 | DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
|
---|
754 | MeasureAtoms(periode, mol, configuration);
|
---|
755 | }
|
---|
756 | break;
|
---|
757 |
|
---|
758 | case 'r': // remove atom
|
---|
759 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
760 | if ((*ListRunner)->ActiveFlag) {
|
---|
761 | mol = *ListRunner;
|
---|
762 | DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
|
---|
763 | RemoveAtoms(mol);
|
---|
764 | }
|
---|
765 | break;
|
---|
766 |
|
---|
767 | case 't': // turn/rotate atom
|
---|
768 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
769 | if ((*ListRunner)->ActiveFlag) {
|
---|
770 | mol = *ListRunner;
|
---|
771 | DoLog(0) && (Log() << Verbose(0) << "Turning atom around another bond - first is atom to turn, second (central) and third specify bond" << endl);
|
---|
772 | first = mol->AskAtom("Enter turning atom: ");
|
---|
773 | second = mol->AskAtom("Enter central atom: ");
|
---|
774 | third = mol->AskAtom("Enter bond atom: ");
|
---|
775 | cout << Verbose(0) << "Enter new angle in degrees: ";
|
---|
776 | double tmp = 0.;
|
---|
777 | cin >> tmp;
|
---|
778 | // calculate old angle
|
---|
779 | x.CopyVector((const Vector *)&first->x);
|
---|
780 | x.SubtractVector((const Vector *)&second->x);
|
---|
781 | y.CopyVector((const Vector *)&third->x);
|
---|
782 | y.SubtractVector((const Vector *)&second->x);
|
---|
783 | double alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
|
---|
784 | cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
|
---|
785 | cout << Verbose(0) << alpha << " degrees" << endl;
|
---|
786 | // rotate
|
---|
787 | z.MakeNormalVector(&x,&y);
|
---|
788 | x.RotateVector(&z,(alpha-tmp)*M_PI/180.);
|
---|
789 | x.AddVector(&second->x);
|
---|
790 | first->x.CopyVector(&x);
|
---|
791 | // check new angle
|
---|
792 | x.CopyVector((const Vector *)&first->x);
|
---|
793 | x.SubtractVector((const Vector *)&second->x);
|
---|
794 | alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
|
---|
795 | cout << Verbose(0) << "new Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
|
---|
796 | cout << Verbose(0) << alpha << " degrees" << endl;
|
---|
797 | }
|
---|
798 | break;
|
---|
799 |
|
---|
800 | case 'u': // change an atom's element
|
---|
801 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
802 | if ((*ListRunner)->ActiveFlag) {
|
---|
803 | int Z;
|
---|
804 | mol = *ListRunner;
|
---|
805 | DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
|
---|
806 | first = NULL;
|
---|
807 | do {
|
---|
808 | DoLog(0) && (Log() << Verbose(0) << "Change the element of which atom: ");
|
---|
809 | cin >> Z;
|
---|
810 | } while ((first = mol->FindAtom(Z)) == NULL);
|
---|
811 | DoLog(0) && (Log() << Verbose(0) << "New element by atomic number Z: ");
|
---|
812 | cin >> Z;
|
---|
813 | first->type = periode->FindElement(Z);
|
---|
814 | DoLog(0) && (Log() << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl);
|
---|
815 | }
|
---|
816 | break;
|
---|
817 | }
|
---|
818 | };
|
---|
819 |
|
---|
820 | /** Submenu for manipulating molecules.
|
---|
821 | * \param *periode periodentafel
|
---|
822 | * \param *molecules list of molecule to manipulate
|
---|
823 | */
|
---|
824 | static void ManipulateMolecules(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
|
---|
825 | {
|
---|
826 | atom *first = NULL;
|
---|
827 | Vector x,y,z,n; // coordinates for absolute point in cell volume
|
---|
828 | int j, axis, count, faktor;
|
---|
829 | char choice; // menu choice char
|
---|
830 | molecule *mol = NULL;
|
---|
831 | element **Elements;
|
---|
832 | Vector **vectors;
|
---|
833 | MoleculeLeafClass *Subgraphs = NULL;
|
---|
834 |
|
---|
835 | DoLog(0) && (Log() << Verbose(0) << "=========MANIPULATE GLOBALLY===================" << endl);
|
---|
836 | DoLog(0) && (Log() << Verbose(0) << "c - scale by unit transformation" << endl);
|
---|
837 | DoLog(0) && (Log() << Verbose(0) << "d - duplicate molecule/periodic cell" << endl);
|
---|
838 | DoLog(0) && (Log() << Verbose(0) << "f - fragment molecule many-body bond order style" << endl);
|
---|
839 | DoLog(0) && (Log() << Verbose(0) << "g - center atoms in box" << endl);
|
---|
840 | DoLog(0) && (Log() << Verbose(0) << "i - realign molecule" << endl);
|
---|
841 | DoLog(0) && (Log() << Verbose(0) << "m - mirror all molecules" << endl);
|
---|
842 | DoLog(0) && (Log() << Verbose(0) << "o - create connection matrix" << endl);
|
---|
843 | DoLog(0) && (Log() << Verbose(0) << "t - translate molecule by vector" << endl);
|
---|
844 | DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
|
---|
845 | DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
|
---|
846 | if (molecules->NumberOfActiveMolecules() > 1)
|
---|
847 | DoeLog(2) && (eLog()<< Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl);
|
---|
848 | DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
|
---|
849 | cin >> choice;
|
---|
850 |
|
---|
851 | switch (choice) {
|
---|
852 | default:
|
---|
853 | DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl);
|
---|
854 | break;
|
---|
855 |
|
---|
856 | case 'd': // duplicate the periodic cell along a given axis, given times
|
---|
857 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
858 | if ((*ListRunner)->ActiveFlag) {
|
---|
859 | mol = *ListRunner;
|
---|
860 | DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
|
---|
861 | DoLog(0) && (Log() << Verbose(0) << "State the axis [(+-)123]: ");
|
---|
862 | cin >> axis;
|
---|
863 | DoLog(0) && (Log() << Verbose(0) << "State the factor: ");
|
---|
864 | cin >> faktor;
|
---|
865 |
|
---|
866 | mol->CountAtoms(); // recount atoms
|
---|
867 | if (mol->getAtomCount() != 0) { // if there is more than none
|
---|
868 | count = mol->getAtomCount(); // is changed becausing of adding, thus has to be stored away beforehand
|
---|
869 | Elements = new element *[count];
|
---|
870 | vectors = new Vector *[count];
|
---|
871 | j = 0;
|
---|
872 | first = mol->start;
|
---|
873 | while (first->next != mol->end) { // make a list of all atoms with coordinates and element
|
---|
874 | first = first->next;
|
---|
875 | Elements[j] = first->type;
|
---|
876 | vectors[j] = &first->x;
|
---|
877 | j++;
|
---|
878 | }
|
---|
879 | if (count != j)
|
---|
880 | DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl);
|
---|
881 | x.Zero();
|
---|
882 | y.Zero();
|
---|
883 | y.x[abs(axis)-1] = World::getInstance().getDomain()[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
|
---|
884 | for (int i=1;i<faktor;i++) { // then add this list with respective translation factor times
|
---|
885 | x.AddVector(&y); // per factor one cell width further
|
---|
886 | for (int k=count;k--;) { // go through every atom of the original cell
|
---|
887 | first = new atom(); // create a new body
|
---|
888 | first->x.CopyVector(vectors[k]); // use coordinate of original atom
|
---|
889 | first->x.AddVector(&x); // translate the coordinates
|
---|
890 | first->type = Elements[k]; // insert original element
|
---|
891 | mol->AddAtom(first); // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
|
---|
892 | }
|
---|
893 | }
|
---|
894 | if (mol->first->next != mol->last) // if connect matrix is present already, redo it
|
---|
895 | mol->CreateAdjacencyList(mol->BondDistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
|
---|
896 | // free memory
|
---|
897 | delete[](Elements);
|
---|
898 | delete[](vectors);
|
---|
899 | // correct cell size
|
---|
900 | if (axis < 0) { // if sign was negative, we have to translate everything
|
---|
901 | x.Zero();
|
---|
902 | x.AddVector(&y);
|
---|
903 | x.Scale(-(faktor-1));
|
---|
904 | mol->Translate(&x);
|
---|
905 | }
|
---|
906 | World::getInstance().getDomain()[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
|
---|
907 | }
|
---|
908 | }
|
---|
909 | break;
|
---|
910 |
|
---|
911 | case 'f':
|
---|
912 | FragmentAtoms(mol, configuration);
|
---|
913 | break;
|
---|
914 |
|
---|
915 | case 'g': // center the atoms
|
---|
916 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
917 | if ((*ListRunner)->ActiveFlag) {
|
---|
918 | mol = *ListRunner;
|
---|
919 | DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
|
---|
920 | CenterAtoms(mol);
|
---|
921 | }
|
---|
922 | break;
|
---|
923 |
|
---|
924 | case 'i': // align all atoms
|
---|
925 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
926 | if ((*ListRunner)->ActiveFlag) {
|
---|
927 | mol = *ListRunner;
|
---|
928 | DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
|
---|
929 | AlignAtoms(periode, mol);
|
---|
930 | }
|
---|
931 | break;
|
---|
932 |
|
---|
933 | case 'm': // mirror atoms along a given axis
|
---|
934 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
935 | if ((*ListRunner)->ActiveFlag) {
|
---|
936 | mol = *ListRunner;
|
---|
937 | DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
|
---|
938 | MirrorAtoms(mol);
|
---|
939 | }
|
---|
940 | break;
|
---|
941 |
|
---|
942 | case 'o': // create the connection matrix
|
---|
943 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
944 | if ((*ListRunner)->ActiveFlag) {
|
---|
945 | mol = *ListRunner;
|
---|
946 | double bonddistance;
|
---|
947 | clock_t start,end;
|
---|
948 | DoLog(0) && (Log() << Verbose(0) << "What's the maximum bond distance: ");
|
---|
949 | cin >> bonddistance;
|
---|
950 | start = clock();
|
---|
951 | mol->CreateAdjacencyList(bonddistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
|
---|
952 | end = clock();
|
---|
953 | DoLog(0) && (Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl);
|
---|
954 | }
|
---|
955 | break;
|
---|
956 |
|
---|
957 | case 't': // translate all atoms
|
---|
958 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
959 | if ((*ListRunner)->ActiveFlag) {
|
---|
960 | mol = *ListRunner;
|
---|
961 | DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
|
---|
962 | DoLog(0) && (Log() << Verbose(0) << "Enter translation vector." << endl);
|
---|
963 | x.AskPosition(World::getInstance().getDomain(),0);
|
---|
964 | mol->Center.AddVector((const Vector *)&x);
|
---|
965 | }
|
---|
966 | break;
|
---|
967 | }
|
---|
968 | // Free all
|
---|
969 | if (Subgraphs != NULL) { // free disconnected subgraph list of DFS analysis was performed
|
---|
970 | while (Subgraphs->next != NULL) {
|
---|
971 | Subgraphs = Subgraphs->next;
|
---|
972 | delete(Subgraphs->previous);
|
---|
973 | }
|
---|
974 | delete(Subgraphs);
|
---|
975 | }
|
---|
976 | };
|
---|
977 |
|
---|
978 |
|
---|
979 | /** Submenu for creating new molecules.
|
---|
980 | * \param *periode periodentafel
|
---|
981 | * \param *molecules list of molecules to add to
|
---|
982 | */
|
---|
983 | static void EditMolecules(periodentafel *periode, MoleculeListClass *molecules)
|
---|
984 | {
|
---|
985 | char choice; // menu choice char
|
---|
986 | Vector center;
|
---|
987 | int nr, count;
|
---|
988 | molecule *mol = NULL;
|
---|
989 |
|
---|
990 | DoLog(0) && (Log() << Verbose(0) << "==========EDIT MOLECULES=====================" << endl);
|
---|
991 | DoLog(0) && (Log() << Verbose(0) << "c - create new molecule" << endl);
|
---|
992 | DoLog(0) && (Log() << Verbose(0) << "l - load molecule from xyz file" << endl);
|
---|
993 | DoLog(0) && (Log() << Verbose(0) << "n - change molecule's name" << endl);
|
---|
994 | DoLog(0) && (Log() << Verbose(0) << "N - give molecules filename" << endl);
|
---|
995 | DoLog(0) && (Log() << Verbose(0) << "p - parse atoms in xyz file into molecule" << endl);
|
---|
996 | DoLog(0) && (Log() << Verbose(0) << "r - remove a molecule" << endl);
|
---|
997 | DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
|
---|
998 | DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
|
---|
999 | DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
|
---|
1000 | cin >> choice;
|
---|
1001 |
|
---|
1002 | switch (choice) {
|
---|
1003 | default:
|
---|
1004 | DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl);
|
---|
1005 | break;
|
---|
1006 | case 'c':
|
---|
1007 | mol = World::getInstance().createMolecule();
|
---|
1008 | molecules->insert(mol);
|
---|
1009 | break;
|
---|
1010 |
|
---|
1011 | case 'l': // load from XYZ file
|
---|
1012 | {
|
---|
1013 | char filename[MAXSTRINGSIZE];
|
---|
1014 | DoLog(0) && (Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl);
|
---|
1015 | mol = World::getInstance().createMolecule();
|
---|
1016 | do {
|
---|
1017 | DoLog(0) && (Log() << Verbose(0) << "Enter file name: ");
|
---|
1018 | cin >> filename;
|
---|
1019 | } while (!mol->AddXYZFile(filename));
|
---|
1020 | mol->SetNameFromFilename(filename);
|
---|
1021 | // center at set box dimensions
|
---|
1022 | mol->CenterEdge(¢er);
|
---|
1023 | double * const cell_size = World::getInstance().getDomain();
|
---|
1024 | cell_size[0] = center.x[0];
|
---|
1025 | cell_size[1] = 0;
|
---|
1026 | cell_size[2] = center.x[1];
|
---|
1027 | cell_size[3] = 0;
|
---|
1028 | cell_size[4] = 0;
|
---|
1029 | cell_size[5] = center.x[2];
|
---|
1030 | molecules->insert(mol);
|
---|
1031 | }
|
---|
1032 | break;
|
---|
1033 |
|
---|
1034 | case 'n':
|
---|
1035 | {
|
---|
1036 | char filename[MAXSTRINGSIZE];
|
---|
1037 | do {
|
---|
1038 | DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: ");
|
---|
1039 | cin >> nr;
|
---|
1040 | mol = molecules->ReturnIndex(nr);
|
---|
1041 | } while (mol == NULL);
|
---|
1042 | DoLog(0) && (Log() << Verbose(0) << "Enter name: ");
|
---|
1043 | cin >> filename;
|
---|
1044 | strcpy(mol->name, filename);
|
---|
1045 | }
|
---|
1046 | break;
|
---|
1047 |
|
---|
1048 | case 'N':
|
---|
1049 | {
|
---|
1050 | char filename[MAXSTRINGSIZE];
|
---|
1051 | do {
|
---|
1052 | DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: ");
|
---|
1053 | cin >> nr;
|
---|
1054 | mol = molecules->ReturnIndex(nr);
|
---|
1055 | } while (mol == NULL);
|
---|
1056 | DoLog(0) && (Log() << Verbose(0) << "Enter name: ");
|
---|
1057 | cin >> filename;
|
---|
1058 | mol->SetNameFromFilename(filename);
|
---|
1059 | }
|
---|
1060 | break;
|
---|
1061 |
|
---|
1062 | case 'p': // parse XYZ file
|
---|
1063 | {
|
---|
1064 | char filename[MAXSTRINGSIZE];
|
---|
1065 | mol = NULL;
|
---|
1066 | do {
|
---|
1067 | DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: ");
|
---|
1068 | cin >> nr;
|
---|
1069 | mol = molecules->ReturnIndex(nr);
|
---|
1070 | } while (mol == NULL);
|
---|
1071 | DoLog(0) && (Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl);
|
---|
1072 | do {
|
---|
1073 | DoLog(0) && (Log() << Verbose(0) << "Enter file name: ");
|
---|
1074 | cin >> filename;
|
---|
1075 | } while (!mol->AddXYZFile(filename));
|
---|
1076 | mol->SetNameFromFilename(filename);
|
---|
1077 | }
|
---|
1078 | break;
|
---|
1079 |
|
---|
1080 | case 'r':
|
---|
1081 | DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: ");
|
---|
1082 | cin >> nr;
|
---|
1083 | count = 1;
|
---|
1084 | for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
1085 | if (nr == (*ListRunner)->IndexNr) {
|
---|
1086 | mol = *ListRunner;
|
---|
1087 | molecules->ListOfMolecules.erase(ListRunner);
|
---|
1088 | delete(mol);
|
---|
1089 | break;
|
---|
1090 | }
|
---|
1091 | break;
|
---|
1092 | }
|
---|
1093 | };
|
---|
1094 |
|
---|
1095 |
|
---|
1096 | /** Submenu for merging molecules.
|
---|
1097 | * \param *periode periodentafel
|
---|
1098 | * \param *molecules list of molecules to add to
|
---|
1099 | */
|
---|
1100 | static void MergeMolecules(periodentafel *periode, MoleculeListClass *molecules)
|
---|
1101 | {
|
---|
1102 | char choice; // menu choice char
|
---|
1103 |
|
---|
1104 | DoLog(0) && (Log() << Verbose(0) << "===========MERGE MOLECULES=====================" << endl);
|
---|
1105 | DoLog(0) && (Log() << Verbose(0) << "a - simple add of one molecule to another" << endl);
|
---|
1106 | DoLog(0) && (Log() << Verbose(0) << "b - count the number of bonds of two elements" << endl);
|
---|
1107 | DoLog(0) && (Log() << Verbose(0) << "B - count the number of bonds of three elements " << endl);
|
---|
1108 | DoLog(0) && (Log() << Verbose(0) << "e - embedding merge of two molecules" << endl);
|
---|
1109 | DoLog(0) && (Log() << Verbose(0) << "h - count the number of hydrogen bonds" << endl);
|
---|
1110 | DoLog(0) && (Log() << Verbose(0) << "b - count the number of hydrogen bonds" << endl);
|
---|
1111 | DoLog(0) && (Log() << Verbose(0) << "m - multi-merge of all molecules" << endl);
|
---|
1112 | DoLog(0) && (Log() << Verbose(0) << "s - scatter merge of two molecules" << endl);
|
---|
1113 | DoLog(0) && (Log() << Verbose(0) << "t - simple merge of two molecules" << endl);
|
---|
1114 | DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
|
---|
1115 | DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
|
---|
1116 | DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
|
---|
1117 | cin >> choice;
|
---|
1118 |
|
---|
1119 | switch (choice) {
|
---|
1120 | default:
|
---|
1121 | DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl);
|
---|
1122 | break;
|
---|
1123 |
|
---|
1124 | case 'a':
|
---|
1125 | {
|
---|
1126 | int src, dest;
|
---|
1127 | molecule *srcmol = NULL, *destmol = NULL;
|
---|
1128 | {
|
---|
1129 | do {
|
---|
1130 | DoLog(0) && (Log() << Verbose(0) << "Enter index of destination molecule: ");
|
---|
1131 | cin >> dest;
|
---|
1132 | destmol = molecules->ReturnIndex(dest);
|
---|
1133 | } while ((destmol == NULL) && (dest != -1));
|
---|
1134 | do {
|
---|
1135 | DoLog(0) && (Log() << Verbose(0) << "Enter index of source molecule to add from: ");
|
---|
1136 | cin >> src;
|
---|
1137 | srcmol = molecules->ReturnIndex(src);
|
---|
1138 | } while ((srcmol == NULL) && (src != -1));
|
---|
1139 | if ((src != -1) && (dest != -1))
|
---|
1140 | molecules->SimpleAdd(srcmol, destmol);
|
---|
1141 | }
|
---|
1142 | }
|
---|
1143 | break;
|
---|
1144 |
|
---|
1145 | case 'b':
|
---|
1146 | {
|
---|
1147 | const int nr = 2;
|
---|
1148 | char *names[nr] = {"first", "second"};
|
---|
1149 | int Z[nr];
|
---|
1150 | element *elements[nr];
|
---|
1151 | for (int i=0;i<nr;i++) {
|
---|
1152 | Z[i] = 0;
|
---|
1153 | do {
|
---|
1154 | cout << "Enter " << names[i] << " element: ";
|
---|
1155 | cin >> Z[i];
|
---|
1156 | } while ((Z[i] <= 0) && (Z[i] > MAX_ELEMENTS));
|
---|
1157 | elements[i] = periode->FindElement(Z[i]);
|
---|
1158 | }
|
---|
1159 | const int count = CountBondsOfTwo(molecules, elements[0], elements[1]);
|
---|
1160 | cout << endl << "There are " << count << " ";
|
---|
1161 | for (int i=0;i<nr;i++) {
|
---|
1162 | if (i==0)
|
---|
1163 | cout << elements[i]->symbol;
|
---|
1164 | else
|
---|
1165 | cout << "-" << elements[i]->symbol;
|
---|
1166 | }
|
---|
1167 | cout << " bonds." << endl;
|
---|
1168 | }
|
---|
1169 | break;
|
---|
1170 |
|
---|
1171 | case 'B':
|
---|
1172 | {
|
---|
1173 | const int nr = 3;
|
---|
1174 | char *names[nr] = {"first", "second", "third"};
|
---|
1175 | int Z[nr];
|
---|
1176 | element *elements[nr];
|
---|
1177 | for (int i=0;i<nr;i++) {
|
---|
1178 | Z[i] = 0;
|
---|
1179 | do {
|
---|
1180 | cout << "Enter " << names[i] << " element: ";
|
---|
1181 | cin >> Z[i];
|
---|
1182 | } while ((Z[i] <= 0) && (Z[i] > MAX_ELEMENTS));
|
---|
1183 | elements[i] = periode->FindElement(Z[i]);
|
---|
1184 | }
|
---|
1185 | const int count = CountBondsOfThree(molecules, elements[0], elements[1], elements[2]);
|
---|
1186 | cout << endl << "There are " << count << " ";
|
---|
1187 | for (int i=0;i<nr;i++) {
|
---|
1188 | if (i==0)
|
---|
1189 | cout << elements[i]->symbol;
|
---|
1190 | else
|
---|
1191 | cout << "-" << elements[i]->symbol;
|
---|
1192 | }
|
---|
1193 | cout << " bonds." << endl;
|
---|
1194 | }
|
---|
1195 | break;
|
---|
1196 |
|
---|
1197 | case 'e':
|
---|
1198 | {
|
---|
1199 | int src, dest;
|
---|
1200 | molecule *srcmol = NULL, *destmol = NULL;
|
---|
1201 | do {
|
---|
1202 | DoLog(0) && (Log() << Verbose(0) << "Enter index of matrix molecule (the variable one): ");
|
---|
1203 | cin >> src;
|
---|
1204 | srcmol = molecules->ReturnIndex(src);
|
---|
1205 | } while ((srcmol == NULL) && (src != -1));
|
---|
1206 | do {
|
---|
1207 | DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule to merge into (the fixed one): ");
|
---|
1208 | cin >> dest;
|
---|
1209 | destmol = molecules->ReturnIndex(dest);
|
---|
1210 | } while ((destmol == NULL) && (dest != -1));
|
---|
1211 | if ((src != -1) && (dest != -1))
|
---|
1212 | molecules->EmbedMerge(destmol, srcmol);
|
---|
1213 | }
|
---|
1214 | break;
|
---|
1215 |
|
---|
1216 | case 'h':
|
---|
1217 | {
|
---|
1218 | int Z;
|
---|
1219 | cout << "Please enter interface element: ";
|
---|
1220 | cin >> Z;
|
---|
1221 | element * const InterfaceElement = periode->FindElement(Z);
|
---|
1222 | cout << endl << "There are " << CountHydrogenBridgeBonds(molecules, InterfaceElement) << " hydrogen bridges with connections to " << (InterfaceElement != 0 ? InterfaceElement->name : "None") << "." << endl;
|
---|
1223 | }
|
---|
1224 | break;
|
---|
1225 |
|
---|
1226 | case 'm':
|
---|
1227 | {
|
---|
1228 | int nr;
|
---|
1229 | molecule *mol = NULL;
|
---|
1230 | do {
|
---|
1231 | DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule to merge into: ");
|
---|
1232 | cin >> nr;
|
---|
1233 | mol = molecules->ReturnIndex(nr);
|
---|
1234 | } while ((mol == NULL) && (nr != -1));
|
---|
1235 | if (nr != -1) {
|
---|
1236 | int N = molecules->ListOfMolecules.size()-1;
|
---|
1237 | int *src = new int(N);
|
---|
1238 | for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
1239 | if ((*ListRunner)->IndexNr != nr)
|
---|
1240 | src[N++] = (*ListRunner)->IndexNr;
|
---|
1241 | molecules->SimpleMultiMerge(mol, src, N);
|
---|
1242 | delete[](src);
|
---|
1243 | }
|
---|
1244 | }
|
---|
1245 | break;
|
---|
1246 |
|
---|
1247 | case 's':
|
---|
1248 | DoLog(0) && (Log() << Verbose(0) << "Not implemented yet." << endl);
|
---|
1249 | break;
|
---|
1250 |
|
---|
1251 | case 't':
|
---|
1252 | {
|
---|
1253 | int src, dest;
|
---|
1254 | molecule *srcmol = NULL, *destmol = NULL;
|
---|
1255 | {
|
---|
1256 | do {
|
---|
1257 | DoLog(0) && (Log() << Verbose(0) << "Enter index of destination molecule: ");
|
---|
1258 | cin >> dest;
|
---|
1259 | destmol = molecules->ReturnIndex(dest);
|
---|
1260 | } while ((destmol == NULL) && (dest != -1));
|
---|
1261 | do {
|
---|
1262 | DoLog(0) && (Log() << Verbose(0) << "Enter index of source molecule to merge into: ");
|
---|
1263 | cin >> src;
|
---|
1264 | srcmol = molecules->ReturnIndex(src);
|
---|
1265 | } while ((srcmol == NULL) && (src != -1));
|
---|
1266 | if ((src != -1) && (dest != -1))
|
---|
1267 | molecules->SimpleMerge(srcmol, destmol);
|
---|
1268 | }
|
---|
1269 | }
|
---|
1270 | break;
|
---|
1271 | }
|
---|
1272 | };
|
---|
1273 |
|
---|
1274 | /********************************************** Test routine **************************************/
|
---|
1275 |
|
---|
1276 | /** Is called always as option 'T' in the menu.
|
---|
1277 | * \param *molecules list of molecules
|
---|
1278 | */
|
---|
1279 | static void testroutine(MoleculeListClass *molecules)
|
---|
1280 | {
|
---|
1281 | // the current test routine checks the functionality of the KeySet&Graph concept:
|
---|
1282 | // We want to have a multiindex (the KeySet) describing a unique subgraph
|
---|
1283 | int i, comp, counter=0;
|
---|
1284 |
|
---|
1285 | // create a clone
|
---|
1286 | molecule *mol = NULL;
|
---|
1287 | if (molecules->ListOfMolecules.size() != 0) // clone
|
---|
1288 | mol = (molecules->ListOfMolecules.front())->CopyMolecule();
|
---|
1289 | else {
|
---|
1290 | DoeLog(0) && (eLog()<< Verbose(0) << "I don't have anything to test on ... ");
|
---|
1291 | performCriticalExit();
|
---|
1292 | return;
|
---|
1293 | }
|
---|
1294 | atom *Walker = mol->start;
|
---|
1295 |
|
---|
1296 | // generate some KeySets
|
---|
1297 | DoLog(0) && (Log() << Verbose(0) << "Generating KeySets." << endl);
|
---|
1298 | KeySet TestSets[mol->getAtomCount()+1];
|
---|
1299 | i=1;
|
---|
1300 | while (Walker->next != mol->end) {
|
---|
1301 | Walker = Walker->next;
|
---|
1302 | for (int j=0;j<i;j++) {
|
---|
1303 | TestSets[j].insert(Walker->nr);
|
---|
1304 | }
|
---|
1305 | i++;
|
---|
1306 | }
|
---|
1307 | DoLog(0) && (Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl);
|
---|
1308 | KeySetTestPair test;
|
---|
1309 | test = TestSets[mol->getAtomCount()-1].insert(Walker->nr);
|
---|
1310 | if (test.second) {
|
---|
1311 | DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl);
|
---|
1312 | } else {
|
---|
1313 | DoLog(1) && (Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl);
|
---|
1314 | }
|
---|
1315 | TestSets[mol->getAtomCount()].insert(mol->end->previous->nr);
|
---|
1316 | TestSets[mol->getAtomCount()].insert(mol->end->previous->previous->previous->nr);
|
---|
1317 |
|
---|
1318 | // constructing Graph structure
|
---|
1319 | DoLog(0) && (Log() << Verbose(0) << "Generating Subgraph class." << endl);
|
---|
1320 | Graph Subgraphs;
|
---|
1321 |
|
---|
1322 | // insert KeySets into Subgraphs
|
---|
1323 | DoLog(0) && (Log() << Verbose(0) << "Inserting KeySets into Subgraph class." << endl);
|
---|
1324 | for (int j=0;j<mol->getAtomCount();j++) {
|
---|
1325 | Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
|
---|
1326 | }
|
---|
1327 | DoLog(0) && (Log() << Verbose(0) << "Testing insertion of already present item in Subgraph." << endl);
|
---|
1328 | GraphTestPair test2;
|
---|
1329 | test2 = Subgraphs.insert(GraphPair (TestSets[mol->getAtomCount()],pair<int, double>(counter++, 1.)));
|
---|
1330 | if (test2.second) {
|
---|
1331 | DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl);
|
---|
1332 | } else {
|
---|
1333 | DoLog(1) && (Log() << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl);
|
---|
1334 | }
|
---|
1335 |
|
---|
1336 | // show graphs
|
---|
1337 | DoLog(0) && (Log() << Verbose(0) << "Showing Subgraph's contents, checking that it's sorted." << endl);
|
---|
1338 | Graph::iterator A = Subgraphs.begin();
|
---|
1339 | while (A != Subgraphs.end()) {
|
---|
1340 | DoLog(0) && (Log() << Verbose(0) << (*A).second.first << ": ");
|
---|
1341 | KeySet::iterator key = (*A).first.begin();
|
---|
1342 | comp = -1;
|
---|
1343 | while (key != (*A).first.end()) {
|
---|
1344 | if ((*key) > comp)
|
---|
1345 | DoLog(0) && (Log() << Verbose(0) << (*key) << " ");
|
---|
1346 | else
|
---|
1347 | DoLog(0) && (Log() << Verbose(0) << (*key) << "! ");
|
---|
1348 | comp = (*key);
|
---|
1349 | key++;
|
---|
1350 | }
|
---|
1351 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
1352 | A++;
|
---|
1353 | }
|
---|
1354 | delete(mol);
|
---|
1355 | };
|
---|
1356 |
|
---|
1357 | #endif
|
---|
1358 |
|
---|
1359 | /** Tries given filename or standard on saving the config file.
|
---|
1360 | * \param *ConfigFileName name of file
|
---|
1361 | * \param *configuration pointer to configuration structure with all the values
|
---|
1362 | * \param *periode pointer to periodentafel structure with all the elements
|
---|
1363 | * \param *molecules list of molecules structure with all the atoms and coordinates
|
---|
1364 | */
|
---|
1365 | static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, MoleculeListClass *molecules)
|
---|
1366 | {
|
---|
1367 | char filename[MAXSTRINGSIZE];
|
---|
1368 | ofstream output;
|
---|
1369 | molecule *mol = World::getInstance().createMolecule();
|
---|
1370 | mol->SetNameFromFilename(ConfigFileName);
|
---|
1371 |
|
---|
1372 | if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
|
---|
1373 | DoeLog(2) && (eLog()<< Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl);
|
---|
1374 | }
|
---|
1375 |
|
---|
1376 |
|
---|
1377 | // first save as PDB data
|
---|
1378 | if (ConfigFileName != NULL)
|
---|
1379 | strcpy(filename, ConfigFileName);
|
---|
1380 | if (output == NULL)
|
---|
1381 | strcpy(filename,"main_pcp_linux");
|
---|
1382 | DoLog(0) && (Log() << Verbose(0) << "Saving as pdb input ");
|
---|
1383 | if (configuration->SavePDB(filename, molecules))
|
---|
1384 | DoLog(0) && (Log() << Verbose(0) << "done." << endl);
|
---|
1385 | else
|
---|
1386 | DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
|
---|
1387 |
|
---|
1388 | // then save as tremolo data file
|
---|
1389 | if (ConfigFileName != NULL)
|
---|
1390 | strcpy(filename, ConfigFileName);
|
---|
1391 | if (output == NULL)
|
---|
1392 | strcpy(filename,"main_pcp_linux");
|
---|
1393 | DoLog(0) && (Log() << Verbose(0) << "Saving as tremolo data input ");
|
---|
1394 | if (configuration->SaveTREMOLO(filename, molecules))
|
---|
1395 | DoLog(0) && (Log() << Verbose(0) << "done." << endl);
|
---|
1396 | else
|
---|
1397 | DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
|
---|
1398 |
|
---|
1399 | // translate each to its center and merge all molecules in MoleculeListClass into this molecule
|
---|
1400 | int N = molecules->ListOfMolecules.size();
|
---|
1401 | int *src = new int[N];
|
---|
1402 | N=0;
|
---|
1403 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
|
---|
1404 | src[N++] = (*ListRunner)->IndexNr;
|
---|
1405 | (*ListRunner)->Translate(&(*ListRunner)->Center);
|
---|
1406 | }
|
---|
1407 | molecules->SimpleMultiAdd(mol, src, N);
|
---|
1408 | delete[](src);
|
---|
1409 |
|
---|
1410 | // ... and translate back
|
---|
1411 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
|
---|
1412 | (*ListRunner)->Center.Scale(-1.);
|
---|
1413 | (*ListRunner)->Translate(&(*ListRunner)->Center);
|
---|
1414 | (*ListRunner)->Center.Scale(-1.);
|
---|
1415 | }
|
---|
1416 |
|
---|
1417 | DoLog(0) && (Log() << Verbose(0) << "Storing configuration ... " << endl);
|
---|
1418 | // get correct valence orbitals
|
---|
1419 | mol->CalculateOrbitals(*configuration);
|
---|
1420 | configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble;
|
---|
1421 | if (ConfigFileName != NULL) { // test the file name
|
---|
1422 | strcpy(filename, ConfigFileName);
|
---|
1423 | output.open(filename, ios::trunc);
|
---|
1424 | } else if (strlen(configuration->configname) != 0) {
|
---|
1425 | strcpy(filename, configuration->configname);
|
---|
1426 | output.open(configuration->configname, ios::trunc);
|
---|
1427 | } else {
|
---|
1428 | strcpy(filename, DEFAULTCONFIG);
|
---|
1429 | output.open(DEFAULTCONFIG, ios::trunc);
|
---|
1430 | }
|
---|
1431 | output.close();
|
---|
1432 | output.clear();
|
---|
1433 | DoLog(0) && (Log() << Verbose(0) << "Saving of config file ");
|
---|
1434 | if (configuration->Save(filename, periode, mol))
|
---|
1435 | DoLog(0) && (Log() << Verbose(0) << "successful." << endl);
|
---|
1436 | else
|
---|
1437 | DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
|
---|
1438 |
|
---|
1439 | // and save to xyz file
|
---|
1440 | if (ConfigFileName != NULL) {
|
---|
1441 | strcpy(filename, ConfigFileName);
|
---|
1442 | strcat(filename, ".xyz");
|
---|
1443 | output.open(filename, ios::trunc);
|
---|
1444 | }
|
---|
1445 | if (output == NULL) {
|
---|
1446 | strcpy(filename,"main_pcp_linux");
|
---|
1447 | strcat(filename, ".xyz");
|
---|
1448 | output.open(filename, ios::trunc);
|
---|
1449 | }
|
---|
1450 | DoLog(0) && (Log() << Verbose(0) << "Saving of XYZ file ");
|
---|
1451 | if (mol->MDSteps <= 1) {
|
---|
1452 | if (mol->OutputXYZ(&output))
|
---|
1453 | DoLog(0) && (Log() << Verbose(0) << "successful." << endl);
|
---|
1454 | else
|
---|
1455 | DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
|
---|
1456 | } else {
|
---|
1457 | if (mol->OutputTrajectoriesXYZ(&output))
|
---|
1458 | DoLog(0) && (Log() << Verbose(0) << "successful." << endl);
|
---|
1459 | else
|
---|
1460 | DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
|
---|
1461 | }
|
---|
1462 | output.close();
|
---|
1463 | output.clear();
|
---|
1464 |
|
---|
1465 | // and save as MPQC configuration
|
---|
1466 | if (ConfigFileName != NULL)
|
---|
1467 | strcpy(filename, ConfigFileName);
|
---|
1468 | if (output == NULL)
|
---|
1469 | strcpy(filename,"main_pcp_linux");
|
---|
1470 | DoLog(0) && (Log() << Verbose(0) << "Saving as mpqc input ");
|
---|
1471 | if (configuration->SaveMPQC(filename, mol))
|
---|
1472 | DoLog(0) && (Log() << Verbose(0) << "done." << endl);
|
---|
1473 | else
|
---|
1474 | DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
|
---|
1475 |
|
---|
1476 | if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
|
---|
1477 | DoeLog(2) && (eLog()<< Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl);
|
---|
1478 | }
|
---|
1479 |
|
---|
1480 | World::getInstance().destroyMolecule(mol);
|
---|
1481 | };
|
---|
1482 |
|
---|
1483 | /** Parses the command line options.
|
---|
1484 | * Note that this function is from now on transitional. All commands that are not passed
|
---|
1485 | * here are handled by CommandLineParser and the actions of CommandLineUIFactory.
|
---|
1486 | * \param argc argument count
|
---|
1487 | * \param **argv arguments array
|
---|
1488 | * \param *molecules list of molecules structure
|
---|
1489 | * \param *periode elements structure
|
---|
1490 | * \param configuration config file structure
|
---|
1491 | * \param *ConfigFileName pointer to config file name in **argv
|
---|
1492 | * \param *PathToDatabases pointer to db's path in **argv
|
---|
1493 | * \param &ArgcList list of arguments that we do not parse here
|
---|
1494 | * \return exit code (0 - successful, all else - something's wrong)
|
---|
1495 | */
|
---|
1496 | static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode,
|
---|
1497 | config& configuration, char **ConfigFileName, set<int> &ArgcList)
|
---|
1498 | {
|
---|
1499 | Vector x,y,z,n; // coordinates for absolute point in cell volume
|
---|
1500 | double *factor; // unit factor if desired
|
---|
1501 | ifstream test;
|
---|
1502 | ofstream output;
|
---|
1503 | string line;
|
---|
1504 | atom *first;
|
---|
1505 | bool SaveFlag = false;
|
---|
1506 | int ExitFlag = 0;
|
---|
1507 | int j;
|
---|
1508 | double volume = 0.;
|
---|
1509 | enum ConfigStatus configPresent = absent;
|
---|
1510 | clock_t start,end;
|
---|
1511 | double MaxDistance = -1;
|
---|
1512 | int argptr;
|
---|
1513 | molecule *mol = NULL;
|
---|
1514 | string BondGraphFileName("\n");
|
---|
1515 | bool DatabasePathGiven = false;
|
---|
1516 |
|
---|
1517 | if (argc > 1) { // config file specified as option
|
---|
1518 | // 1. : Parse options that just set variables or print help
|
---|
1519 | argptr = 1;
|
---|
1520 | do {
|
---|
1521 | if (argv[argptr][0] == '-') {
|
---|
1522 | DoLog(0) && (Log() << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n");
|
---|
1523 | argptr++;
|
---|
1524 | switch(argv[argptr-1][1]) {
|
---|
1525 | case 'h':
|
---|
1526 | case 'H':
|
---|
1527 | case '?':
|
---|
1528 | ArgcList.insert(argptr-1);
|
---|
1529 | return(1);
|
---|
1530 | break;
|
---|
1531 | case 'v':
|
---|
1532 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
---|
1533 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying verbosity: -v <level>" << endl);
|
---|
1534 | performCriticalExit();
|
---|
1535 | } else {
|
---|
1536 | setVerbosity(atoi(argv[argptr]));
|
---|
1537 | ArgcList.insert(argptr-1);
|
---|
1538 | ArgcList.insert(argptr);
|
---|
1539 | argptr++;
|
---|
1540 | }
|
---|
1541 | break;
|
---|
1542 | case 'V':
|
---|
1543 | ArgcList.insert(argptr-1);
|
---|
1544 | return(1);
|
---|
1545 | break;
|
---|
1546 | case 'B':
|
---|
1547 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1548 | if ((argptr+5 >= argc)) {
|
---|
1549 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for setting Box: -B <xx> <<xy> <<xz> <yy> <yz> <zz>" << endl);
|
---|
1550 | performCriticalExit();
|
---|
1551 | } else {
|
---|
1552 | ArgcList.insert(argptr-1);
|
---|
1553 | ArgcList.insert(argptr);
|
---|
1554 | ArgcList.insert(argptr+1);
|
---|
1555 | ArgcList.insert(argptr+2);
|
---|
1556 | ArgcList.insert(argptr+3);
|
---|
1557 | ArgcList.insert(argptr+4);
|
---|
1558 | ArgcList.insert(argptr+5);
|
---|
1559 | argptr+=6;
|
---|
1560 | }
|
---|
1561 | break;
|
---|
1562 | case 'e':
|
---|
1563 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
---|
1564 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl);
|
---|
1565 | performCriticalExit();
|
---|
1566 | } else {
|
---|
1567 | ArgcList.insert(argptr-1);
|
---|
1568 | ArgcList.insert(argptr);
|
---|
1569 | argptr+=1;
|
---|
1570 | }
|
---|
1571 | break;
|
---|
1572 | case 'g':
|
---|
1573 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
---|
1574 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying bond length table: -g <table file>" << endl);
|
---|
1575 | performCriticalExit();
|
---|
1576 | } else {
|
---|
1577 | ArgcList.insert(argptr-1);
|
---|
1578 | ArgcList.insert(argptr);
|
---|
1579 | argptr+=1;
|
---|
1580 | }
|
---|
1581 | break;
|
---|
1582 | case 'M':
|
---|
1583 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
---|
1584 | ExitFlag = 255;
|
---|
1585 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting MPQC basis: -M <basis name>" << endl);
|
---|
1586 | performCriticalExit();
|
---|
1587 | } else {
|
---|
1588 | ArgcList.insert(argptr-1);
|
---|
1589 | ArgcList.insert(argptr);
|
---|
1590 | argptr+=1;
|
---|
1591 | }
|
---|
1592 | break;
|
---|
1593 | case 'n':
|
---|
1594 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
---|
1595 | ExitFlag = 255;
|
---|
1596 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting fast-parsing: -n <0/1>" << endl);
|
---|
1597 | performCriticalExit();
|
---|
1598 | } else {
|
---|
1599 | ArgcList.insert(argptr-1);
|
---|
1600 | ArgcList.insert(argptr);
|
---|
1601 | argptr+=1;
|
---|
1602 | }
|
---|
1603 | break;
|
---|
1604 | case 'X':
|
---|
1605 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
---|
1606 | ExitFlag = 255;
|
---|
1607 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting default molecule name: -X <name>" << endl);
|
---|
1608 | performCriticalExit();
|
---|
1609 | } else {
|
---|
1610 | ArgcList.insert(argptr-1);
|
---|
1611 | ArgcList.insert(argptr);
|
---|
1612 | argptr+=1;
|
---|
1613 | }
|
---|
1614 | break;
|
---|
1615 | default: // no match? Step on
|
---|
1616 | argptr++;
|
---|
1617 | break;
|
---|
1618 | }
|
---|
1619 | } else
|
---|
1620 | argptr++;
|
---|
1621 | } while (argptr < argc);
|
---|
1622 |
|
---|
1623 | // 3b. Find config file name and parse if possible, also BondGraphFileName
|
---|
1624 | if (argv[1][0] != '-') {
|
---|
1625 | // simply create a new molecule, wherein the config file is loaded and the manipulation takes place
|
---|
1626 | DoLog(0) && (Log() << Verbose(0) << "Config file given." << endl);
|
---|
1627 | test.open(argv[1], ios::in);
|
---|
1628 | if (test == NULL) {
|
---|
1629 | //return (1);
|
---|
1630 | output.open(argv[1], ios::out);
|
---|
1631 | if (output == NULL) {
|
---|
1632 | DoLog(1) && (Log() << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl);
|
---|
1633 | configPresent = absent;
|
---|
1634 | } else {
|
---|
1635 | DoLog(0) && (Log() << Verbose(0) << "Empty configuration file." << endl);
|
---|
1636 | strcpy(*ConfigFileName, argv[1]);
|
---|
1637 | configPresent = empty;
|
---|
1638 | output.close();
|
---|
1639 | }
|
---|
1640 | } else {
|
---|
1641 | test.close();
|
---|
1642 | strcpy(*ConfigFileName, argv[1]);
|
---|
1643 | DoLog(1) && (Log() << Verbose(1) << "Specified config file found, parsing ... ");
|
---|
1644 | switch (configuration.TestSyntax(*ConfigFileName, periode)) {
|
---|
1645 | case 1:
|
---|
1646 | DoLog(0) && (Log() << Verbose(0) << "new syntax." << endl);
|
---|
1647 | configuration.Load(*ConfigFileName, BondGraphFileName, periode, molecules);
|
---|
1648 | configPresent = present;
|
---|
1649 | break;
|
---|
1650 | case 0:
|
---|
1651 | DoLog(0) && (Log() << Verbose(0) << "old syntax." << endl);
|
---|
1652 | configuration.LoadOld(*ConfigFileName, BondGraphFileName, periode, molecules);
|
---|
1653 | configPresent = present;
|
---|
1654 | break;
|
---|
1655 | default:
|
---|
1656 | DoLog(0) && (Log() << Verbose(0) << "Unknown syntax or empty, yet present file." << endl);
|
---|
1657 | configPresent = empty;
|
---|
1658 | }
|
---|
1659 | }
|
---|
1660 | } else
|
---|
1661 | configPresent = absent;
|
---|
1662 | // set mol to first active molecule
|
---|
1663 | if (molecules->ListOfMolecules.size() != 0) {
|
---|
1664 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
---|
1665 | if ((*ListRunner)->ActiveFlag) {
|
---|
1666 | mol = *ListRunner;
|
---|
1667 | break;
|
---|
1668 | }
|
---|
1669 | }
|
---|
1670 | if (mol == NULL) {
|
---|
1671 | mol = World::getInstance().createMolecule();
|
---|
1672 | mol->ActiveFlag = true;
|
---|
1673 | if (*ConfigFileName != NULL)
|
---|
1674 | mol->SetNameFromFilename(*ConfigFileName);
|
---|
1675 | molecules->insert(mol);
|
---|
1676 | }
|
---|
1677 |
|
---|
1678 | // 4. parse again through options, now for those depending on elements db and config presence
|
---|
1679 | argptr = 1;
|
---|
1680 | do {
|
---|
1681 | DoLog(0) && (Log() << Verbose(0) << "Current Command line argument: " << argv[argptr] << "." << endl);
|
---|
1682 | if (argv[argptr][0] == '-') {
|
---|
1683 | argptr++;
|
---|
1684 | if ((configPresent == present) || (configPresent == empty)) {
|
---|
1685 | switch(argv[argptr-1][1]) {
|
---|
1686 | case 'p':
|
---|
1687 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1688 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
---|
1689 | ExitFlag = 255;
|
---|
1690 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough arguments for parsing: -p <xyz file>" << endl);
|
---|
1691 | performCriticalExit();
|
---|
1692 | } else {
|
---|
1693 | SaveFlag = true;
|
---|
1694 | DoLog(1) && (Log() << Verbose(1) << "Parsing xyz file for new atoms." << endl);
|
---|
1695 | if (!mol->AddXYZFile(argv[argptr]))
|
---|
1696 | DoLog(2) && (Log() << Verbose(2) << "File not found." << endl);
|
---|
1697 | else {
|
---|
1698 | DoLog(2) && (Log() << Verbose(2) << "File found and parsed." << endl);
|
---|
1699 | configPresent = present;
|
---|
1700 | }
|
---|
1701 | }
|
---|
1702 | break;
|
---|
1703 | case 'a':
|
---|
1704 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1705 | if ((argptr+4 >= argc) || (argv[argptr][0] == '-')) {
|
---|
1706 | ExitFlag = 255;
|
---|
1707 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough arguments for adding atom: -a <Z> --position <x> <y> <z>" << endl);
|
---|
1708 | performCriticalExit();
|
---|
1709 | } else {
|
---|
1710 | ArgcList.insert(argptr-1);
|
---|
1711 | ArgcList.insert(argptr);
|
---|
1712 | ArgcList.insert(argptr+1);
|
---|
1713 | ArgcList.insert(argptr+2);
|
---|
1714 | ArgcList.insert(argptr+3);
|
---|
1715 | ArgcList.insert(argptr+4);
|
---|
1716 | argptr+=5;
|
---|
1717 | }
|
---|
1718 | break;
|
---|
1719 | default: // no match? Don't step on (this is done in next switch's default)
|
---|
1720 | break;
|
---|
1721 | }
|
---|
1722 | }
|
---|
1723 | if (configPresent == present) {
|
---|
1724 | switch(argv[argptr-1][1]) {
|
---|
1725 | case 'D':
|
---|
1726 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1727 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
---|
1728 | ExitFlag = 255;
|
---|
1729 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for depth-first-search analysis: -D <max. bond distance>" << endl);
|
---|
1730 | performCriticalExit();
|
---|
1731 | } else {
|
---|
1732 | ArgcList.insert(argptr-1);
|
---|
1733 | ArgcList.insert(argptr);
|
---|
1734 | argptr+=1;
|
---|
1735 | }
|
---|
1736 | break;
|
---|
1737 | case 'I':
|
---|
1738 | DoLog(1) && (Log() << Verbose(1) << "Dissecting molecular system into a set of disconnected subgraphs ... " << endl);
|
---|
1739 | ArgcList.insert(argptr-1);
|
---|
1740 | argptr+=0;
|
---|
1741 | break;
|
---|
1742 | case 'C':
|
---|
1743 | {
|
---|
1744 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1745 | if ((argptr+11 >= argc)) {
|
---|
1746 | ExitFlag = 255;
|
---|
1747 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C[p] <type: E/P/S> [more params] <output> <bin output> <BinStart> <BinEnd>" << endl);
|
---|
1748 | performCriticalExit();
|
---|
1749 | } else {
|
---|
1750 | switch(argv[argptr][0]) {
|
---|
1751 | case 'E':
|
---|
1752 | ArgcList.insert(argptr-1);
|
---|
1753 | ArgcList.insert(argptr);
|
---|
1754 | ArgcList.insert(argptr+1);
|
---|
1755 | ArgcList.insert(argptr+2);
|
---|
1756 | ArgcList.insert(argptr+3);
|
---|
1757 | ArgcList.insert(argptr+4);
|
---|
1758 | ArgcList.insert(argptr+5);
|
---|
1759 | ArgcList.insert(argptr+6);
|
---|
1760 | ArgcList.insert(argptr+7);
|
---|
1761 | ArgcList.insert(argptr+8);
|
---|
1762 | ArgcList.insert(argptr+9);
|
---|
1763 | ArgcList.insert(argptr+10);
|
---|
1764 | ArgcList.insert(argptr+11);
|
---|
1765 | argptr+=12;
|
---|
1766 | break;
|
---|
1767 |
|
---|
1768 | case 'P':
|
---|
1769 | ArgcList.insert(argptr-1);
|
---|
1770 | ArgcList.insert(argptr);
|
---|
1771 | ArgcList.insert(argptr+1);
|
---|
1772 | ArgcList.insert(argptr+2);
|
---|
1773 | ArgcList.insert(argptr+3);
|
---|
1774 | ArgcList.insert(argptr+4);
|
---|
1775 | ArgcList.insert(argptr+5);
|
---|
1776 | ArgcList.insert(argptr+6);
|
---|
1777 | ArgcList.insert(argptr+7);
|
---|
1778 | ArgcList.insert(argptr+8);
|
---|
1779 | ArgcList.insert(argptr+9);
|
---|
1780 | ArgcList.insert(argptr+10);
|
---|
1781 | ArgcList.insert(argptr+11);
|
---|
1782 | ArgcList.insert(argptr+12);
|
---|
1783 | ArgcList.insert(argptr+13);
|
---|
1784 | ArgcList.insert(argptr+14);
|
---|
1785 | argptr+=15;
|
---|
1786 | break;
|
---|
1787 |
|
---|
1788 | case 'S':
|
---|
1789 | ArgcList.insert(argptr-1);
|
---|
1790 | ArgcList.insert(argptr);
|
---|
1791 | ArgcList.insert(argptr+1);
|
---|
1792 | ArgcList.insert(argptr+2);
|
---|
1793 | ArgcList.insert(argptr+3);
|
---|
1794 | ArgcList.insert(argptr+4);
|
---|
1795 | ArgcList.insert(argptr+5);
|
---|
1796 | ArgcList.insert(argptr+6);
|
---|
1797 | ArgcList.insert(argptr+7);
|
---|
1798 | ArgcList.insert(argptr+8);
|
---|
1799 | ArgcList.insert(argptr+9);
|
---|
1800 | ArgcList.insert(argptr+10);
|
---|
1801 | ArgcList.insert(argptr+11);
|
---|
1802 | ArgcList.insert(argptr+12);
|
---|
1803 | ArgcList.insert(argptr+13);
|
---|
1804 | ArgcList.insert(argptr+14);
|
---|
1805 | argptr+=15;
|
---|
1806 | break;
|
---|
1807 |
|
---|
1808 | default:
|
---|
1809 | ExitFlag = 255;
|
---|
1810 | DoeLog(0) && (eLog()<< Verbose(0) << "Invalid type given for pair correlation analysis: -C <type: E/P/S> [more params] <output> <bin output>" << endl);
|
---|
1811 | performCriticalExit();
|
---|
1812 | break;
|
---|
1813 | }
|
---|
1814 | }
|
---|
1815 | break;
|
---|
1816 | }
|
---|
1817 | case 'E':
|
---|
1818 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1819 | if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr]))) {
|
---|
1820 | ExitFlag = 255;
|
---|
1821 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for changing element: -E <atom nr.> --element <Z>" << endl);
|
---|
1822 | performCriticalExit();
|
---|
1823 | } else {
|
---|
1824 | ArgcList.insert(argptr-1);
|
---|
1825 | ArgcList.insert(argptr);
|
---|
1826 | ArgcList.insert(argptr+1);
|
---|
1827 | ArgcList.insert(argptr+2);
|
---|
1828 | argptr+=3;
|
---|
1829 | }
|
---|
1830 | break;
|
---|
1831 | case 'F':
|
---|
1832 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1833 | if ((argptr+12 >= argc) || (argv[argptr][0] == '-')) {
|
---|
1834 | ExitFlag = 255;
|
---|
1835 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for filling with molecule: -F <filler xyz file> --MaxDistance <distance or -1> --distances <x> <y> <z> --lengths <surface> <randatm> <randmol> --DoRotate <0/1>" << endl);
|
---|
1836 | performCriticalExit();
|
---|
1837 | } else {
|
---|
1838 | ArgcList.insert(argptr-1);
|
---|
1839 | ArgcList.insert(argptr);
|
---|
1840 | ArgcList.insert(argptr+1);
|
---|
1841 | ArgcList.insert(argptr+2);
|
---|
1842 | ArgcList.insert(argptr+3);
|
---|
1843 | ArgcList.insert(argptr+4);
|
---|
1844 | ArgcList.insert(argptr+5);
|
---|
1845 | ArgcList.insert(argptr+6);
|
---|
1846 | ArgcList.insert(argptr+7);
|
---|
1847 | ArgcList.insert(argptr+8);
|
---|
1848 | ArgcList.insert(argptr+9);
|
---|
1849 | ArgcList.insert(argptr+10);
|
---|
1850 | ArgcList.insert(argptr+11);
|
---|
1851 | ArgcList.insert(argptr+12);
|
---|
1852 | argptr+=13;
|
---|
1853 | }
|
---|
1854 | break;
|
---|
1855 | case 'A':
|
---|
1856 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1857 | if ((argptr+2 >= argc) || (argv[argptr][0] == '-')) {
|
---|
1858 | ExitFlag =255;
|
---|
1859 | DoeLog(0) && (eLog()<< Verbose(0) << "Missing source file for bonds in molecule: -A <bond sourcefile> --molecule-by-id <molecule_id>" << endl);
|
---|
1860 | performCriticalExit();
|
---|
1861 | } else {
|
---|
1862 | ArgcList.insert(argptr-1);
|
---|
1863 | ArgcList.insert(argptr);
|
---|
1864 | ArgcList.insert(argptr+1);
|
---|
1865 | ArgcList.insert(argptr+2);
|
---|
1866 | argptr+=3;
|
---|
1867 | }
|
---|
1868 | break;
|
---|
1869 |
|
---|
1870 | case 'J':
|
---|
1871 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1872 | if ((argptr+2 >= argc) || (argv[argptr][0] == '-')) {
|
---|
1873 | ExitFlag =255;
|
---|
1874 | DoeLog(0) && (eLog()<< Verbose(0) << "Missing path of adjacency file: -J <path> --molecule-by-id <molecule_id>" << endl);
|
---|
1875 | performCriticalExit();
|
---|
1876 | } else {
|
---|
1877 | ArgcList.insert(argptr-1);
|
---|
1878 | ArgcList.insert(argptr);
|
---|
1879 | ArgcList.insert(argptr+1);
|
---|
1880 | ArgcList.insert(argptr+2);
|
---|
1881 | argptr+=3;
|
---|
1882 | }
|
---|
1883 | break;
|
---|
1884 |
|
---|
1885 | case 'j':
|
---|
1886 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1887 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
---|
1888 | ExitFlag =255;
|
---|
1889 | DoeLog(0) && (eLog()<< Verbose(0) << "Missing path of bonds file: -j <path> --molecule-by-id <molecule_id>" << endl);
|
---|
1890 | performCriticalExit();
|
---|
1891 | } else {
|
---|
1892 | ArgcList.insert(argptr-1);
|
---|
1893 | ArgcList.insert(argptr);
|
---|
1894 | ArgcList.insert(argptr+1);
|
---|
1895 | ArgcList.insert(argptr+2);
|
---|
1896 | argptr+=3;
|
---|
1897 | }
|
---|
1898 | break;
|
---|
1899 |
|
---|
1900 | case 'N':
|
---|
1901 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1902 | if ((argptr+4 >= argc) || (argv[argptr][0] == '-')){
|
---|
1903 | ExitFlag = 255;
|
---|
1904 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for non-convex envelope: -N <molecule_id> --sphere-radius <radius> --nonconvex-file <output prefix>" << endl);
|
---|
1905 | performCriticalExit();
|
---|
1906 | } else {
|
---|
1907 | ArgcList.insert(argptr-1);
|
---|
1908 | ArgcList.insert(argptr);
|
---|
1909 | ArgcList.insert(argptr+1);
|
---|
1910 | ArgcList.insert(argptr+2);
|
---|
1911 | ArgcList.insert(argptr+3);
|
---|
1912 | ArgcList.insert(argptr+4);
|
---|
1913 | argptr+=5;
|
---|
1914 | }
|
---|
1915 | break;
|
---|
1916 | case 'S':
|
---|
1917 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1918 | if ((argptr+2 >= argc) || (argv[argptr][0] == '-')) {
|
---|
1919 | ExitFlag = 255;
|
---|
1920 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for storing tempature: -S <temperature file> --molecule-by-id 0" << endl);
|
---|
1921 | performCriticalExit();
|
---|
1922 | } else {
|
---|
1923 | ArgcList.insert(argptr-1);
|
---|
1924 | ArgcList.insert(argptr);
|
---|
1925 | ArgcList.insert(argptr+1);
|
---|
1926 | ArgcList.insert(argptr+2);
|
---|
1927 | argptr+=3;
|
---|
1928 | }
|
---|
1929 | break;
|
---|
1930 | case 'L':
|
---|
1931 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1932 | if ((argptr+8 >= argc) || (argv[argptr][0] == '-')) {
|
---|
1933 | ExitFlag = 255;
|
---|
1934 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for linear interpolation: -L <prefix> --start-step <step0> --end-step <step1> --molecule-by-id 0 --id-mapping <0/1>" << endl);
|
---|
1935 | performCriticalExit();
|
---|
1936 | } else {
|
---|
1937 | ArgcList.insert(argptr-1);
|
---|
1938 | ArgcList.insert(argptr);
|
---|
1939 | ArgcList.insert(argptr+1);
|
---|
1940 | ArgcList.insert(argptr+2);
|
---|
1941 | ArgcList.insert(argptr+3);
|
---|
1942 | ArgcList.insert(argptr+4);
|
---|
1943 | ArgcList.insert(argptr+5);
|
---|
1944 | ArgcList.insert(argptr+6);
|
---|
1945 | ArgcList.insert(argptr+7);
|
---|
1946 | ArgcList.insert(argptr+8);
|
---|
1947 | argptr+=9;
|
---|
1948 | }
|
---|
1949 | break;
|
---|
1950 | case 'P':
|
---|
1951 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1952 | if ((argptr+2 >= argc) || (argv[argptr][0] == '-')) {
|
---|
1953 | ExitFlag = 255;
|
---|
1954 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file> --molecule-by-id <molecule_id>" << endl);
|
---|
1955 | performCriticalExit();
|
---|
1956 | } else {
|
---|
1957 | ArgcList.insert(argptr-1);
|
---|
1958 | ArgcList.insert(argptr);
|
---|
1959 | ArgcList.insert(argptr+1);
|
---|
1960 | ArgcList.insert(argptr+2);
|
---|
1961 | argptr+=3;
|
---|
1962 | }
|
---|
1963 | break;
|
---|
1964 | case 'R':
|
---|
1965 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1966 | if ((argptr+4 >= argc) || (argv[argptr][0] == '-')) {
|
---|
1967 | ExitFlag = 255;
|
---|
1968 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -R <distance> --position <x> <y> <z>" << endl);
|
---|
1969 | performCriticalExit();
|
---|
1970 | } else {
|
---|
1971 | ArgcList.insert(argptr-1);
|
---|
1972 | ArgcList.insert(argptr);
|
---|
1973 | ArgcList.insert(argptr+1);
|
---|
1974 | ArgcList.insert(argptr+2);
|
---|
1975 | ArgcList.insert(argptr+3);
|
---|
1976 | ArgcList.insert(argptr+4);
|
---|
1977 | argptr+=5;
|
---|
1978 | }
|
---|
1979 | break;
|
---|
1980 | case 't':
|
---|
1981 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
1982 | if ((argptr+4 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
|
---|
1983 | ExitFlag = 255;
|
---|
1984 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for translation: -t <x> <y> <z> --molecule-by-id <molecule_id> --periodic <0/1>" << endl);
|
---|
1985 | performCriticalExit();
|
---|
1986 | } else {
|
---|
1987 | ArgcList.insert(argptr-1);
|
---|
1988 | ArgcList.insert(argptr);
|
---|
1989 | ArgcList.insert(argptr+1);
|
---|
1990 | ArgcList.insert(argptr+2);
|
---|
1991 | ArgcList.insert(argptr+3);
|
---|
1992 | ArgcList.insert(argptr+4);
|
---|
1993 | ArgcList.insert(argptr+5);
|
---|
1994 | ArgcList.insert(argptr+6);
|
---|
1995 | argptr+=7;
|
---|
1996 | }
|
---|
1997 | break;
|
---|
1998 | case 's':
|
---|
1999 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
2000 | if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
|
---|
2001 | ExitFlag = 255;
|
---|
2002 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for scaling: -s <factor_x> [factor_y] [factor_z]" << endl);
|
---|
2003 | performCriticalExit();
|
---|
2004 | } else {
|
---|
2005 | ArgcList.insert(argptr-1);
|
---|
2006 | ArgcList.insert(argptr);
|
---|
2007 | ArgcList.insert(argptr+1);
|
---|
2008 | ArgcList.insert(argptr+2);
|
---|
2009 | argptr+=3;
|
---|
2010 | }
|
---|
2011 | break;
|
---|
2012 | case 'b':
|
---|
2013 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
2014 | if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
|
---|
2015 | ExitFlag = 255;
|
---|
2016 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl);
|
---|
2017 | performCriticalExit();
|
---|
2018 | } else {
|
---|
2019 | ArgcList.insert(argptr-1);
|
---|
2020 | ArgcList.insert(argptr);
|
---|
2021 | ArgcList.insert(argptr+1);
|
---|
2022 | ArgcList.insert(argptr+2);
|
---|
2023 | ArgcList.insert(argptr+3);
|
---|
2024 | ArgcList.insert(argptr+4);
|
---|
2025 | ArgcList.insert(argptr+5);
|
---|
2026 | argptr+=6;
|
---|
2027 | }
|
---|
2028 | break;
|
---|
2029 | case 'B':
|
---|
2030 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
2031 | if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
|
---|
2032 | ExitFlag = 255;
|
---|
2033 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl);
|
---|
2034 | performCriticalExit();
|
---|
2035 | } else {
|
---|
2036 | ArgcList.insert(argptr-1);
|
---|
2037 | ArgcList.insert(argptr);
|
---|
2038 | ArgcList.insert(argptr+1);
|
---|
2039 | ArgcList.insert(argptr+2);
|
---|
2040 | ArgcList.insert(argptr+3);
|
---|
2041 | ArgcList.insert(argptr+4);
|
---|
2042 | ArgcList.insert(argptr+5);
|
---|
2043 | argptr+=6;
|
---|
2044 | }
|
---|
2045 | break;
|
---|
2046 | case 'c':
|
---|
2047 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
2048 | if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
|
---|
2049 | ExitFlag = 255;
|
---|
2050 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl);
|
---|
2051 | performCriticalExit();
|
---|
2052 | } else {
|
---|
2053 | ArgcList.insert(argptr-1);
|
---|
2054 | ArgcList.insert(argptr);
|
---|
2055 | ArgcList.insert(argptr+1);
|
---|
2056 | ArgcList.insert(argptr+2);
|
---|
2057 | argptr+=3;
|
---|
2058 | }
|
---|
2059 | break;
|
---|
2060 | case 'O':
|
---|
2061 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
2062 | ArgcList.insert(argptr-1);
|
---|
2063 | argptr+=0;
|
---|
2064 | break;
|
---|
2065 | case 'r':
|
---|
2066 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
2067 | if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr]))) {
|
---|
2068 | ExitFlag = 255;
|
---|
2069 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -r <id>" << endl);
|
---|
2070 | performCriticalExit();
|
---|
2071 | } else {
|
---|
2072 | ArgcList.insert(argptr-1);
|
---|
2073 | ArgcList.insert(argptr);
|
---|
2074 | argptr+=1;
|
---|
2075 | }
|
---|
2076 | break;
|
---|
2077 | case 'f':
|
---|
2078 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
2079 | if ((argptr+1 >= argc) || (argv[argptr][0] == '-')) {
|
---|
2080 | ExitFlag = 255;
|
---|
2081 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl);
|
---|
2082 | performCriticalExit();
|
---|
2083 | } else {
|
---|
2084 | ArgcList.insert(argptr-1);
|
---|
2085 | ArgcList.insert(argptr);
|
---|
2086 | ArgcList.insert(argptr+1);
|
---|
2087 | ArgcList.insert(argptr+2);
|
---|
2088 | ArgcList.insert(argptr+3);
|
---|
2089 | ArgcList.insert(argptr+4);
|
---|
2090 | argptr+=5;
|
---|
2091 | }
|
---|
2092 | break;
|
---|
2093 | case 'm':
|
---|
2094 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
2095 | j = atoi(argv[argptr++]);
|
---|
2096 | if ((j<0) || (j>1)) {
|
---|
2097 | DoeLog(1) && (eLog()<< Verbose(1) << "Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl);
|
---|
2098 | j = 0;
|
---|
2099 | }
|
---|
2100 | if (j) {
|
---|
2101 | SaveFlag = true;
|
---|
2102 | DoLog(0) && (Log() << Verbose(0) << "Converting to prinicipal axis system." << endl);
|
---|
2103 | mol->PrincipalAxisSystem((bool)j);
|
---|
2104 | } else
|
---|
2105 | ArgcList.insert(argptr-1);
|
---|
2106 | argptr+=0;
|
---|
2107 | break;
|
---|
2108 | case 'o':
|
---|
2109 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
2110 | if ((argptr+4 >= argc) || (argv[argptr][0] == '-')){
|
---|
2111 | ExitFlag = 255;
|
---|
2112 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for convex envelope: -o <molecule_id> --output-file <output file> --output-file <binned output file>" << endl);
|
---|
2113 | performCriticalExit();
|
---|
2114 | } else {
|
---|
2115 | ArgcList.insert(argptr-1);
|
---|
2116 | ArgcList.insert(argptr);
|
---|
2117 | ArgcList.insert(argptr+1);
|
---|
2118 | ArgcList.insert(argptr+2);
|
---|
2119 | ArgcList.insert(argptr+3);
|
---|
2120 | ArgcList.insert(argptr+4);
|
---|
2121 | argptr+=5;
|
---|
2122 | }
|
---|
2123 | break;
|
---|
2124 | case 'U':
|
---|
2125 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
2126 | if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
|
---|
2127 | ExitFlag = 255;
|
---|
2128 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl);
|
---|
2129 | performCriticalExit();
|
---|
2130 | } else {
|
---|
2131 | volume = atof(argv[argptr++]);
|
---|
2132 | DoLog(0) && (Log() << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl);
|
---|
2133 | }
|
---|
2134 | case 'u':
|
---|
2135 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
2136 | if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) ) {
|
---|
2137 | if (volume != -1)
|
---|
2138 | ExitFlag = 255;
|
---|
2139 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl);
|
---|
2140 | performCriticalExit();
|
---|
2141 | } else {
|
---|
2142 | ArgcList.insert(argptr-1);
|
---|
2143 | ArgcList.insert(argptr);
|
---|
2144 | argptr+=1;
|
---|
2145 | }
|
---|
2146 | break;
|
---|
2147 | case 'd':
|
---|
2148 | if (ExitFlag == 0) ExitFlag = 1;
|
---|
2149 | if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
|
---|
2150 | ExitFlag = 255;
|
---|
2151 | DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl);
|
---|
2152 | performCriticalExit();
|
---|
2153 | } else {
|
---|
2154 | ArgcList.insert(argptr-1);
|
---|
2155 | ArgcList.insert(argptr);
|
---|
2156 | ArgcList.insert(argptr+1);
|
---|
2157 | ArgcList.insert(argptr+2);
|
---|
2158 | argptr+=3;
|
---|
2159 | }
|
---|
2160 | break;
|
---|
2161 | default: // no match? Step on
|
---|
2162 | if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step!
|
---|
2163 | argptr++;
|
---|
2164 | break;
|
---|
2165 | }
|
---|
2166 | }
|
---|
2167 | } else argptr++;
|
---|
2168 | } while (argptr < argc);
|
---|
2169 | if (SaveFlag)
|
---|
2170 | configuration.SaveAll(*ConfigFileName, periode, molecules);
|
---|
2171 | } else { // no arguments, hence scan the elements db
|
---|
2172 | if (periode->LoadPeriodentafel(configuration.databasepath))
|
---|
2173 | DoLog(0) && (Log() << Verbose(0) << "Element list loaded successfully." << endl);
|
---|
2174 | else
|
---|
2175 | DoLog(0) && (Log() << Verbose(0) << "Element list loading failed." << endl);
|
---|
2176 | configuration.RetrieveConfigPathAndName("main_pcp_linux");
|
---|
2177 | }
|
---|
2178 | return(ExitFlag);
|
---|
2179 | };
|
---|
2180 |
|
---|
2181 | /********************************************** Main routine **************************************/
|
---|
2182 |
|
---|
2183 | void cleanUp(){
|
---|
2184 | World::purgeInstance();
|
---|
2185 | logger::purgeInstance();
|
---|
2186 | errorLogger::purgeInstance();
|
---|
2187 | UIFactory::purgeInstance();
|
---|
2188 | MapOfActions::purgeInstance();
|
---|
2189 | CommandLineParser::purgeInstance();
|
---|
2190 | ActionRegistry::purgeInstance();
|
---|
2191 | ActionHistory::purgeInstance();
|
---|
2192 | Memory::getState();
|
---|
2193 | }
|
---|
2194 |
|
---|
2195 | int main(int argc, char **argv)
|
---|
2196 | {
|
---|
2197 | config *configuration = World::getInstance().getConfig();
|
---|
2198 | // while we are non interactive, we want to abort from asserts
|
---|
2199 | //ASSERT_DO(Assert::Abort);
|
---|
2200 | molecule *mol = NULL;
|
---|
2201 | Vector x, y, z, n;
|
---|
2202 | ifstream test;
|
---|
2203 | ofstream output;
|
---|
2204 | string line;
|
---|
2205 | char **Arguments = NULL;
|
---|
2206 | int ArgcSize = 0;
|
---|
2207 | int ExitFlag = 0;
|
---|
2208 | bool ArgumentsCopied = false;
|
---|
2209 | char *ConfigFileName = new char[MAXSTRINGSIZE];
|
---|
2210 |
|
---|
2211 | // print version check whether arguments are present at all
|
---|
2212 | cout << ESPACKVersion << endl;
|
---|
2213 | if (argc < 2) {
|
---|
2214 | cout << "Obtain help with " << argv[0] << " -h." << endl;
|
---|
2215 | cleanUp();
|
---|
2216 | Memory::getState();
|
---|
2217 | return(1);
|
---|
2218 | }
|
---|
2219 |
|
---|
2220 |
|
---|
2221 | setVerbosity(0);
|
---|
2222 | // need to init the history before any action is created
|
---|
2223 | ActionHistory::init();
|
---|
2224 |
|
---|
2225 | // In the interactive mode, we can leave the user the choice in case of error
|
---|
2226 | ASSERT_DO(Assert::Ask);
|
---|
2227 |
|
---|
2228 | // from this moment on, we need to be sure to deeinitialize in the correct order
|
---|
2229 | // this is handled by the cleanup function
|
---|
2230 | atexit(cleanUp);
|
---|
2231 |
|
---|
2232 | // Parse command line options and if present create respective UI
|
---|
2233 | {
|
---|
2234 | set<int> ArgcList;
|
---|
2235 | ArgcList.insert(0); // push back program!
|
---|
2236 | ArgcList.insert(1); // push back config file name
|
---|
2237 | // handle arguments by ParseCommandLineOptions()
|
---|
2238 | ExitFlag = ParseCommandLineOptions(argc,argv,World::getInstance().getMolecules(),World::getInstance().getPeriode(),*World::getInstance().getConfig(), &ConfigFileName, ArgcList);
|
---|
2239 | World::getInstance().setExitFlag(ExitFlag);
|
---|
2240 | // copy all remaining arguments to a new argv
|
---|
2241 | Arguments = new char *[ArgcList.size()];
|
---|
2242 | cout << "The following arguments are handled by CommandLineParser: ";
|
---|
2243 | for (set<int>::iterator ArgcRunner = ArgcList.begin(); ArgcRunner != ArgcList.end(); ++ArgcRunner) {
|
---|
2244 | Arguments[ArgcSize] = new char[strlen(argv[*ArgcRunner])+2];
|
---|
2245 | strcpy(Arguments[ArgcSize], argv[*ArgcRunner]);
|
---|
2246 | cout << " " << argv[*ArgcRunner];
|
---|
2247 | ArgcSize++;
|
---|
2248 | }
|
---|
2249 | cout << endl;
|
---|
2250 | ArgumentsCopied = true;
|
---|
2251 | // handle remaining arguments by CommandLineParser
|
---|
2252 | MapOfActions::getInstance().AddOptionsToParser();
|
---|
2253 | map <std::string, std::string> ShortFormToActionMap = MapOfActions::getInstance().getShortFormToActionMap();
|
---|
2254 | CommandLineParser::getInstance().Run(ArgcSize,Arguments, ShortFormToActionMap);
|
---|
2255 | if (!CommandLineParser::getInstance().isEmpty()) {
|
---|
2256 | DoLog(0) && (Log() << Verbose(0) << "Setting UI to CommandLine." << endl);
|
---|
2257 | UIFactory::makeUserInterface(UIFactory::CommandLine);
|
---|
2258 | } else {
|
---|
2259 | DoLog(0) && (Log() << Verbose(0) << "Setting UI to Text." << endl);
|
---|
2260 | UIFactory::makeUserInterface(UIFactory::Text);
|
---|
2261 | }
|
---|
2262 | }
|
---|
2263 |
|
---|
2264 | {
|
---|
2265 | MainWindow *mainWindow = UIFactory::getInstance().makeMainWindow();
|
---|
2266 | mainWindow->display();
|
---|
2267 | delete mainWindow;
|
---|
2268 | }
|
---|
2269 |
|
---|
2270 | Log() << Verbose(0) << "Saving to " << ConfigFileName << "." << endl;
|
---|
2271 | World::getInstance().getConfig()->SaveAll(ConfigFileName, World::getInstance().getPeriode(), World::getInstance().getMolecules());
|
---|
2272 |
|
---|
2273 | // free the new argv
|
---|
2274 | if (ArgumentsCopied) {
|
---|
2275 | for (int i=0; i<ArgcSize;i++)
|
---|
2276 | delete[](Arguments[i]);
|
---|
2277 | delete[](Arguments);
|
---|
2278 | }
|
---|
2279 | delete[](ConfigFileName);
|
---|
2280 |
|
---|
2281 | ExitFlag = World::getInstance().getExitFlag();
|
---|
2282 | return (ExitFlag == 1 ? 0 : ExitFlag);
|
---|
2283 | }
|
---|
2284 |
|
---|
2285 | /********************************************** E N D **************************************************/
|
---|