Changes in src/builder.cpp [4fc93f:a67d19]
- File:
-
- 1 edited
-
src/builder.cpp (modified) (134 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/builder.cpp
r4fc93f ra67d19 50 50 using namespace std; 51 51 52 #include <cstring> 53 54 #include "analysis_bonds.hpp" 52 55 #include "analysis_correlation.hpp" 53 56 #include "atom.hpp" … … 65 68 #include "molecule.hpp" 66 69 #include "periodentafel.hpp" 70 #include "version.h" 71 #include "World.hpp" 67 72 68 73 /********************************************* Subsubmenu routine ************************************/ … … 81 86 bool valid; 82 87 83 Log()<< Verbose(0) << "===========ADD ATOM============================" << endl;84 Log()<< Verbose(0) << " a - state absolute coordinates of atom" << endl;85 Log()<< Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;86 Log()<< Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;87 Log()<< Verbose(0) << " d - state two atoms, two angles and a distance" << endl;88 Log()<< Verbose(0) << " e - least square distance position to a set of atoms" << endl;89 Log()<< Verbose(0) << "all else - go back" << endl;90 Log()<< Verbose(0) << "===============================================" << endl;91 Log()<< Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;92 Log()<< Verbose(0) << "INPUT: ";88 cout << Verbose(0) << "===========ADD ATOM============================" << endl; 89 cout << Verbose(0) << " a - state absolute coordinates of atom" << endl; 90 cout << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl; 91 cout << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl; 92 cout << Verbose(0) << " d - state two atoms, two angles and a distance" << endl; 93 cout << Verbose(0) << " e - least square distance position to a set of atoms" << endl; 94 cout << Verbose(0) << "all else - go back" << endl; 95 cout << Verbose(0) << "===============================================" << endl; 96 cout << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl; 97 cout << Verbose(0) << "INPUT: "; 93 98 cin >> choice; 94 99 95 100 switch (choice) { 96 101 default: 97 eLog() << Verbose(2) << "Not a valid choice." << endl;102 DoeLog(2) && (eLog()<< Verbose(2) << "Not a valid choice." << endl); 98 103 break; 99 104 case 'a': // absolute coordinates of atom 100 Log()<< Verbose(0) << "Enter absolute coordinates." << endl;105 cout << Verbose(0) << "Enter absolute coordinates." << endl; 101 106 first = new atom; 102 first->x.AskPosition( mol->cell_size, false);107 first->x.AskPosition(World::get()->cell_size, false); 103 108 first->type = periode->AskElement(); // give type 104 109 mol->AddAtom(first); // add to molecule … … 109 114 valid = true; 110 115 do { 111 if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;112 Log()<< Verbose(0) << "Enter reference coordinates." << endl;113 x.AskPosition( mol->cell_size, true);114 Log()<< Verbose(0) << "Enter relative coordinates." << endl;115 first->x.AskPosition( mol->cell_size, false);116 if (!valid) DoeLog(2) && (eLog()<< Verbose(2) << "Resulting position out of cell." << endl); 117 cout << Verbose(0) << "Enter reference coordinates." << endl; 118 x.AskPosition(World::get()->cell_size, true); 119 cout << Verbose(0) << "Enter relative coordinates." << endl; 120 first->x.AskPosition(World::get()->cell_size, false); 116 121 first->x.AddVector((const Vector *)&x); 117 Log()<< Verbose(0) << "\n";122 cout << Verbose(0) << "\n"; 118 123 } while (!(valid = mol->CheckBounds((const Vector *)&first->x))); 119 124 first->type = periode->AskElement(); // give type … … 125 130 valid = true; 126 131 do { 127 if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;132 if (!valid) DoeLog(2) && (eLog()<< Verbose(2) << "Resulting position out of cell." << endl); 128 133 second = mol->AskAtom("Enter atom number: "); 129 Log() << Verbose(0) << "Enter relative coordinates." << endl;130 first->x.AskPosition( mol->cell_size, false);134 DoLog(0) && (Log() << Verbose(0) << "Enter relative coordinates." << endl); 135 first->x.AskPosition(World::get()->cell_size, false); 131 136 for (int i=NDIM;i--;) { 132 137 first->x.x[i] += second->x.x[i]; … … 142 147 do { 143 148 if (!valid) { 144 eLog() << Verbose(2) << "Resulting coordinates out of cell - " << first->x << endl;149 DoeLog(2) && (eLog()<< Verbose(2) << "Resulting coordinates out of cell - " << first->x << endl); 145 150 } 146 Log()<< Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;151 cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl; 147 152 second = mol->AskAtom("Enter central atom: "); 148 153 third = mol->AskAtom("Enter second atom (specifying the axis for first angle): "); … … 155 160 c *= M_PI/180.; 156 161 bound(&c, -M_PI, M_PI); 157 Log()<< Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;162 cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl; 158 163 /* 159 164 second->Output(1,1,(ofstream *)&cout); … … 167 172 168 173 if (!z.SolveSystem(&x,&y,&n, b, c, a)) { 169 Log() << Verbose(0) << "Failure solving self-dependent linear system!" << endl;174 coutg() << Verbose(0) << "Failure solving self-dependent linear system!" << endl; 170 175 continue; 171 176 } 172 Log() << Verbose(0) << "resulting relative coordinates: ";177 DoLog(0) && (Log() << Verbose(0) << "resulting relative coordinates: "); 173 178 z.Output(); 174 Log() << Verbose(0) << endl;179 DoLog(0) && (Log() << Verbose(0) << endl); 175 180 */ 176 181 // calc axis vector … … 180 185 Log() << Verbose(0) << "x: ", 181 186 x.Output(); 182 Log() << Verbose(0) << endl;187 DoLog(0) && (Log() << Verbose(0) << endl); 183 188 z.MakeNormalVector(&second->x,&third->x,&fourth->x); 184 189 Log() << Verbose(0) << "z: ", 185 190 z.Output(); 186 Log() << Verbose(0) << endl;191 DoLog(0) && (Log() << Verbose(0) << endl); 187 192 y.MakeNormalVector(&x,&z); 188 193 Log() << Verbose(0) << "y: ", 189 194 y.Output(); 190 Log() << Verbose(0) << endl;195 DoLog(0) && (Log() << Verbose(0) << endl); 191 196 192 197 // rotate vector around first angle … … 195 200 Log() << Verbose(0) << "Rotated vector: ", 196 201 first->x.Output(); 197 Log() << Verbose(0) << endl;202 DoLog(0) && (Log() << Verbose(0) << endl); 198 203 // remove the projection onto the rotation plane of the second angle 199 204 n.CopyVector(&y); … … 201 206 Log() << Verbose(0) << "N1: ", 202 207 n.Output(); 203 Log() << Verbose(0) << endl;208 DoLog(0) && (Log() << Verbose(0) << endl); 204 209 first->x.SubtractVector(&n); 205 210 Log() << Verbose(0) << "Subtracted vector: ", 206 211 first->x.Output(); 207 Log() << Verbose(0) << endl;212 DoLog(0) && (Log() << Verbose(0) << endl); 208 213 n.CopyVector(&z); 209 214 n.Scale(first->x.ScalarProduct(&z)); 210 215 Log() << Verbose(0) << "N2: ", 211 216 n.Output(); 212 Log() << Verbose(0) << endl;217 DoLog(0) && (Log() << Verbose(0) << endl); 213 218 first->x.SubtractVector(&n); 214 219 Log() << Verbose(0) << "2nd subtracted vector: ", 215 220 first->x.Output(); 216 Log() << Verbose(0) << endl;221 DoLog(0) && (Log() << Verbose(0) << endl); 217 222 218 223 // rotate another vector around second angle … … 221 226 Log() << Verbose(0) << "2nd Rotated vector: ", 222 227 n.Output(); 223 Log() << Verbose(0) << endl;228 DoLog(0) && (Log() << Verbose(0) << endl); 224 229 225 230 // add the two linear independent vectors … … 229 234 first->x.AddVector(&second->x); 230 235 231 Log() << Verbose(0) << "resulting coordinates: ";236 DoLog(0) && (Log() << Verbose(0) << "resulting coordinates: "); 232 237 first->x.Output(); 233 Log() << Verbose(0) << endl;238 DoLog(0) && (Log() << Verbose(0) << endl); 234 239 } while (!(valid = mol->CheckBounds((const Vector *)&first->x))); 235 240 first->type = periode->AskElement(); // give type … … 244 249 atoms[i] = NULL; 245 250 int i=0, j=0; 246 Log()<< Verbose(0) << "Now we need at least three molecules.\n";251 cout << Verbose(0) << "Now we need at least three molecules.\n"; 247 252 do { 248 Log()<< Verbose(0) << "Enter " << i+1 << "th atom: ";253 cout << Verbose(0) << "Enter " << i+1 << "th atom: "; 249 254 cin >> j; 250 255 if (j != -1) { … … 261 266 } else { 262 267 delete first; 263 Log()<< Verbose(0) << "Please enter at least two vectors!\n";268 cout << Verbose(0) << "Please enter at least two vectors!\n"; 264 269 } 265 270 break; … … 275 280 char choice; // menu choice char 276 281 277 Log()<< Verbose(0) << "===========CENTER ATOMS=========================" << endl;278 Log()<< Verbose(0) << " a - on origin" << endl;279 Log()<< Verbose(0) << " b - on center of gravity" << endl;280 Log()<< Verbose(0) << " c - within box with additional boundary" << endl;281 Log()<< Verbose(0) << " d - within given simulation box" << endl;282 Log()<< Verbose(0) << "all else - go back" << endl;283 Log()<< Verbose(0) << "===============================================" << endl;284 Log()<< Verbose(0) << "INPUT: ";282 cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl; 283 cout << Verbose(0) << " a - on origin" << endl; 284 cout << Verbose(0) << " b - on center of gravity" << endl; 285 cout << Verbose(0) << " c - within box with additional boundary" << endl; 286 cout << Verbose(0) << " d - within given simulation box" << endl; 287 cout << Verbose(0) << "all else - go back" << endl; 288 cout << Verbose(0) << "===============================================" << endl; 289 cout << Verbose(0) << "INPUT: "; 285 290 cin >> choice; 286 291 287 292 switch (choice) { 288 293 default: 289 Log()<< Verbose(0) << "Not a valid choice." << endl;294 cout << Verbose(0) << "Not a valid choice." << endl; 290 295 break; 291 296 case 'a': 292 Log()<< Verbose(0) << "Centering atoms in config file on origin." << endl;297 cout << Verbose(0) << "Centering atoms in config file on origin." << endl; 293 298 mol->CenterOrigin(); 294 299 break; 295 300 case 'b': 296 Log()<< Verbose(0) << "Centering atoms in config file on center of gravity." << endl;301 cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl; 297 302 mol->CenterPeriodic(); 298 303 break; 299 304 case 'c': 300 Log()<< Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;305 cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl; 301 306 for (int i=0;i<NDIM;i++) { 302 Log()<< Verbose(0) << "Enter axis " << i << " boundary: ";307 cout << Verbose(0) << "Enter axis " << i << " boundary: "; 303 308 cin >> y.x[i]; 304 309 } … … 311 316 break; 312 317 case 'd': 313 Log()<< Verbose(1) << "Centering atoms in config file within given simulation box." << endl;318 cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl; 314 319 for (int i=0;i<NDIM;i++) { 315 Log()<< Verbose(0) << "Enter axis " << i << " boundary: ";320 cout << Verbose(0) << "Enter axis " << i << " boundary: "; 316 321 cin >> x.x[i]; 317 322 } … … 334 339 char choice; // menu choice char 335 340 336 Log()<< Verbose(0) << "===========ALIGN ATOMS=========================" << endl;337 Log()<< Verbose(0) << " a - state three atoms defining align plane" << endl;338 Log()<< Verbose(0) << " b - state alignment vector" << endl;339 Log()<< Verbose(0) << " c - state two atoms in alignment direction" << endl;340 Log()<< Verbose(0) << " d - align automatically by least square fit" << endl;341 Log()<< Verbose(0) << "all else - go back" << endl;342 Log()<< Verbose(0) << "===============================================" << endl;343 Log()<< Verbose(0) << "INPUT: ";341 cout << Verbose(0) << "===========ALIGN ATOMS=========================" << endl; 342 cout << Verbose(0) << " a - state three atoms defining align plane" << endl; 343 cout << Verbose(0) << " b - state alignment vector" << endl; 344 cout << Verbose(0) << " c - state two atoms in alignment direction" << endl; 345 cout << Verbose(0) << " d - align automatically by least square fit" << endl; 346 cout << Verbose(0) << "all else - go back" << endl; 347 cout << Verbose(0) << "===============================================" << endl; 348 cout << Verbose(0) << "INPUT: "; 344 349 cin >> choice; 345 350 … … 354 359 break; 355 360 case 'b': // normal vector of mirror plane 356 Log()<< Verbose(0) << "Enter normal vector of mirror plane." << endl;357 n.AskPosition( mol->cell_size,0);361 cout << Verbose(0) << "Enter normal vector of mirror plane." << endl; 362 n.AskPosition(World::get()->cell_size,0); 358 363 n.Normalize(); 359 364 break; … … 374 379 fscanf(stdin, "%3s", shorthand); 375 380 } while ((param.type = periode->FindElement(shorthand)) == NULL); 376 Log()<< Verbose(0) << "Element is " << param.type->name << endl;381 cout << Verbose(0) << "Element is " << param.type->name << endl; 377 382 mol->GetAlignvector(¶m); 378 383 for (int i=NDIM;i--;) { … … 381 386 } 382 387 gsl_vector_free(param.x); 383 Log()<< Verbose(0) << "Offset vector: ";388 cout << Verbose(0) << "Offset vector: "; 384 389 x.Output(); 385 Log() << Verbose(0) << endl;390 DoLog(0) && (Log() << Verbose(0) << endl); 386 391 n.Normalize(); 387 392 break; 388 393 }; 389 Log() << Verbose(0) << "Alignment vector: ";394 DoLog(0) && (Log() << Verbose(0) << "Alignment vector: "); 390 395 n.Output(); 391 Log() << Verbose(0) << endl;396 DoLog(0) && (Log() << Verbose(0) << endl); 392 397 mol->Align(&n); 393 398 }; … … 402 407 char choice; // menu choice char 403 408 404 Log() << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;405 Log() << Verbose(0) << " a - state three atoms defining mirror plane" << endl;406 Log() << Verbose(0) << " b - state normal vector of mirror plane" << endl;407 Log() << Verbose(0) << " c - state two atoms in normal direction" << endl;408 Log() << Verbose(0) << "all else - go back" << endl;409 Log() << Verbose(0) << "===============================================" << endl;410 Log() << Verbose(0) << "INPUT: ";409 DoLog(0) && (Log() << Verbose(0) << "===========MIRROR ATOMS=========================" << endl); 410 DoLog(0) && (Log() << Verbose(0) << " a - state three atoms defining mirror plane" << endl); 411 DoLog(0) && (Log() << Verbose(0) << " b - state normal vector of mirror plane" << endl); 412 DoLog(0) && (Log() << Verbose(0) << " c - state two atoms in normal direction" << endl); 413 DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl); 414 DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl); 415 DoLog(0) && (Log() << Verbose(0) << "INPUT: "); 411 416 cin >> choice; 412 417 … … 421 426 break; 422 427 case 'b': // normal vector of mirror plane 423 Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;424 n.AskPosition( mol->cell_size,0);428 DoLog(0) && (Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl); 429 n.AskPosition(World::get()->cell_size,0); 425 430 n.Normalize(); 426 431 break; … … 434 439 break; 435 440 }; 436 Log() << Verbose(0) << "Normal vector: ";441 DoLog(0) && (Log() << Verbose(0) << "Normal vector: "); 437 442 n.Output(); 438 Log() << Verbose(0) << endl;443 DoLog(0) && (Log() << Verbose(0) << endl); 439 444 mol->Mirror((const Vector *)&n); 440 445 }; … … 450 455 char choice; // menu choice char 451 456 452 Log() << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;453 Log() << Verbose(0) << " a - state atom for removal by number" << endl;454 Log() << Verbose(0) << " b - keep only in radius around atom" << endl;455 Log() << Verbose(0) << " c - remove this with one axis greater value" << endl;456 Log() << Verbose(0) << "all else - go back" << endl;457 Log() << Verbose(0) << "===============================================" << endl;458 Log() << Verbose(0) << "INPUT: ";457 DoLog(0) && (Log() << Verbose(0) << "===========REMOVE ATOMS=========================" << endl); 458 DoLog(0) && (Log() << Verbose(0) << " a - state atom for removal by number" << endl); 459 DoLog(0) && (Log() << Verbose(0) << " b - keep only in radius around atom" << endl); 460 DoLog(0) && (Log() << Verbose(0) << " c - remove this with one axis greater value" << endl); 461 DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl); 462 DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl); 463 DoLog(0) && (Log() << Verbose(0) << "INPUT: "); 459 464 cin >> choice; 460 465 … … 463 468 case 'a': 464 469 if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: "))) 465 Log() << Verbose(1) << "Atom removed." << endl;470 DoLog(1) && (Log() << Verbose(1) << "Atom removed." << endl); 466 471 else 467 Log() << Verbose(1) << "Atom not found." << endl;472 DoLog(1) && (Log() << Verbose(1) << "Atom not found." << endl); 468 473 break; 469 474 case 'b': 470 475 second = mol->AskAtom("Enter number of atom as reference point: "); 471 Log() << Verbose(0) << "Enter radius: ";476 DoLog(0) && (Log() << Verbose(0) << "Enter radius: "); 472 477 cin >> tmp1; 473 478 first = mol->start; … … 481 486 break; 482 487 case 'c': 483 Log() << Verbose(0) << "Which axis is it: ";488 DoLog(0) && (Log() << Verbose(0) << "Which axis is it: "); 484 489 cin >> axis; 485 Log() << Verbose(0) << "Lower boundary: ";490 DoLog(0) && (Log() << Verbose(0) << "Lower boundary: "); 486 491 cin >> tmp1; 487 Log() << Verbose(0) << "Upper boundary: ";492 DoLog(0) && (Log() << Verbose(0) << "Upper boundary: "); 488 493 cin >> tmp2; 489 494 first = mol->start; … … 515 520 char choice; // menu choice char 516 521 517 Log() << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;518 Log() << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;519 Log() << Verbose(0) << " b - calculate bond length between two atoms" << endl;520 Log() << Verbose(0) << " c - calculate bond angle" << endl;521 Log() << Verbose(0) << " d - calculate principal axis of the system" << endl;522 Log() << Verbose(0) << " e - calculate volume of the convex envelope" << endl;523 Log() << Verbose(0) << " f - calculate temperature from current velocity" << endl;524 Log() << Verbose(0) << " g - output all temperatures per step from velocities" << endl;525 Log() << Verbose(0) << "all else - go back" << endl;526 Log() << Verbose(0) << "===============================================" << endl;527 Log() << Verbose(0) << "INPUT: ";522 DoLog(0) && (Log() << Verbose(0) << "===========MEASURE ATOMS=========================" << endl); 523 DoLog(0) && (Log() << Verbose(0) << " a - calculate bond length between one atom and all others" << endl); 524 DoLog(0) && (Log() << Verbose(0) << " b - calculate bond length between two atoms" << endl); 525 DoLog(0) && (Log() << Verbose(0) << " c - calculate bond angle" << endl); 526 DoLog(0) && (Log() << Verbose(0) << " d - calculate principal axis of the system" << endl); 527 DoLog(0) && (Log() << Verbose(0) << " e - calculate volume of the convex envelope" << endl); 528 DoLog(0) && (Log() << Verbose(0) << " f - calculate temperature from current velocity" << endl); 529 DoLog(0) && (Log() << Verbose(0) << " g - output all temperatures per step from velocities" << endl); 530 DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl); 531 DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl); 532 DoLog(0) && (Log() << Verbose(0) << "INPUT: "); 528 533 cin >> choice; 529 534 530 535 switch(choice) { 531 536 default: 532 Log() << Verbose(1) << "Not a valid choice." << endl;537 DoLog(1) && (Log() << Verbose(1) << "Not a valid choice." << endl); 533 538 break; 534 539 case 'a': … … 562 567 x.SubtractVector((const Vector *)&second->x); 563 568 tmp1 = x.Norm(); 564 Log() << Verbose(1) << "Distance vector is ";569 DoLog(1) && (Log() << Verbose(1) << "Distance vector is "); 565 570 x.Output(); 566 Log() << Verbose(0) << "." << endl << "Norm of distance is " << tmp1 << "." << endl;571 DoLog(0) && (Log() << Verbose(0) << "." << endl << "Norm of distance is " << tmp1 << "." << endl); 567 572 break; 568 573 569 574 case 'c': 570 Log() << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl;575 DoLog(0) && (Log() << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl); 571 576 first = mol->AskAtom("Enter first atom: "); 572 577 second = mol->AskAtom("Enter central atom: "); … … 577 582 y.CopyVector((const Vector *)&third->x); 578 583 y.SubtractVector((const Vector *)&second->x); 579 Log() << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";580 Log() << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;584 DoLog(0) && (Log() << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": "); 585 DoLog(0) && (Log() << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl); 581 586 break; 582 587 case 'd': 583 Log() << Verbose(0) << "Evaluating prinicipal axis." << endl;584 Log() << Verbose(0) << "Shall we rotate? [0/1]: ";588 DoLog(0) && (Log() << Verbose(0) << "Evaluating prinicipal axis." << endl); 589 DoLog(0) && (Log() << Verbose(0) << "Shall we rotate? [0/1]: "); 585 590 cin >> Z; 586 591 if ((Z >=0) && (Z <=1)) … … 591 596 case 'e': 592 597 { 593 Log() << Verbose(0) << "Evaluating volume of the convex envelope.";598 DoLog(0) && (Log() << Verbose(0) << "Evaluating volume of the convex envelope."); 594 599 class Tesselation *TesselStruct = NULL; 595 600 const LinkedCell *LCList = NULL; … … 597 602 FindConvexBorder(mol, TesselStruct, LCList, NULL); 598 603 double clustervolume = VolumeOfConvexEnvelope(TesselStruct, configuration); 599 Log() << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl;\604 DoLog(0) && (Log() << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl);\ 600 605 delete(LCList); 601 606 delete(TesselStruct); … … 608 613 { 609 614 char filename[255]; 610 Log() << Verbose(0) << "Please enter filename: " << endl;615 DoLog(0) && (Log() << Verbose(0) << "Please enter filename: " << endl); 611 616 cin >> filename; 612 Log() << Verbose(1) << "Storing temperatures in " << filename << "." << endl;617 DoLog(1) && (Log() << Verbose(1) << "Storing temperatures in " << filename << "." << endl); 613 618 ofstream *output = new ofstream(filename, ios::trunc); 614 619 if (!mol->OutputTemperatureFromTrajectories(output, 0, mol->MDSteps)) 615 Log() << Verbose(2) << "File could not be written." << endl;620 DoLog(2) && (Log() << Verbose(2) << "File could not be written." << endl); 616 621 else 617 Log() << Verbose(2) << "File stored." << endl;622 DoLog(2) && (Log() << Verbose(2) << "File stored." << endl); 618 623 output->close(); 619 624 delete(output); … … 632 637 clock_t start, end; 633 638 634 Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;635 Log() << Verbose(0) << "What's the desired bond order: ";639 DoLog(0) && (Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl); 640 DoLog(0) && (Log() << Verbose(0) << "What's the desired bond order: "); 636 641 cin >> Order1; 637 642 if (mol->first->next != mol->last) { // there are bonds … … 639 644 mol->FragmentMolecule(Order1, configuration); 640 645 end = clock(); 641 Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;646 DoLog(0) && (Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl); 642 647 } else 643 Log() << Verbose(0) << "Connection matrix has not yet been generated!" << endl;648 DoLog(0) && (Log() << Verbose(0) << "Connection matrix has not yet been generated!" << endl); 644 649 }; 645 650 … … 652 657 static void ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration) 653 658 { 654 atom *first, *second ;659 atom *first, *second, *third; 655 660 molecule *mol = NULL; 656 661 Vector x,y,z,n; // coordinates for absolute point in cell volume … … 660 665 bool valid; 661 666 662 Log() << Verbose(0) << "=========MANIPULATE ATOMS======================" << endl; 663 Log() << Verbose(0) << "a - add an atom" << endl; 664 Log() << Verbose(0) << "r - remove an atom" << endl; 665 Log() << Verbose(0) << "b - scale a bond between atoms" << endl; 666 Log() << Verbose(0) << "u - change an atoms element" << endl; 667 Log() << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl; 668 Log() << Verbose(0) << "all else - go back" << endl; 669 Log() << Verbose(0) << "===============================================" << endl; 667 DoLog(0) && (Log() << Verbose(0) << "=========MANIPULATE ATOMS======================" << endl); 668 DoLog(0) && (Log() << Verbose(0) << "a - add an atom" << endl); 669 DoLog(0) && (Log() << Verbose(0) << "r - remove an atom" << endl); 670 DoLog(0) && (Log() << Verbose(0) << "b - scale a bond between atoms" << endl); 671 DoLog(0) && (Log() << Verbose(0) << "t - turn an atom round another bond" << endl); 672 DoLog(0) && (Log() << Verbose(0) << "u - change an atoms element" << endl); 673 DoLog(0) && (Log() << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl); 674 DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl); 675 DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl); 670 676 if (molecules->NumberOfActiveMolecules() > 1) 671 eLog() << Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl;672 Log() << Verbose(0) << "INPUT: ";677 DoeLog(2) && (eLog()<< Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl); 678 DoLog(0) && (Log() << Verbose(0) << "INPUT: "); 673 679 cin >> choice; 674 680 675 681 switch (choice) { 676 682 default: 677 Log() << Verbose(0) << "Not a valid choice." << endl;683 DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl); 678 684 break; 679 685 … … 682 688 if ((*ListRunner)->ActiveFlag) { 683 689 mol = *ListRunner; 684 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;690 DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl); 685 691 AddAtoms(periode, mol); 686 692 } … … 691 697 if ((*ListRunner)->ActiveFlag) { 692 698 mol = *ListRunner; 693 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;694 Log() << Verbose(0) << "Scaling bond length between two atoms." << endl;699 DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl); 700 DoLog(0) && (Log() << Verbose(0) << "Scaling bond length between two atoms." << endl); 695 701 first = mol->AskAtom("Enter first (fixed) atom: "); 696 702 second = mol->AskAtom("Enter second (shifting) atom: "); … … 699 705 minBond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]); 700 706 minBond = sqrt(minBond); 701 Log() << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << minBond << " a.u." << endl;702 Log() << Verbose(0) << "Enter new bond length [a.u.]: ";707 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); 708 DoLog(0) && (Log() << Verbose(0) << "Enter new bond length [a.u.]: "); 703 709 cin >> bond; 704 710 for (int i=NDIM;i--;) { … … 714 720 if ((*ListRunner)->ActiveFlag) { 715 721 mol = *ListRunner; 716 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;717 Log() << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;718 Log() << Verbose(0) << "Enter three factors: ";722 DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl); 723 DoLog(0) && (Log() << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl); 724 DoLog(0) && (Log() << Verbose(0) << "Enter three factors: "); 719 725 factor = new double[NDIM]; 720 726 cin >> factor[0]; … … 731 737 if ((*ListRunner)->ActiveFlag) { 732 738 mol = *ListRunner; 733 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;739 DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl); 734 740 MeasureAtoms(periode, mol, configuration); 735 741 } … … 740 746 if ((*ListRunner)->ActiveFlag) { 741 747 mol = *ListRunner; 742 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;748 DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl); 743 749 RemoveAtoms(mol); 744 750 } 751 break; 752 753 case 't': // turn/rotate atom 754 for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) 755 if ((*ListRunner)->ActiveFlag) { 756 mol = *ListRunner; 757 DoLog(0) && (Log() << Verbose(0) << "Turning atom around another bond - first is atom to turn, second (central) and third specify bond" << endl); 758 first = mol->AskAtom("Enter turning atom: "); 759 second = mol->AskAtom("Enter central atom: "); 760 third = mol->AskAtom("Enter bond atom: "); 761 cout << Verbose(0) << "Enter new angle in degrees: "; 762 double tmp = 0.; 763 cin >> tmp; 764 // calculate old angle 765 x.CopyVector((const Vector *)&first->x); 766 x.SubtractVector((const Vector *)&second->x); 767 y.CopyVector((const Vector *)&third->x); 768 y.SubtractVector((const Vector *)&second->x); 769 double alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.); 770 cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": "; 771 cout << Verbose(0) << alpha << " degrees" << endl; 772 // rotate 773 z.MakeNormalVector(&x,&y); 774 x.RotateVector(&z,(alpha-tmp)*M_PI/180.); 775 x.AddVector(&second->x); 776 first->x.CopyVector(&x); 777 // check new angle 778 x.CopyVector((const Vector *)&first->x); 779 x.SubtractVector((const Vector *)&second->x); 780 alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.); 781 cout << Verbose(0) << "new Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": "; 782 cout << Verbose(0) << alpha << " degrees" << endl; 783 } 745 784 break; 746 785 … … 750 789 int Z; 751 790 mol = *ListRunner; 752 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;791 DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl); 753 792 first = NULL; 754 793 do { 755 Log() << Verbose(0) << "Change the element of which atom: ";794 DoLog(0) && (Log() << Verbose(0) << "Change the element of which atom: "); 756 795 cin >> Z; 757 796 } while ((first = mol->FindAtom(Z)) == NULL); 758 Log() << Verbose(0) << "New element by atomic number Z: ";797 DoLog(0) && (Log() << Verbose(0) << "New element by atomic number Z: "); 759 798 cin >> Z; 760 799 first->type = periode->FindElement(Z); 761 Log() << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;800 DoLog(0) && (Log() << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl); 762 801 } 763 802 break; … … 780 819 MoleculeLeafClass *Subgraphs = NULL; 781 820 782 Log() << Verbose(0) << "=========MANIPULATE GLOBALLY===================" << endl;783 Log() << Verbose(0) << "c - scale by unit transformation" << endl;784 Log() << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;785 Log() << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;786 Log() << Verbose(0) << "g - center atoms in box" << endl;787 Log() << Verbose(0) << "i - realign molecule" << endl;788 Log() << Verbose(0) << "m - mirror all molecules" << endl;789 Log() << Verbose(0) << "o - create connection matrix" << endl;790 Log() << Verbose(0) << "t - translate molecule by vector" << endl;791 Log() << Verbose(0) << "all else - go back" << endl;792 Log() << Verbose(0) << "===============================================" << endl;821 DoLog(0) && (Log() << Verbose(0) << "=========MANIPULATE GLOBALLY===================" << endl); 822 DoLog(0) && (Log() << Verbose(0) << "c - scale by unit transformation" << endl); 823 DoLog(0) && (Log() << Verbose(0) << "d - duplicate molecule/periodic cell" << endl); 824 DoLog(0) && (Log() << Verbose(0) << "f - fragment molecule many-body bond order style" << endl); 825 DoLog(0) && (Log() << Verbose(0) << "g - center atoms in box" << endl); 826 DoLog(0) && (Log() << Verbose(0) << "i - realign molecule" << endl); 827 DoLog(0) && (Log() << Verbose(0) << "m - mirror all molecules" << endl); 828 DoLog(0) && (Log() << Verbose(0) << "o - create connection matrix" << endl); 829 DoLog(0) && (Log() << Verbose(0) << "t - translate molecule by vector" << endl); 830 DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl); 831 DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl); 793 832 if (molecules->NumberOfActiveMolecules() > 1) 794 eLog() << Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl;795 Log() << Verbose(0) << "INPUT: ";833 DoeLog(2) && (eLog()<< Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl); 834 DoLog(0) && (Log() << Verbose(0) << "INPUT: "); 796 835 cin >> choice; 797 836 798 837 switch (choice) { 799 838 default: 800 Log() << Verbose(0) << "Not a valid choice." << endl;839 DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl); 801 840 break; 802 841 … … 805 844 if ((*ListRunner)->ActiveFlag) { 806 845 mol = *ListRunner; 807 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;808 Log() << Verbose(0) << "State the axis [(+-)123]: ";846 DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl); 847 DoLog(0) && (Log() << Verbose(0) << "State the axis [(+-)123]: "); 809 848 cin >> axis; 810 Log() << Verbose(0) << "State the factor: ";849 DoLog(0) && (Log() << Verbose(0) << "State the factor: "); 811 850 cin >> faktor; 812 851 … … 825 864 } 826 865 if (count != j) 827 eLog() << Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;866 DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl); 828 867 x.Zero(); 829 868 y.Zero(); 830 y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude869 y.x[abs(axis)-1] = World::get()->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude 831 870 for (int i=1;i<faktor;i++) { // then add this list with respective translation factor times 832 871 x.AddVector(&y); // per factor one cell width further … … 851 890 mol->Translate(&x); 852 891 } 853 mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;892 World::get()->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor; 854 893 } 855 894 } … … 864 903 if ((*ListRunner)->ActiveFlag) { 865 904 mol = *ListRunner; 866 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;905 DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl); 867 906 CenterAtoms(mol); 868 907 } … … 873 912 if ((*ListRunner)->ActiveFlag) { 874 913 mol = *ListRunner; 875 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;914 DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl); 876 915 AlignAtoms(periode, mol); 877 916 } … … 882 921 if ((*ListRunner)->ActiveFlag) { 883 922 mol = *ListRunner; 884 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;923 DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl); 885 924 MirrorAtoms(mol); 886 925 } … … 893 932 double bonddistance; 894 933 clock_t start,end; 895 Log() << Verbose(0) << "What's the maximum bond distance: ";934 DoLog(0) && (Log() << Verbose(0) << "What's the maximum bond distance: "); 896 935 cin >> bonddistance; 897 936 start = clock(); 898 937 mol->CreateAdjacencyList(bonddistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL); 899 938 end = clock(); 900 Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;939 DoLog(0) && (Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl); 901 940 } 902 941 break; … … 906 945 if ((*ListRunner)->ActiveFlag) { 907 946 mol = *ListRunner; 908 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;909 Log() << Verbose(0) << "Enter translation vector." << endl;910 x.AskPosition( mol->cell_size,0);947 DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl); 948 DoLog(0) && (Log() << Verbose(0) << "Enter translation vector." << endl); 949 x.AskPosition(World::get()->cell_size,0); 911 950 mol->Center.AddVector((const Vector *)&x); 912 951 } … … 935 974 molecule *mol = NULL; 936 975 937 Log() << Verbose(0) << "==========EDIT MOLECULES=====================" << endl;938 Log() << Verbose(0) << "c - create new molecule" << endl;939 Log() << Verbose(0) << "l - load molecule from xyz file" << endl;940 Log() << Verbose(0) << "n - change molecule's name" << endl;941 Log() << Verbose(0) << "N - give molecules filename" << endl;942 Log() << Verbose(0) << "p - parse atoms in xyz file into molecule" << endl;943 Log() << Verbose(0) << "r - remove a molecule" << endl;944 Log() << Verbose(0) << "all else - go back" << endl;945 Log() << Verbose(0) << "===============================================" << endl;946 Log() << Verbose(0) << "INPUT: ";976 DoLog(0) && (Log() << Verbose(0) << "==========EDIT MOLECULES=====================" << endl); 977 DoLog(0) && (Log() << Verbose(0) << "c - create new molecule" << endl); 978 DoLog(0) && (Log() << Verbose(0) << "l - load molecule from xyz file" << endl); 979 DoLog(0) && (Log() << Verbose(0) << "n - change molecule's name" << endl); 980 DoLog(0) && (Log() << Verbose(0) << "N - give molecules filename" << endl); 981 DoLog(0) && (Log() << Verbose(0) << "p - parse atoms in xyz file into molecule" << endl); 982 DoLog(0) && (Log() << Verbose(0) << "r - remove a molecule" << endl); 983 DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl); 984 DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl); 985 DoLog(0) && (Log() << Verbose(0) << "INPUT: "); 947 986 cin >> choice; 948 987 949 988 switch (choice) { 950 989 default: 951 Log() << Verbose(0) << "Not a valid choice." << endl;990 DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl); 952 991 break; 953 992 case 'c': … … 959 998 { 960 999 char filename[MAXSTRINGSIZE]; 961 Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;1000 DoLog(0) && (Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl); 962 1001 mol = new molecule(periode); 963 1002 do { 964 Log() << Verbose(0) << "Enter file name: ";1003 DoLog(0) && (Log() << Verbose(0) << "Enter file name: "); 965 1004 cin >> filename; 966 1005 } while (!mol->AddXYZFile(filename)); … … 968 1007 // center at set box dimensions 969 1008 mol->CenterEdge(¢er); 970 mol->cell_size[0] = center.x[0]; 971 mol->cell_size[1] = 0; 972 mol->cell_size[2] = center.x[1]; 973 mol->cell_size[3] = 0; 974 mol->cell_size[4] = 0; 975 mol->cell_size[5] = center.x[2]; 1009 double * const cell_size = World::get()->cell_size; 1010 cell_size[0] = center.x[0]; 1011 cell_size[1] = 0; 1012 cell_size[2] = center.x[1]; 1013 cell_size[3] = 0; 1014 cell_size[4] = 0; 1015 cell_size[5] = center.x[2]; 976 1016 molecules->insert(mol); 977 1017 } … … 982 1022 char filename[MAXSTRINGSIZE]; 983 1023 do { 984 Log() << Verbose(0) << "Enter index of molecule: ";1024 DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: "); 985 1025 cin >> nr; 986 1026 mol = molecules->ReturnIndex(nr); 987 1027 } while (mol == NULL); 988 Log() << Verbose(0) << "Enter name: ";1028 DoLog(0) && (Log() << Verbose(0) << "Enter name: "); 989 1029 cin >> filename; 990 1030 strcpy(mol->name, filename); … … 996 1036 char filename[MAXSTRINGSIZE]; 997 1037 do { 998 Log() << Verbose(0) << "Enter index of molecule: ";1038 DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: "); 999 1039 cin >> nr; 1000 1040 mol = molecules->ReturnIndex(nr); 1001 1041 } while (mol == NULL); 1002 Log() << Verbose(0) << "Enter name: ";1042 DoLog(0) && (Log() << Verbose(0) << "Enter name: "); 1003 1043 cin >> filename; 1004 1044 mol->SetNameFromFilename(filename); … … 1011 1051 mol = NULL; 1012 1052 do { 1013 Log() << Verbose(0) << "Enter index of molecule: ";1053 DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: "); 1014 1054 cin >> nr; 1015 1055 mol = molecules->ReturnIndex(nr); 1016 1056 } while (mol == NULL); 1017 Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;1057 DoLog(0) && (Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl); 1018 1058 do { 1019 Log() << Verbose(0) << "Enter file name: ";1059 DoLog(0) && (Log() << Verbose(0) << "Enter file name: "); 1020 1060 cin >> filename; 1021 1061 } while (!mol->AddXYZFile(filename)); … … 1025 1065 1026 1066 case 'r': 1027 Log() << Verbose(0) << "Enter index of molecule: ";1067 DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: "); 1028 1068 cin >> nr; 1029 1069 count = 1; … … 1048 1088 char choice; // menu choice char 1049 1089 1050 Log() << Verbose(0) << "===========MERGE MOLECULES=====================" << endl; 1051 Log() << Verbose(0) << "a - simple add of one molecule to another" << endl; 1052 Log() << Verbose(0) << "e - embedding merge of two molecules" << endl; 1053 Log() << Verbose(0) << "m - multi-merge of all molecules" << endl; 1054 Log() << Verbose(0) << "s - scatter merge of two molecules" << endl; 1055 Log() << Verbose(0) << "t - simple merge of two molecules" << endl; 1056 Log() << Verbose(0) << "all else - go back" << endl; 1057 Log() << Verbose(0) << "===============================================" << endl; 1058 Log() << Verbose(0) << "INPUT: "; 1090 DoLog(0) && (Log() << Verbose(0) << "===========MERGE MOLECULES=====================" << endl); 1091 DoLog(0) && (Log() << Verbose(0) << "a - simple add of one molecule to another" << endl); 1092 DoLog(0) && (Log() << Verbose(0) << "b - count the number of bonds of two elements" << endl); 1093 DoLog(0) && (Log() << Verbose(0) << "B - count the number of bonds of three elements " << endl); 1094 DoLog(0) && (Log() << Verbose(0) << "e - embedding merge of two molecules" << endl); 1095 DoLog(0) && (Log() << Verbose(0) << "h - count the number of hydrogen bonds" << endl); 1096 DoLog(0) && (Log() << Verbose(0) << "b - count the number of hydrogen bonds" << endl); 1097 DoLog(0) && (Log() << Verbose(0) << "m - multi-merge of all molecules" << endl); 1098 DoLog(0) && (Log() << Verbose(0) << "s - scatter merge of two molecules" << endl); 1099 DoLog(0) && (Log() << Verbose(0) << "t - simple merge of two molecules" << endl); 1100 DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl); 1101 DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl); 1102 DoLog(0) && (Log() << Verbose(0) << "INPUT: "); 1059 1103 cin >> choice; 1060 1104 1061 1105 switch (choice) { 1062 1106 default: 1063 Log() << Verbose(0) << "Not a valid choice." << endl;1107 DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl); 1064 1108 break; 1065 1109 … … 1070 1114 { 1071 1115 do { 1072 Log() << Verbose(0) << "Enter index of destination molecule: ";1116 DoLog(0) && (Log() << Verbose(0) << "Enter index of destination molecule: "); 1073 1117 cin >> dest; 1074 1118 destmol = molecules->ReturnIndex(dest); 1075 1119 } while ((destmol == NULL) && (dest != -1)); 1076 1120 do { 1077 Log() << Verbose(0) << "Enter index of source molecule to add from: ";1121 DoLog(0) && (Log() << Verbose(0) << "Enter index of source molecule to add from: "); 1078 1122 cin >> src; 1079 1123 srcmol = molecules->ReturnIndex(src); … … 1085 1129 break; 1086 1130 1131 case 'b': 1132 { 1133 const int nr = 2; 1134 char *names[nr] = {"first", "second"}; 1135 int Z[nr]; 1136 element *elements[nr]; 1137 for (int i=0;i<nr;i++) { 1138 Z[i] = 0; 1139 do { 1140 cout << "Enter " << names[i] << " element: "; 1141 cin >> Z[i]; 1142 } while ((Z[i] <= 0) && (Z[i] > MAX_ELEMENTS)); 1143 elements[i] = periode->FindElement(Z[i]); 1144 } 1145 const int count = CountBondsOfTwo(molecules, elements[0], elements[1]); 1146 cout << endl << "There are " << count << " "; 1147 for (int i=0;i<nr;i++) { 1148 if (i==0) 1149 cout << elements[i]->symbol; 1150 else 1151 cout << "-" << elements[i]->symbol; 1152 } 1153 cout << " bonds." << endl; 1154 } 1155 break; 1156 1157 case 'B': 1158 { 1159 const int nr = 3; 1160 char *names[nr] = {"first", "second", "third"}; 1161 int Z[nr]; 1162 element *elements[nr]; 1163 for (int i=0;i<nr;i++) { 1164 Z[i] = 0; 1165 do { 1166 cout << "Enter " << names[i] << " element: "; 1167 cin >> Z[i]; 1168 } while ((Z[i] <= 0) && (Z[i] > MAX_ELEMENTS)); 1169 elements[i] = periode->FindElement(Z[i]); 1170 } 1171 const int count = CountBondsOfThree(molecules, elements[0], elements[1], elements[2]); 1172 cout << endl << "There are " << count << " "; 1173 for (int i=0;i<nr;i++) { 1174 if (i==0) 1175 cout << elements[i]->symbol; 1176 else 1177 cout << "-" << elements[i]->symbol; 1178 } 1179 cout << " bonds." << endl; 1180 } 1181 break; 1182 1087 1183 case 'e': 1088 1184 { … … 1090 1186 molecule *srcmol = NULL, *destmol = NULL; 1091 1187 do { 1092 Log() << Verbose(0) << "Enter index of matrix molecule (the variable one): ";1188 DoLog(0) && (Log() << Verbose(0) << "Enter index of matrix molecule (the variable one): "); 1093 1189 cin >> src; 1094 1190 srcmol = molecules->ReturnIndex(src); 1095 1191 } while ((srcmol == NULL) && (src != -1)); 1096 1192 do { 1097 Log() << Verbose(0) << "Enter index of molecule to merge into (the fixed one): ";1193 DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule to merge into (the fixed one): "); 1098 1194 cin >> dest; 1099 1195 destmol = molecules->ReturnIndex(dest); … … 1104 1200 break; 1105 1201 1202 case 'h': 1203 { 1204 int Z; 1205 cout << "Please enter interface element: "; 1206 cin >> Z; 1207 element * const InterfaceElement = periode->FindElement(Z); 1208 cout << endl << "There are " << CountHydrogenBridgeBonds(molecules, InterfaceElement) << " hydrogen bridges with connections to " << (InterfaceElement != 0 ? InterfaceElement->name : "None") << "." << endl; 1209 } 1210 break; 1211 1106 1212 case 'm': 1107 1213 { … … 1109 1215 molecule *mol = NULL; 1110 1216 do { 1111 Log() << Verbose(0) << "Enter index of molecule to merge into: ";1217 DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule to merge into: "); 1112 1218 cin >> nr; 1113 1219 mol = molecules->ReturnIndex(nr); … … 1126 1232 1127 1233 case 's': 1128 Log() << Verbose(0) << "Not implemented yet." << endl;1234 DoLog(0) && (Log() << Verbose(0) << "Not implemented yet." << endl); 1129 1235 break; 1130 1236 … … 1135 1241 { 1136 1242 do { 1137 Log() << Verbose(0) << "Enter index of destination molecule: ";1243 DoLog(0) && (Log() << Verbose(0) << "Enter index of destination molecule: "); 1138 1244 cin >> dest; 1139 1245 destmol = molecules->ReturnIndex(dest); 1140 1246 } while ((destmol == NULL) && (dest != -1)); 1141 1247 do { 1142 Log() << Verbose(0) << "Enter index of source molecule to merge into: ";1248 DoLog(0) && (Log() << Verbose(0) << "Enter index of source molecule to merge into: "); 1143 1249 cin >> src; 1144 1250 srcmol = molecules->ReturnIndex(src); … … 1169 1275 mol = (molecules->ListOfMolecules.front())->CopyMolecule(); 1170 1276 else { 1171 eLog() << Verbose(0) << "I don't have anything to test on ... ";1277 DoeLog(0) && (eLog()<< Verbose(0) << "I don't have anything to test on ... "); 1172 1278 performCriticalExit(); 1173 1279 return; … … 1176 1282 1177 1283 // generate some KeySets 1178 Log() << Verbose(0) << "Generating KeySets." << endl;1284 DoLog(0) && (Log() << Verbose(0) << "Generating KeySets." << endl); 1179 1285 KeySet TestSets[mol->AtomCount+1]; 1180 1286 i=1; … … 1186 1292 i++; 1187 1293 } 1188 Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl;1294 DoLog(0) && (Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl); 1189 1295 KeySetTestPair test; 1190 1296 test = TestSets[mol->AtomCount-1].insert(Walker->nr); 1191 1297 if (test.second) { 1192 Log() << Verbose(1) << "Insertion worked?!" << endl;1298 DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl); 1193 1299 } else { 1194 Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;1300 DoLog(1) && (Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl); 1195 1301 } 1196 1302 TestSets[mol->AtomCount].insert(mol->end->previous->nr); … … 1198 1304 1199 1305 // constructing Graph structure 1200 Log() << Verbose(0) << "Generating Subgraph class." << endl;1306 DoLog(0) && (Log() << Verbose(0) << "Generating Subgraph class." << endl); 1201 1307 Graph Subgraphs; 1202 1308 1203 1309 // insert KeySets into Subgraphs 1204 Log() << Verbose(0) << "Inserting KeySets into Subgraph class." << endl;1310 DoLog(0) && (Log() << Verbose(0) << "Inserting KeySets into Subgraph class." << endl); 1205 1311 for (int j=0;j<mol->AtomCount;j++) { 1206 1312 Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.))); 1207 1313 } 1208 Log() << Verbose(0) << "Testing insertion of already present item in Subgraph." << endl;1314 DoLog(0) && (Log() << Verbose(0) << "Testing insertion of already present item in Subgraph." << endl); 1209 1315 GraphTestPair test2; 1210 1316 test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.))); 1211 1317 if (test2.second) { 1212 Log() << Verbose(1) << "Insertion worked?!" << endl;1318 DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl); 1213 1319 } else { 1214 Log() << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;1320 DoLog(1) && (Log() << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl); 1215 1321 } 1216 1322 1217 1323 // show graphs 1218 Log() << Verbose(0) << "Showing Subgraph's contents, checking that it's sorted." << endl;1324 DoLog(0) && (Log() << Verbose(0) << "Showing Subgraph's contents, checking that it's sorted." << endl); 1219 1325 Graph::iterator A = Subgraphs.begin(); 1220 1326 while (A != Subgraphs.end()) { 1221 Log() << Verbose(0) << (*A).second.first << ": ";1327 DoLog(0) && (Log() << Verbose(0) << (*A).second.first << ": "); 1222 1328 KeySet::iterator key = (*A).first.begin(); 1223 1329 comp = -1; 1224 1330 while (key != (*A).first.end()) { 1225 1331 if ((*key) > comp) 1226 Log() << Verbose(0) << (*key) << " ";1332 DoLog(0) && (Log() << Verbose(0) << (*key) << " "); 1227 1333 else 1228 Log() << Verbose(0) << (*key) << "! ";1334 DoLog(0) && (Log() << Verbose(0) << (*key) << "! "); 1229 1335 comp = (*key); 1230 1336 key++; 1231 1337 } 1232 Log() << Verbose(0) << endl;1338 DoLog(0) && (Log() << Verbose(0) << endl); 1233 1339 A++; 1234 1340 } … … 1250 1356 1251 1357 if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) { 1252 eLog() << Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl;1358 DoeLog(2) && (eLog()<< Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl); 1253 1359 } 1254 1360 … … 1259 1365 if (output == NULL) 1260 1366 strcpy(filename,"main_pcp_linux"); 1261 Log() << Verbose(0) << "Saving as pdb input ";1367 DoLog(0) && (Log() << Verbose(0) << "Saving as pdb input "); 1262 1368 if (configuration->SavePDB(filename, molecules)) 1263 Log() << Verbose(0) << "done." << endl;1369 DoLog(0) && (Log() << Verbose(0) << "done." << endl); 1264 1370 else 1265 Log() << Verbose(0) << "failed." << endl;1371 DoLog(0) && (Log() << Verbose(0) << "failed." << endl); 1266 1372 1267 1373 // then save as tremolo data file … … 1270 1376 if (output == NULL) 1271 1377 strcpy(filename,"main_pcp_linux"); 1272 Log() << Verbose(0) << "Saving as tremolo data input ";1378 DoLog(0) && (Log() << Verbose(0) << "Saving as tremolo data input "); 1273 1379 if (configuration->SaveTREMOLO(filename, molecules)) 1274 Log() << Verbose(0) << "done." << endl;1380 DoLog(0) && (Log() << Verbose(0) << "done." << endl); 1275 1381 else 1276 Log() << Verbose(0) << "failed." << endl;1382 DoLog(0) && (Log() << Verbose(0) << "failed." << endl); 1277 1383 1278 1384 // translate each to its center and merge all molecules in MoleculeListClass into this molecule … … 1294 1400 } 1295 1401 1296 Log() << Verbose(0) << "Storing configuration ... " << endl;1402 DoLog(0) && (Log() << Verbose(0) << "Storing configuration ... " << endl); 1297 1403 // get correct valence orbitals 1298 1404 mol->CalculateOrbitals(*configuration); … … 1310 1416 output.close(); 1311 1417 output.clear(); 1312 Log() << Verbose(0) << "Saving of config file ";1418 DoLog(0) && (Log() << Verbose(0) << "Saving of config file "); 1313 1419 if (configuration->Save(filename, periode, mol)) 1314 Log() << Verbose(0) << "successful." << endl;1420 DoLog(0) && (Log() << Verbose(0) << "successful." << endl); 1315 1421 else 1316 Log() << Verbose(0) << "failed." << endl;1422 DoLog(0) && (Log() << Verbose(0) << "failed." << endl); 1317 1423 1318 1424 // and save to xyz file … … 1327 1433 output.open(filename, ios::trunc); 1328 1434 } 1329 Log() << Verbose(0) << "Saving of XYZ file ";1435 DoLog(0) && (Log() << Verbose(0) << "Saving of XYZ file "); 1330 1436 if (mol->MDSteps <= 1) { 1331 1437 if (mol->OutputXYZ(&output)) 1332 Log() << Verbose(0) << "successful." << endl;1438 DoLog(0) && (Log() << Verbose(0) << "successful." << endl); 1333 1439 else 1334 Log() << Verbose(0) << "failed." << endl;1440 DoLog(0) && (Log() << Verbose(0) << "failed." << endl); 1335 1441 } else { 1336 1442 if (mol->OutputTrajectoriesXYZ(&output)) 1337 Log() << Verbose(0) << "successful." << endl;1443 DoLog(0) && (Log() << Verbose(0) << "successful." << endl); 1338 1444 else 1339 Log() << Verbose(0) << "failed." << endl;1445 DoLog(0) && (Log() << Verbose(0) << "failed." << endl); 1340 1446 } 1341 1447 output.close(); … … 1347 1453 if (output == NULL) 1348 1454 strcpy(filename,"main_pcp_linux"); 1349 Log() << Verbose(0) << "Saving as mpqc input ";1455 DoLog(0) && (Log() << Verbose(0) << "Saving as mpqc input "); 1350 1456 if (configuration->SaveMPQC(filename, mol)) 1351 Log() << Verbose(0) << "done." << endl;1457 DoLog(0) && (Log() << Verbose(0) << "done." << endl); 1352 1458 else 1353 Log() << Verbose(0) << "failed." << endl;1459 DoLog(0) && (Log() << Verbose(0) << "failed." << endl); 1354 1460 1355 1461 if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) { 1356 eLog() << Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl;1462 DoeLog(2) && (eLog()<< Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl); 1357 1463 } 1358 1464 … … 1384 1490 enum ConfigStatus configPresent = absent; 1385 1491 clock_t start,end; 1492 double MaxDistance = -1; 1386 1493 int argptr; 1387 1494 molecule *mol = NULL; … … 1395 1502 do { 1396 1503 if (argv[argptr][0] == '-') { 1397 Log() << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n";1504 DoLog(0) && (Log() << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n"); 1398 1505 argptr++; 1399 1506 switch(argv[argptr-1][1]) { … … 1401 1508 case 'H': 1402 1509 case '?': 1403 Log() << Verbose(0) << "MoleCuilder suite" << endl << "==================" << endl << endl; 1404 Log() << Verbose(0) << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl; 1405 Log() << Verbose(0) << "or simply " << argv[0] << " without arguments for interactive session." << endl; 1406 Log() << Verbose(0) << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl; 1407 Log() << Verbose(0) << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl; 1408 Log() << Verbose(0) << "\t-b xx xy xz yy yz zz\tCenter atoms in domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl; 1409 Log() << Verbose(0) << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl; 1410 Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl; 1411 Log() << Verbose(0) << "\t-C\tPair Correlation analysis." << endl; 1412 Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl; 1413 Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl; 1414 Log() << Verbose(0) << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl; 1415 Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl; 1416 Log() << Verbose(0) << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl; 1417 Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl; 1418 Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl; 1419 Log() << Verbose(0) << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl; 1420 Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl; 1421 Log() << Verbose(0) << "\t-M <basis>\tSetting basis to store to MPQC config files." << endl; 1422 Log() << Verbose(0) << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl; 1423 Log() << Verbose(0) << "\t-N <radius> <file>\tGet non-convex-envelope." << endl; 1424 Log() << Verbose(0) << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl; 1425 Log() << Verbose(0) << "\t-O\tCenter atoms in origin." << endl; 1426 Log() << Verbose(0) << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl; 1427 Log() << Verbose(0) << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl; 1428 Log() << Verbose(0) << "\t-r <id>\t\tRemove an atom with given id." << endl; 1429 Log() << Verbose(0) << "\t-R <id> <radius>\t\tRemove all atoms out of sphere around a given one." << endl; 1430 Log() << Verbose(0) << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl; 1431 Log() << Verbose(0) << "\t-S <file> Store temperatures from the config file in <file>." << endl; 1432 Log() << Verbose(0) << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl; 1433 Log() << Verbose(0) << "\t-T x1 x2 x3\tTranslate periodically all atoms by this vector (x1,x2,x3)." << endl; 1434 Log() << Verbose(0) << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl; 1435 Log() << Verbose(0) << "\t-v\t\tsets verbosity (more is more)." << endl; 1436 Log() << Verbose(0) << "\t-V\t\tGives version information." << endl; 1437 Log() << Verbose(0) << "Note: config files must not begin with '-' !" << endl; 1510 DoLog(0) && (Log() << Verbose(0) << "MoleCuilder suite" << endl << "==================" << endl << endl); 1511 DoLog(0) && (Log() << Verbose(0) << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl); 1512 DoLog(0) && (Log() << Verbose(0) << "or simply " << argv[0] << " without arguments for interactive session." << endl); 1513 DoLog(0) && (Log() << Verbose(0) << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl); 1514 DoLog(0) && (Log() << Verbose(0) << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl); 1515 DoLog(0) && (Log() << Verbose(0) << "\t-b xx xy xz yy yz zz\tCenter atoms in domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl); 1516 DoLog(0) && (Log() << Verbose(0) << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl); 1517 DoLog(0) && (Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl); 1518 DoLog(0) && (Log() << Verbose(0) << "\t-C <type> [params] <output> <bin output> <BinWidth> <BinStart> <BinEnd>\tPair Correlation analysis." << endl); 1519 DoLog(0) && (Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl); 1520 DoLog(0) && (Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl); 1521 DoLog(0) && (Log() << Verbose(0) << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl); 1522 DoLog(0) && (Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl); 1523 DoLog(0) && (Log() << Verbose(0) << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl); 1524 DoLog(0) && (Log() << Verbose(0) << "\t-F <xyz of filler> <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl); 1525 DoLog(0) && (Log() << Verbose(0) << "\t-FF <MaxDistance> <xyz of filler> <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl); 1526 DoLog(0) && (Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl); 1527 DoLog(0) && (Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl); 1528 DoLog(0) && (Log() << Verbose(0) << "\t-I\t Dissect current system of molecules into a set of disconnected (subgraphs of) molecules." << endl); 1529 DoLog(0) && (Log() << Verbose(0) << "\t-j\t<path> Store all bonds to file." << endl); 1530 DoLog(0) && (Log() << Verbose(0) << "\t-J\t<path> Store adjacency per atom to file." << endl); 1531 DoLog(0) && (Log() << Verbose(0) << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl); 1532 DoLog(0) && (Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl); 1533 DoLog(0) && (Log() << Verbose(0) << "\t-M <basis>\tSetting basis to store to MPQC config files." << endl); 1534 DoLog(0) && (Log() << Verbose(0) << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl); 1535 DoLog(0) && (Log() << Verbose(0) << "\t-N <radius> <file>\tGet non-convex-envelope." << endl); 1536 DoLog(0) && (Log() << Verbose(0) << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl); 1537 DoLog(0) && (Log() << Verbose(0) << "\t-O\tCenter atoms in origin." << endl); 1538 DoLog(0) && (Log() << Verbose(0) << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl); 1539 DoLog(0) && (Log() << Verbose(0) << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl); 1540 DoLog(0) && (Log() << Verbose(0) << "\t-r <id>\t\tRemove an atom with given id." << endl); 1541 DoLog(0) && (Log() << Verbose(0) << "\t-R <id> <radius>\t\tRemove all atoms out of sphere around a given one." << endl); 1542 DoLog(0) && (Log() << Verbose(0) << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl); 1543 DoLog(0) && (Log() << Verbose(0) << "\t-S <file> Store temperatures from the config file in <file>." << endl); 1544 DoLog(0) && (Log() << Verbose(0) << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl); 1545 DoLog(0) && (Log() << Verbose(0) << "\t-T x1 x2 x3\tTranslate periodically all atoms by this vector (x1,x2,x3)." << endl); 1546 DoLog(0) && (Log() << Verbose(0) << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl); 1547 DoLog(0) && (Log() << Verbose(0) << "\t-v\t\tsets verbosity (more is more)." << endl); 1548 DoLog(0) && (Log() << Verbose(0) << "\t-V\t\tGives version information." << endl); 1549 DoLog(0) && (Log() << Verbose(0) << "\t-X\t\tset default name of a molecule." << endl); 1550 DoLog(0) && (Log() << Verbose(0) << "Note: config files must not begin with '-' !" << endl); 1438 1551 return (1); 1439 1552 break; … … 1443 1556 } 1444 1557 setVerbosity(verbosity); 1445 Log() << Verbose(0) << "Setting verbosity to " << verbosity << "." << endl;1558 DoLog(0) && (Log() << Verbose(0) << "Setting verbosity to " << verbosity << "." << endl); 1446 1559 break; 1447 1560 case 'V': 1448 Log() << Verbose(0) << argv[0] << " " << VERSIONSTRING << endl;1449 Log() << Verbose(0) << "Build your own molecule position set." << endl;1561 DoLog(0) && (Log() << Verbose(0) << argv[0] << " " << VERSIONSTRING << endl); 1562 DoLog(0) && (Log() << Verbose(0) << "Build your own molecule position set." << endl); 1450 1563 return (1); 1564 break; 1565 case 'B': 1566 if (ExitFlag == 0) ExitFlag = 1; 1567 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])) ) { 1568 ExitFlag = 255; 1569 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl); 1570 performCriticalExit(); 1571 } else { 1572 SaveFlag = true; 1573 j = -1; 1574 DoLog(1) && (Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl); 1575 double * const cell_size = World::get()->cell_size; 1576 for (int i=0;i<6;i++) { 1577 cell_size[i] = atof(argv[argptr+i]); 1578 } 1579 argptr+=6; 1580 } 1451 1581 break; 1452 1582 case 'e': 1453 1583 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1454 eLog() << Verbose(0) << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl;1584 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl); 1455 1585 performCriticalExit(); 1456 1586 } else { 1457 Log() << Verbose(0) << "Using " << argv[argptr] << " as elements database." << endl;1587 DoLog(0) && (Log() << Verbose(0) << "Using " << argv[argptr] << " as elements database." << endl); 1458 1588 strncpy (configuration.databasepath, argv[argptr], MAXSTRINGSIZE-1); 1459 1589 argptr+=1; … … 1462 1592 case 'g': 1463 1593 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1464 eLog() << Verbose(0) << "Not enough or invalid arguments for specifying bond length table: -g <table file>" << endl;1594 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying bond length table: -g <table file>" << endl); 1465 1595 performCriticalExit(); 1466 1596 } else { 1467 1597 BondGraphFileName = argv[argptr]; 1468 Log() << Verbose(0) << "Using " << BondGraphFileName << " as bond length table." << endl;1598 DoLog(0) && (Log() << Verbose(0) << "Using " << BondGraphFileName << " as bond length table." << endl); 1469 1599 argptr+=1; 1470 1600 } 1471 1601 break; 1472 1602 case 'n': 1473 Log() << Verbose(0) << "I won't parse trajectories." << endl;1603 DoLog(0) && (Log() << Verbose(0) << "I won't parse trajectories." << endl); 1474 1604 configuration.FastParsing = true; 1605 break; 1606 case 'X': 1607 { 1608 char **name = &(World::get()->DefaultName); 1609 delete[](*name); 1610 const int length = strlen(argv[argptr]); 1611 *name = new char[length+2]; 1612 strncpy(*name, argv[argptr], length); 1613 DoLog(0) && (Log() << Verbose(0) << "Default name of new molecules set to " << *name << "." << endl); 1614 } 1475 1615 break; 1476 1616 default: // no match? Step on … … 1484 1624 // 3a. Parse the element database 1485 1625 if (periode->LoadPeriodentafel(configuration.databasepath)) { 1486 Log() << Verbose(0) << "Element list loaded successfully." << endl;1626 DoLog(0) && (Log() << Verbose(0) << "Element list loaded successfully." << endl); 1487 1627 //periode->Output(); 1488 1628 } else { 1489 Log() << Verbose(0) << "Element list loading failed." << endl;1629 DoLog(0) && (Log() << Verbose(0) << "Element list loading failed." << endl); 1490 1630 return 1; 1491 1631 } … … 1493 1633 if (argv[1][0] != '-') { 1494 1634 // simply create a new molecule, wherein the config file is loaded and the manipulation takes place 1495 Log() << Verbose(0) << "Config file given." << endl;1635 DoLog(0) && (Log() << Verbose(0) << "Config file given." << endl); 1496 1636 test.open(argv[1], ios::in); 1497 1637 if (test == NULL) { … … 1499 1639 output.open(argv[1], ios::out); 1500 1640 if (output == NULL) { 1501 Log() << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl;1641 DoLog(1) && (Log() << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl); 1502 1642 configPresent = absent; 1503 1643 } else { 1504 Log() << Verbose(0) << "Empty configuration file." << endl;1644 DoLog(0) && (Log() << Verbose(0) << "Empty configuration file." << endl); 1505 1645 ConfigFileName = argv[1]; 1506 1646 configPresent = empty; … … 1510 1650 test.close(); 1511 1651 ConfigFileName = argv[1]; 1512 Log() << Verbose(1) << "Specified config file found, parsing ... ";1652 DoLog(1) && (Log() << Verbose(1) << "Specified config file found, parsing ... "); 1513 1653 switch (configuration.TestSyntax(ConfigFileName, periode)) { 1514 1654 case 1: 1515 Log() << Verbose(0) << "new syntax." << endl;1655 DoLog(0) && (Log() << Verbose(0) << "new syntax." << endl); 1516 1656 configuration.Load(ConfigFileName, BondGraphFileName, periode, molecules); 1517 1657 configPresent = present; 1518 1658 break; 1519 1659 case 0: 1520 Log() << Verbose(0) << "old syntax." << endl;1660 DoLog(0) && (Log() << Verbose(0) << "old syntax." << endl); 1521 1661 configuration.LoadOld(ConfigFileName, BondGraphFileName, periode, molecules); 1522 1662 configPresent = present; 1523 1663 break; 1524 1664 default: 1525 Log() << Verbose(0) << "Unknown syntax or empty, yet present file." << endl;1665 DoLog(0) && (Log() << Verbose(0) << "Unknown syntax or empty, yet present file." << endl); 1526 1666 configPresent = empty; 1527 1667 } … … 1546 1686 if (configuration.BG == NULL) { 1547 1687 configuration.BG = new BondGraph(configuration.GetIsAngstroem()); 1548 if (( BondGraphFileName.empty()) && (configuration.BG->LoadBondLengthTable(BondGraphFileName))) {1549 Log() << Verbose(0) << "Bond length table loaded successfully." << endl;1688 if ((!BondGraphFileName.empty()) && (configuration.BG->LoadBondLengthTable(BondGraphFileName))) { 1689 DoLog(0) && (Log() << Verbose(0) << "Bond length table loaded successfully." << endl); 1550 1690 } else { 1551 eLog() << Verbose(1) << "Bond length table loading failed." << endl;1691 DoeLog(1) && (eLog()<< Verbose(1) << "Bond length table loading failed." << endl); 1552 1692 } 1553 1693 } … … 1556 1696 argptr = 1; 1557 1697 do { 1558 Log() << Verbose(0) << "Current Command line argument: " << argv[argptr] << "." << endl;1698 DoLog(0) && (Log() << Verbose(0) << "Current Command line argument: " << argv[argptr] << "." << endl); 1559 1699 if (argv[argptr][0] == '-') { 1560 1700 argptr++; … … 1565 1705 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1566 1706 ExitFlag = 255; 1567 eLog() << Verbose(0) << "Not enough arguments for parsing: -p <xyz file>" << endl;1707 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough arguments for parsing: -p <xyz file>" << endl); 1568 1708 performCriticalExit(); 1569 1709 } else { 1570 1710 SaveFlag = true; 1571 Log() << Verbose(1) << "Parsing xyz file for new atoms." << endl;1711 DoLog(1) && (Log() << Verbose(1) << "Parsing xyz file for new atoms." << endl); 1572 1712 if (!mol->AddXYZFile(argv[argptr])) 1573 Log() << Verbose(2) << "File not found." << endl;1713 DoLog(2) && (Log() << Verbose(2) << "File not found." << endl); 1574 1714 else { 1575 Log() << Verbose(2) << "File found and parsed." << endl; 1576 // @TODO rather do the dissection afterwards 1577 // mol->SetNameFromFilename(argv[argptr]); 1578 // molecules->ListOfMolecules.remove(mol); 1579 // molecules->DissectMoleculeIntoConnectedSubgraphs(mol,&configuration); 1580 // delete(mol); 1581 // if (molecules->ListOfMolecules.size() != 0) { 1582 // for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) 1583 // if ((*ListRunner)->ActiveFlag) { 1584 // mol = *ListRunner; 1585 // break; 1586 // } 1587 // } 1715 DoLog(2) && (Log() << Verbose(2) << "File found and parsed." << endl); 1588 1716 configPresent = present; 1589 1717 } … … 1594 1722 if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3]))) { 1595 1723 ExitFlag = 255; 1596 eLog() << Verbose(0) << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl;1724 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl); 1597 1725 performCriticalExit(); 1598 1726 } else { 1599 1727 SaveFlag = true; 1600 Log() << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";1728 DoLog(1) && (Log() << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), "); 1601 1729 first = new atom; 1602 1730 first->type = periode->FindElement(atoi(argv[argptr])); 1603 1731 if (first->type != NULL) 1604 Log() << Verbose(2) << "found element " << first->type->name << endl;1732 DoLog(2) && (Log() << Verbose(2) << "found element " << first->type->name << endl); 1605 1733 for (int i=NDIM;i--;) 1606 1734 first->x.x[i] = atof(argv[argptr+1+i]); … … 1610 1738 configPresent = present; 1611 1739 } else 1612 eLog() << Verbose(1) << "Could not find the specified element." << endl;1740 DoeLog(1) && (eLog()<< Verbose(1) << "Could not find the specified element." << endl); 1613 1741 argptr+=4; 1614 1742 } … … 1623 1751 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1624 1752 ExitFlag = 255; 1625 eLog() << Verbose(0) << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl;1753 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl); 1626 1754 performCriticalExit(); 1627 1755 } else { 1628 1756 configuration.basis = argv[argptr]; 1629 Log() << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl;1757 DoLog(1) && (Log() << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl); 1630 1758 argptr+=1; 1631 1759 } … … 1634 1762 if (ExitFlag == 0) ExitFlag = 1; 1635 1763 { 1636 Log() << Verbose(1) << "Depth-First-Search Analysis." << endl;1764 DoLog(1) && (Log() << Verbose(1) << "Depth-First-Search Analysis." << endl); 1637 1765 MoleculeLeafClass *Subgraphs = NULL; // list of subgraphs from DFS analysis 1638 1766 int *MinimumRingSize = new int[mol->AtomCount]; … … 1664 1792 //argptr+=1; 1665 1793 break; 1794 case 'I': 1795 DoLog(1) && (Log() << Verbose(1) << "Dissecting molecular system into a set of disconnected subgraphs ... " << endl); 1796 // @TODO rather do the dissection afterwards 1797 molecules->DissectMoleculeIntoConnectedSubgraphs(periode, &configuration); 1798 mol = NULL; 1799 if (molecules->ListOfMolecules.size() != 0) { 1800 for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) 1801 if ((*ListRunner)->ActiveFlag) { 1802 mol = *ListRunner; 1803 break; 1804 } 1805 } 1806 if ((mol == NULL) && (!molecules->ListOfMolecules.empty())) { 1807 mol = *(molecules->ListOfMolecules.begin()); 1808 if (mol != NULL) 1809 mol->ActiveFlag = true; 1810 } 1811 break; 1666 1812 case 'C': 1667 if (ExitFlag == 0) ExitFlag = 1; 1668 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr][0] == '-') || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-')) { 1669 ExitFlag = 255; 1670 eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C <Z> <output> <bin output>" << endl; 1671 performCriticalExit(); 1672 } else { 1673 SaveFlag = false; 1674 ofstream output(argv[argptr+1]); 1675 ofstream binoutput(argv[argptr+2]); 1676 const double radius = 5.; 1677 1678 // get the boundary 1679 class molecule *Boundary = NULL; 1680 class Tesselation *TesselStruct = NULL; 1681 const LinkedCell *LCList = NULL; 1682 // find biggest molecule 1683 int counter = 0; 1684 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) { 1685 if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) { 1686 Boundary = *BigFinder; 1813 { 1814 int ranges[3] = {1, 1, 1}; 1815 bool periodic = (argv[argptr-1][2] =='p'); 1816 if (ExitFlag == 0) ExitFlag = 1; 1817 if ((argptr >= argc)) { 1818 ExitFlag = 255; 1819 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); 1820 performCriticalExit(); 1821 } else { 1822 switch(argv[argptr][0]) { 1823 case 'E': 1824 { 1825 if ((argptr+6 >= argc) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (!IsValidNumber(argv[argptr+2])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-')) { 1826 ExitFlag = 255; 1827 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C E <Z1> <Z2> <output> <bin output>" << endl); 1828 performCriticalExit(); 1829 } else { 1830 ofstream output(argv[argptr+3]); 1831 ofstream binoutput(argv[argptr+4]); 1832 const double BinStart = atof(argv[argptr+5]); 1833 const double BinEnd = atof(argv[argptr+6]); 1834 1835 element *elemental = periode->FindElement((const int) atoi(argv[argptr+1])); 1836 element *elemental2 = periode->FindElement((const int) atoi(argv[argptr+2])); 1837 PairCorrelationMap *correlationmap = NULL; 1838 if (periodic) 1839 correlationmap = PeriodicPairCorrelation(molecules, elemental, elemental2, ranges); 1840 else 1841 correlationmap = PairCorrelation(molecules, elemental, elemental2); 1842 //OutputCorrelationToSurface(&output, correlationmap); 1843 BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd ); 1844 OutputCorrelation ( &binoutput, binmap ); 1845 output.close(); 1846 binoutput.close(); 1847 delete(binmap); 1848 delete(correlationmap); 1849 argptr+=7; 1850 } 1851 } 1852 break; 1853 1854 case 'P': 1855 { 1856 if ((argptr+8 >= argc) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+7])) || (!IsValidNumber(argv[argptr+8])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-') || (argv[argptr+5][0] == '-') || (argv[argptr+6][0] == '-')) { 1857 ExitFlag = 255; 1858 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C P <Z1> <x> <y> <z> <output> <bin output>" << endl); 1859 performCriticalExit(); 1860 } else { 1861 ofstream output(argv[argptr+5]); 1862 ofstream binoutput(argv[argptr+6]); 1863 const double BinStart = atof(argv[argptr+7]); 1864 const double BinEnd = atof(argv[argptr+8]); 1865 1866 element *elemental = periode->FindElement((const int) atoi(argv[argptr+1])); 1867 Vector *Point = new Vector((const double) atof(argv[argptr+1]),(const double) atof(argv[argptr+2]),(const double) atof(argv[argptr+3])); 1868 CorrelationToPointMap *correlationmap = NULL; 1869 if (periodic) 1870 correlationmap = PeriodicCorrelationToPoint(molecules, elemental, Point, ranges); 1871 else 1872 correlationmap = CorrelationToPoint(molecules, elemental, Point); 1873 //OutputCorrelationToSurface(&output, correlationmap); 1874 BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd ); 1875 OutputCorrelation ( &binoutput, binmap ); 1876 output.close(); 1877 binoutput.close(); 1878 delete(Point); 1879 delete(binmap); 1880 delete(correlationmap); 1881 argptr+=9; 1882 } 1883 } 1884 break; 1885 1886 case 'S': 1887 { 1888 if ((argptr+6 >= argc) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-')) { 1889 ExitFlag = 255; 1890 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C S <Z> <output> <bin output> <BinWidth> <BinStart> <BinEnd>" << endl); 1891 performCriticalExit(); 1892 } else { 1893 ofstream output(argv[argptr+2]); 1894 ofstream binoutput(argv[argptr+3]); 1895 const double radius = 4.; 1896 const double BinWidth = atof(argv[argptr+4]); 1897 const double BinStart = atof(argv[argptr+5]); 1898 const double BinEnd = atof(argv[argptr+6]); 1899 double LCWidth = 20.; 1900 if (BinEnd > 0) { 1901 if (BinEnd > 2.*radius) 1902 LCWidth = BinEnd; 1903 else 1904 LCWidth = 2.*radius; 1905 } 1906 1907 // get the boundary 1908 class molecule *Boundary = NULL; 1909 class Tesselation *TesselStruct = NULL; 1910 const LinkedCell *LCList = NULL; 1911 // find biggest molecule 1912 int counter = 0; 1913 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) { 1914 if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) { 1915 Boundary = *BigFinder; 1916 } 1917 counter++; 1918 } 1919 bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives"); 1920 counter = 0; 1921 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) { 1922 Actives[counter++] = (*BigFinder)->ActiveFlag; 1923 (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true; 1924 } 1925 LCList = new LinkedCell(Boundary, LCWidth); 1926 element *elemental = periode->FindElement((const int) atoi(argv[argptr+1])); 1927 FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL); 1928 CorrelationToSurfaceMap *surfacemap = NULL; 1929 if (periodic) 1930 surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges); 1931 else 1932 surfacemap = CorrelationToSurface( molecules, elemental, TesselStruct, LCList); 1933 OutputCorrelationToSurface(&output, surfacemap); 1934 // check whether radius was appropriate 1935 { 1936 double start; double end; 1937 GetMinMax( surfacemap, start, end); 1938 if (LCWidth < end) 1939 DoeLog(1) && (eLog()<< Verbose(1) << "Linked Cell width is smaller than the found range of values! Bins can only be correct up to: " << radius << "." << endl); 1940 } 1941 BinPairMap *binmap = BinData( surfacemap, BinWidth, BinStart, BinEnd ); 1942 OutputCorrelation ( &binoutput, binmap ); 1943 output.close(); 1944 binoutput.close(); 1945 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) 1946 (*BigFinder)->ActiveFlag = Actives[counter++]; 1947 Free(&Actives); 1948 delete(LCList); 1949 delete(TesselStruct); 1950 delete(binmap); 1951 delete(surfacemap); 1952 argptr+=7; 1953 } 1954 } 1955 break; 1956 1957 default: 1958 ExitFlag = 255; 1959 DoeLog(0) && (eLog()<< Verbose(0) << "Invalid type given for pair correlation analysis: -C <type: E/P/S> [more params] <output> <bin output>" << endl); 1960 performCriticalExit(); 1961 break; 1687 1962 } 1688 counter++;1689 1963 } 1690 bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives"); 1691 counter = 0; 1692 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) { 1693 Actives[counter] = (*BigFinder)->ActiveFlag; 1694 (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true; 1695 } 1696 LCList = new LinkedCell(Boundary, 2.*radius); 1697 element *elemental = periode->FindElement((const int) atoi(argv[argptr])); 1698 FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL); 1699 int ranges[NDIM] = {1,1,1}; 1700 CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges ); 1701 BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 0. ); 1702 OutputCorrelation ( &binoutput, binmap ); 1703 output.close(); 1704 binoutput.close(); 1705 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) 1706 (*BigFinder)->ActiveFlag = Actives[counter]; 1707 Free(&Actives); 1708 delete(LCList); 1709 delete(TesselStruct); 1710 argptr+=3; 1711 } 1712 break; 1964 break; 1965 } 1713 1966 case 'E': 1714 1967 if (ExitFlag == 0) ExitFlag = 1; 1715 1968 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) { 1716 1969 ExitFlag = 255; 1717 eLog() << Verbose(0) << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl;1970 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl); 1718 1971 performCriticalExit(); 1719 1972 } else { 1720 1973 SaveFlag = true; 1721 Log() << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl;1974 DoLog(1) && (Log() << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl); 1722 1975 first = mol->FindAtom(atoi(argv[argptr])); 1723 1976 first->type = periode->FindElement(atoi(argv[argptr+1])); … … 1727 1980 case 'F': 1728 1981 if (ExitFlag == 0) ExitFlag = 1; 1729 if (argptr+5 >=argc) { 1982 MaxDistance = -1; 1983 if (argv[argptr-1][2] == 'F') { // option is -FF? 1984 // fetch first argument as max distance to surface 1985 MaxDistance = atof(argv[argptr++]); 1986 DoLog(0) && (Log() << Verbose(0) << "Filling with maximum layer distance of " << MaxDistance << "." << endl); 1987 } 1988 if ((argptr+7 >=argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (!IsValidNumber(argv[argptr+7]))) { 1730 1989 ExitFlag = 255; 1731 eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <randatom> <randmol> <DoRotate>" << endl;1990 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <xyz of filler> <dist_x> <dist_y> <dist_z> <boundary> <randatom> <randmol> <DoRotate>" << endl); 1732 1991 performCriticalExit(); 1733 1992 } else { 1734 1993 SaveFlag = true; 1735 Log() << Verbose(1) << "Filling Box with water molecules." << endl;1994 DoLog(1) && (Log() << Verbose(1) << "Filling Box with water molecules." << endl); 1736 1995 // construct water molecule 1737 molecule *filler = new molecule(periode);; 1996 molecule *filler = new molecule(periode); 1997 if (!filler->AddXYZFile(argv[argptr])) { 1998 DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse filler molecule from " << argv[argptr] << "." << endl); 1999 } 2000 filler->SetNameFromFilename(argv[argptr]); 2001 configuration.BG->ConstructBondGraph(filler); 1738 2002 molecule *Filling = NULL; 1739 atom *second = NULL, *third = NULL;1740 first = new atom();1741 first->type = periode->FindElement(1);1742 first->x.Init(0.441, -0.143, 0.);1743 filler->AddAtom(first);1744 second = new atom();1745 second->type = periode->FindElement(1);1746 second->x.Init(-0.464, 1.137, 0.0);1747 filler->AddAtom(second);1748 third = new atom();1749 third->type = periode->FindElement(8);1750 third->x.Init(-0.464, 0.177, 0.);1751 filler->AddAtom(third);1752 filler->AddBond(first, third, 1);1753 filler->AddBond(second, third, 1);1754 2003 // call routine 1755 2004 double distance[NDIM]; 1756 2005 for (int i=0;i<NDIM;i++) 1757 distance[i] = atof(argv[argptr+i ]);1758 Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atoi(argv[argptr+5]));2006 distance[i] = atof(argv[argptr+i+1]); 2007 Filling = FillBoxWithMolecule(molecules, filler, configuration, MaxDistance, distance, atof(argv[argptr+4]), atof(argv[argptr+5]), atof(argv[argptr+6]), atoi(argv[argptr+7])); 1759 2008 if (Filling != NULL) { 2009 Filling->ActiveFlag = false; 1760 2010 molecules->insert(Filling); 1761 2011 } … … 1768 2018 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1769 2019 ExitFlag =255; 1770 eLog() << Verbose(0) << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;2020 DoeLog(0) && (eLog()<< Verbose(0) << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl); 1771 2021 performCriticalExit(); 1772 2022 } else { 1773 Log() << Verbose(0) << "Parsing bonds from " << argv[argptr] << "." << endl;2023 DoLog(0) && (Log() << Verbose(0) << "Parsing bonds from " << argv[argptr] << "." << endl); 1774 2024 ifstream *input = new ifstream(argv[argptr]); 1775 2025 mol->CreateAdjacencyListFromDbondFile(input); … … 1778 2028 } 1779 2029 break; 2030 2031 case 'J': 2032 if (ExitFlag == 0) ExitFlag = 1; 2033 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 2034 ExitFlag =255; 2035 DoeLog(0) && (eLog()<< Verbose(0) << "Missing path of adjacency file: -j <path>" << endl); 2036 performCriticalExit(); 2037 } else { 2038 DoLog(0) && (Log() << Verbose(0) << "Storing adjacency to path " << argv[argptr] << "." << endl); 2039 configuration.BG->ConstructBondGraph(mol); 2040 mol->StoreAdjacencyToFile(NULL, argv[argptr]); 2041 argptr+=1; 2042 } 2043 break; 2044 2045 case 'j': 2046 if (ExitFlag == 0) ExitFlag = 1; 2047 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 2048 ExitFlag =255; 2049 DoeLog(0) && (eLog()<< Verbose(0) << "Missing path of bonds file: -j <path>" << endl); 2050 performCriticalExit(); 2051 } else { 2052 DoLog(0) && (Log() << Verbose(0) << "Storing bonds to path " << argv[argptr] << "." << endl); 2053 configuration.BG->ConstructBondGraph(mol); 2054 mol->StoreBondsToFile(NULL, argv[argptr]); 2055 argptr+=1; 2056 } 2057 break; 2058 1780 2059 case 'N': 1781 2060 if (ExitFlag == 0) ExitFlag = 1; 1782 2061 if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){ 1783 2062 ExitFlag = 255; 1784 eLog() << Verbose(0) << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;2063 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl); 1785 2064 performCriticalExit(); 1786 2065 } else { … … 1790 2069 //string filename(argv[argptr+1]); 1791 2070 //filename.append(".csv"); 1792 Log() << Verbose(0) << "Evaluating non-convex envelope of biggest molecule.";1793 Log() << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;2071 DoLog(0) && (Log() << Verbose(0) << "Evaluating non-convex envelope of biggest molecule."); 2072 DoLog(1) && (Log() << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl); 1794 2073 // find biggest molecule 1795 2074 int counter = 0; … … 1801 2080 counter++; 1802 2081 } 1803 Log() << Verbose(1) << "Biggest molecule has " << Boundary->AtomCount << " atoms." << endl;2082 DoLog(1) && (Log() << Verbose(1) << "Biggest molecule has " << Boundary->AtomCount << " atoms." << endl); 1804 2083 start = clock(); 1805 2084 LCList = new LinkedCell(Boundary, atof(argv[argptr])*2.); … … 1808 2087 //FindDistributionOfEllipsoids(T, &LCList, N, number, filename.c_str()); 1809 2088 end = clock(); 1810 Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;2089 DoLog(0) && (Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl); 1811 2090 delete(LCList); 1812 2091 delete(T); … … 1818 2097 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1819 2098 ExitFlag = 255; 1820 eLog() << Verbose(0) << "Not enough or invalid arguments given for storing tempature: -S <temperature file>" << endl;2099 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for storing tempature: -S <temperature file>" << endl); 1821 2100 performCriticalExit(); 1822 2101 } else { 1823 Log() << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;2102 DoLog(1) && (Log() << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl); 1824 2103 ofstream *output = new ofstream(argv[argptr], ios::trunc); 1825 2104 if (!mol->OutputTemperatureFromTrajectories(output, 0, mol->MDSteps)) 1826 Log() << Verbose(2) << "File could not be written." << endl;2105 DoLog(2) && (Log() << Verbose(2) << "File could not be written." << endl); 1827 2106 else 1828 Log() << Verbose(2) << "File stored." << endl;2107 DoLog(2) && (Log() << Verbose(2) << "File stored." << endl); 1829 2108 output->close(); 1830 2109 delete(output); … … 1836 2115 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1837 2116 ExitFlag = 255; 1838 eLog() << Verbose(0) << "Not enough or invalid arguments given for storing tempature: -L <step0> <step1> <prefix> <identity mapping?>" << endl;2117 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for storing tempature: -L <step0> <step1> <prefix> <identity mapping?>" << endl); 1839 2118 performCriticalExit(); 1840 2119 } else { 1841 2120 SaveFlag = true; 1842 Log() << Verbose(1) << "Linear interpolation between configuration " << argv[argptr] << " and " << argv[argptr+1] << "." << endl;2121 DoLog(1) && (Log() << Verbose(1) << "Linear interpolation between configuration " << argv[argptr] << " and " << argv[argptr+1] << "." << endl); 1843 2122 if (atoi(argv[argptr+3]) == 1) 1844 Log() << Verbose(1) << "Using Identity for the permutation map." << endl;2123 DoLog(1) && (Log() << Verbose(1) << "Using Identity for the permutation map." << endl); 1845 2124 if (!mol->LinearInterpolationBetweenConfiguration(atoi(argv[argptr]), atoi(argv[argptr+1]), argv[argptr+2], configuration, atoi(argv[argptr+3])) == 1 ? true : false) 1846 Log() << Verbose(2) << "Could not store " << argv[argptr+2] << " files." << endl;2125 DoLog(2) && (Log() << Verbose(2) << "Could not store " << argv[argptr+2] << " files." << endl); 1847 2126 else 1848 Log() << Verbose(2) << "Steps created and " << argv[argptr+2] << " files stored." << endl;2127 DoLog(2) && (Log() << Verbose(2) << "Steps created and " << argv[argptr+2] << " files stored." << endl); 1849 2128 argptr+=4; 1850 2129 } … … 1854 2133 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1855 2134 ExitFlag = 255; 1856 eLog() << Verbose(0) << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl;2135 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl); 1857 2136 performCriticalExit(); 1858 2137 } else { 1859 2138 SaveFlag = true; 1860 Log() << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;2139 DoLog(1) && (Log() << Verbose(1) << "Parsing forces file and Verlet integrating." << endl); 1861 2140 if (!mol->VerletForceIntegration(argv[argptr], configuration)) 1862 Log() << Verbose(2) << "File not found." << endl;2141 DoLog(2) && (Log() << Verbose(2) << "File not found." << endl); 1863 2142 else 1864 Log() << Verbose(2) << "File found and parsed." << endl;2143 DoLog(2) && (Log() << Verbose(2) << "File found and parsed." << endl); 1865 2144 argptr+=1; 1866 2145 } … … 1870 2149 if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) { 1871 2150 ExitFlag = 255; 1872 eLog() << Verbose(0) << "Not enough or invalid arguments given for removing atoms: -R <id> <distance>" << endl;2151 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -R <id> <distance>" << endl); 1873 2152 performCriticalExit(); 1874 2153 } else { 1875 2154 SaveFlag = true; 1876 Log() << Verbose(1) << "Removing atoms around " << argv[argptr] << " with radius " << argv[argptr+1] << "." << endl;2155 DoLog(1) && (Log() << Verbose(1) << "Removing atoms around " << argv[argptr] << " with radius " << argv[argptr+1] << "." << endl); 1877 2156 double tmp1 = atof(argv[argptr+1]); 1878 2157 atom *third = mol->FindAtom(atoi(argv[argptr])); … … 1887 2166 } 1888 2167 } else { 1889 eLog() << Verbose(1) << "Removal failed due to missing atoms on molecule or wrong id." << endl;2168 DoeLog(1) && (eLog()<< Verbose(1) << "Removal failed due to missing atoms on molecule or wrong id." << endl); 1890 2169 } 1891 2170 argptr+=2; … … 1896 2175 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1897 2176 ExitFlag = 255; 1898 eLog() << Verbose(0) << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl;2177 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl); 1899 2178 performCriticalExit(); 1900 2179 } else { 1901 2180 if (ExitFlag == 0) ExitFlag = 1; 1902 2181 SaveFlag = true; 1903 Log() << Verbose(1) << "Translating all ions by given vector." << endl;2182 DoLog(1) && (Log() << Verbose(1) << "Translating all ions by given vector." << endl); 1904 2183 for (int i=NDIM;i--;) 1905 2184 x.x[i] = atof(argv[argptr+i]); … … 1912 2191 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1913 2192 ExitFlag = 255; 1914 eLog() << Verbose(0) << "Not enough or invalid arguments given for periodic translation: -T <x> <y> <z>" << endl;2193 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for periodic translation: -T <x> <y> <z>" << endl); 1915 2194 performCriticalExit(); 1916 2195 } else { 1917 2196 if (ExitFlag == 0) ExitFlag = 1; 1918 2197 SaveFlag = true; 1919 Log() << Verbose(1) << "Translating all ions periodically by given vector." << endl;2198 DoLog(1) && (Log() << Verbose(1) << "Translating all ions periodically by given vector." << endl); 1920 2199 for (int i=NDIM;i--;) 1921 2200 x.x[i] = atof(argv[argptr+i]); … … 1928 2207 if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1929 2208 ExitFlag = 255; 1930 eLog() << Verbose(0) << "Not enough or invalid arguments given for scaling: -s <factor_x> [factor_y] [factor_z]" << endl;2209 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for scaling: -s <factor_x> [factor_y] [factor_z]" << endl); 1931 2210 performCriticalExit(); 1932 2211 } else { 1933 2212 SaveFlag = true; 1934 2213 j = -1; 1935 Log() << Verbose(1) << "Scaling all ion positions by factor." << endl;2214 DoLog(1) && (Log() << Verbose(1) << "Scaling all ion positions by factor." << endl); 1936 2215 factor = new double[NDIM]; 1937 2216 factor[0] = atof(argv[argptr]); … … 1939 2218 factor[2] = atof(argv[argptr+2]); 1940 2219 mol->Scale((const double ** const)&factor); 2220 double * const cell_size = World::get()->cell_size; 1941 2221 for (int i=0;i<NDIM;i++) { 1942 2222 j += i+1; 1943 2223 x.x[i] = atof(argv[NDIM+i]); 1944 mol->cell_size[j]*=factor[i];2224 cell_size[j]*=factor[i]; 1945 2225 } 1946 2226 delete[](factor); … … 1952 2232 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])) ) { 1953 2233 ExitFlag = 255; 1954 eLog() << Verbose(0) << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl;2234 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl); 1955 2235 performCriticalExit(); 1956 2236 } else { 1957 2237 SaveFlag = true; 1958 2238 j = -1; 1959 Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl; 2239 DoLog(1) && (Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl); 2240 double * const cell_size = World::get()->cell_size; 1960 2241 for (int i=0;i<6;i++) { 1961 mol->cell_size[i] = atof(argv[argptr+i]);2242 cell_size[i] = atof(argv[argptr+i]); 1962 2243 } 1963 2244 // center … … 1970 2251 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])) ) { 1971 2252 ExitFlag = 255; 1972 eLog() << Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl;2253 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl); 1973 2254 performCriticalExit(); 1974 2255 } else { 1975 2256 SaveFlag = true; 1976 2257 j = -1; 1977 Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl; 2258 DoLog(1) && (Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl); 2259 double * const cell_size = World::get()->cell_size; 1978 2260 for (int i=0;i<6;i++) { 1979 mol->cell_size[i] = atof(argv[argptr+i]);2261 cell_size[i] = atof(argv[argptr+i]); 1980 2262 } 1981 2263 // center … … 1988 2270 if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1989 2271 ExitFlag = 255; 1990 eLog() << Verbose(0) << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl;2272 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl); 1991 2273 performCriticalExit(); 1992 2274 } else { 1993 2275 SaveFlag = true; 1994 2276 j = -1; 1995 Log() << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;2277 DoLog(1) && (Log() << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl); 1996 2278 // make every coordinate positive 1997 2279 mol->CenterEdge(&x); … … 1999 2281 mol->SetBoxDimension(&x); 2000 2282 // translate each coordinate by boundary 2283 double * const cell_size = World::get()->cell_size; 2001 2284 j=-1; 2002 2285 for (int i=0;i<NDIM;i++) { 2003 2286 j += i+1; 2004 2287 x.x[i] = atof(argv[argptr+i]); 2005 mol->cell_size[j] += x.x[i]*2.;2288 cell_size[j] += x.x[i]*2.; 2006 2289 } 2007 2290 mol->Translate((const Vector *)&x); … … 2012 2295 if (ExitFlag == 0) ExitFlag = 1; 2013 2296 SaveFlag = true; 2014 Log() << Verbose(1) << "Centering atoms on edge and setting box dimensions." << endl;2297 DoLog(1) && (Log() << Verbose(1) << "Centering atoms on edge and setting box dimensions." << endl); 2015 2298 x.Zero(); 2016 2299 mol->CenterEdge(&x); … … 2022 2305 if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr]))) { 2023 2306 ExitFlag = 255; 2024 eLog() << Verbose(0) << "Not enough or invalid arguments given for removing atoms: -r <id>" << endl;2307 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -r <id>" << endl); 2025 2308 performCriticalExit(); 2026 2309 } else { 2027 2310 SaveFlag = true; 2028 Log() << Verbose(1) << "Removing atom " << argv[argptr] << "." << endl;2311 DoLog(1) && (Log() << Verbose(1) << "Removing atom " << argv[argptr] << "." << endl); 2029 2312 atom *first = mol->FindAtom(atoi(argv[argptr])); 2030 2313 mol->RemoveAtom(first); … … 2036 2319 if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) { 2037 2320 ExitFlag = 255; 2038 eLog() << Verbose(0) << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl;2321 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl); 2039 2322 performCriticalExit(); 2040 2323 } else { 2041 Log() << Verbose(0) << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;2042 Log() << Verbose(0) << "Creating connection matrix..." << endl;2324 DoLog(0) && (Log() << Verbose(0) << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl); 2325 DoLog(0) && (Log() << Verbose(0) << "Creating connection matrix..." << endl); 2043 2326 start = clock(); 2044 2327 mol->CreateAdjacencyList(atof(argv[argptr++]), configuration.GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL); 2045 Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;2328 DoLog(0) && (Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl); 2046 2329 if (mol->first->next != mol->last) { 2047 2330 ExitFlag = mol->FragmentMolecule(atoi(argv[argptr]), &configuration); 2048 2331 } 2049 2332 end = clock(); 2050 Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;2333 DoLog(0) && (Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl); 2051 2334 argptr+=2; 2052 2335 } … … 2056 2339 j = atoi(argv[argptr++]); 2057 2340 if ((j<0) || (j>1)) { 2058 eLog() << Verbose(1) << "Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;2341 DoeLog(1) && (eLog()<< Verbose(1) << "Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl); 2059 2342 j = 0; 2060 2343 } 2061 2344 if (j) { 2062 2345 SaveFlag = true; 2063 Log() << Verbose(0) << "Converting to prinicipal axis system." << endl;2346 DoLog(0) && (Log() << Verbose(0) << "Converting to prinicipal axis system." << endl); 2064 2347 } else 2065 Log() << Verbose(0) << "Evaluating prinicipal axis." << endl;2348 DoLog(0) && (Log() << Verbose(0) << "Evaluating prinicipal axis." << endl); 2066 2349 mol->PrincipalAxisSystem((bool)j); 2067 2350 break; … … 2070 2353 if ((argptr+1 >= argc) || (argv[argptr][0] == '-')){ 2071 2354 ExitFlag = 255; 2072 eLog() << Verbose(0) << "Not enough or invalid arguments given for convex envelope: -o <convex output file> <non-convex output file>" << endl;2355 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for convex envelope: -o <convex output file> <non-convex output file>" << endl); 2073 2356 performCriticalExit(); 2074 2357 } else { 2075 2358 class Tesselation *TesselStruct = NULL; 2076 2359 const LinkedCell *LCList = NULL; 2077 Log() << Verbose(0) << "Evaluating volume of the convex envelope.";2078 Log() << Verbose(1) << "Storing tecplot convex data in " << argv[argptr] << "." << endl;2079 Log() << Verbose(1) << "Storing tecplot non-convex data in " << argv[argptr+1] << "." << endl;2360 DoLog(0) && (Log() << Verbose(0) << "Evaluating volume of the convex envelope."); 2361 DoLog(1) && (Log() << Verbose(1) << "Storing tecplot convex data in " << argv[argptr] << "." << endl); 2362 DoLog(1) && (Log() << Verbose(1) << "Storing tecplot non-convex data in " << argv[argptr+1] << "." << endl); 2080 2363 LCList = new LinkedCell(mol, 10.); 2081 2364 //FindConvexBorder(mol, LCList, argv[argptr]); … … 2084 2367 double volumedifference = ConvexizeNonconvexEnvelope(TesselStruct, mol, argv[argptr]); 2085 2368 double clustervolume = VolumeOfConvexEnvelope(TesselStruct, &configuration); 2086 Log() << Verbose(0) << "The tesselated volume area is " << clustervolume << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;2087 Log() << Verbose(0) << "The non-convex tesselated volume area is " << clustervolume-volumedifference << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;2369 DoLog(0) && (Log() << Verbose(0) << "The tesselated volume area is " << clustervolume << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl); 2370 DoLog(0) && (Log() << Verbose(0) << "The non-convex tesselated volume area is " << clustervolume-volumedifference << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl); 2088 2371 delete(TesselStruct); 2089 2372 delete(LCList); … … 2095 2378 if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) { 2096 2379 ExitFlag = 255; 2097 eLog() << Verbose(0) << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl;2380 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl); 2098 2381 performCriticalExit(); 2099 2382 } else { 2100 2383 volume = atof(argv[argptr++]); 2101 Log() << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;2384 DoLog(0) && (Log() << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl); 2102 2385 } 2103 2386 case 'u': … … 2106 2389 if (volume != -1) 2107 2390 ExitFlag = 255; 2108 eLog() << Verbose(0) << "Not enough arguments given for suspension: -u <density>" << endl;2391 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl); 2109 2392 performCriticalExit(); 2110 2393 } else { 2111 2394 double density; 2112 2395 SaveFlag = true; 2113 Log() << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";2396 DoLog(0) && (Log() << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water."); 2114 2397 density = atof(argv[argptr++]); 2115 2398 if (density < 1.0) { 2116 eLog() << Verbose(1) << "Density must be greater than 1.0g/cm^3 !" << endl;2399 DoeLog(1) && (eLog()<< Verbose(1) << "Density must be greater than 1.0g/cm^3 !" << endl); 2117 2400 density = 1.3; 2118 2401 } … … 2120 2403 // repetition[i] = atoi(argv[argptr++]); 2121 2404 // if (repetition[i] < 1) 2122 // eLog() << Verbose(1) << "repetition value must be greater 1!" << endl;2405 // DoeLog(1) && (eLog()<< Verbose(1) << "repetition value must be greater 1!" << endl); 2123 2406 // repetition[i] = 1; 2124 2407 // } … … 2130 2413 if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 2131 2414 ExitFlag = 255; 2132 eLog() << Verbose(0) << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl;2415 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl); 2133 2416 performCriticalExit(); 2134 2417 } else { 2135 2418 SaveFlag = true; 2419 double * const cell_size = World::get()->cell_size; 2136 2420 for (int axis = 1; axis <= NDIM; axis++) { 2137 2421 int faktor = atoi(argv[argptr++]); … … 2140 2424 Vector ** vectors; 2141 2425 if (faktor < 1) { 2142 eLog() << Verbose(1) << "Repetition factor mus be greater than 1!" << endl;2426 DoeLog(1) && (eLog()<< Verbose(1) << "Repetition factor mus be greater than 1!" << endl); 2143 2427 faktor = 1; 2144 2428 } … … 2157 2441 } 2158 2442 if (count != j) 2159 eLog() << Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;2443 DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl); 2160 2444 x.Zero(); 2161 2445 y.Zero(); 2162 y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude2446 y.x[abs(axis)-1] = cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude 2163 2447 for (int i=1;i<faktor;i++) { // then add this list with respective translation factor times 2164 2448 x.AddVector(&y); // per factor one cell width further … … 2181 2465 mol->Translate(&x); 2182 2466 } 2183 mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;2467 cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor; 2184 2468 } 2185 2469 } … … 2198 2482 } else { // no arguments, hence scan the elements db 2199 2483 if (periode->LoadPeriodentafel(configuration.databasepath)) 2200 Log() << Verbose(0) << "Element list loaded successfully." << endl;2484 DoLog(0) && (Log() << Verbose(0) << "Element list loaded successfully." << endl); 2201 2485 else 2202 Log() << Verbose(0) << "Element list loading failed." << endl;2486 DoLog(0) && (Log() << Verbose(0) << "Element list loading failed." << endl); 2203 2487 configuration.RetrieveConfigPathAndName("main_pcp_linux"); 2204 2488 } … … 2221 2505 char *ConfigFileName = NULL; 2222 2506 int j; 2507 2508 cout << ESPACKVersion << endl; 2509 2510 DoLog(1) && (Log() << Verbose(1) << "test" << endl); 2511 DoLog(3) && (Log() << Verbose(1) << "test"); 2223 2512 2224 2513 setVerbosity(0); … … 2248 2537 if (molecules->ListOfMolecules.size() == 0) { 2249 2538 mol = new molecule(periode); 2250 if (mol->cell_size[0] == 0.) { 2251 Log() << Verbose(0) << "enter lower tridiagonal form of basis matrix" << endl << endl; 2539 double * const cell_size = World::get()->cell_size; 2540 if (cell_size[0] == 0.) { 2541 DoLog(0) && (Log() << Verbose(0) << "enter lower tridiagonal form of basis matrix" << endl << endl); 2252 2542 for (int i=0;i<6;i++) { 2253 Log() << Verbose(1) << "Cell size" << i << ": ";2254 cin >> mol->cell_size[i];2543 DoLog(1) && (Log() << Verbose(1) << "Cell size" << i << ": "); 2544 cin >> cell_size[i]; 2255 2545 } 2256 2546 } … … 2262 2552 2263 2553 // now the main construction loop 2264 Log() << Verbose(0) << endl << "Now comes the real construction..." << endl;2554 DoLog(0) && (Log() << Verbose(0) << endl << "Now comes the real construction..." << endl); 2265 2555 do { 2266 Log() << Verbose(0) << endl << endl;2267 Log() << Verbose(0) << "============Molecule list=======================" << endl;2556 DoLog(0) && (Log() << Verbose(0) << endl << endl); 2557 DoLog(0) && (Log() << Verbose(0) << "============Molecule list=======================" << endl); 2268 2558 molecules->Enumerate((ofstream *)&cout); 2269 Log() << Verbose(0) << "============Menu===============================" << endl;2270 Log() << Verbose(0) << "a - set molecule (in)active" << endl;2271 Log() << Verbose(0) << "e - edit molecules (load, parse, save)" << endl;2272 Log() << Verbose(0) << "g - globally manipulate atoms in molecule" << endl;2273 Log() << Verbose(0) << "M - Merge molecules" << endl;2274 Log() << Verbose(0) << "m - manipulate atoms" << endl;2275 Log() << Verbose(0) << "-----------------------------------------------" << endl;2276 Log() << Verbose(0) << "c - edit the current configuration" << endl;2277 Log() << Verbose(0) << "-----------------------------------------------" << endl;2278 Log() << Verbose(0) << "s - save current setup to config file" << endl;2279 Log() << Verbose(0) << "T - call the current test routine" << endl;2280 Log() << Verbose(0) << "q - quit" << endl;2281 Log() << Verbose(0) << "===============================================" << endl;2282 Log() << Verbose(0) << "Input: ";2559 DoLog(0) && (Log() << Verbose(0) << "============Menu===============================" << endl); 2560 DoLog(0) && (Log() << Verbose(0) << "a - set molecule (in)active" << endl); 2561 DoLog(0) && (Log() << Verbose(0) << "e - edit molecules (load, parse, save)" << endl); 2562 DoLog(0) && (Log() << Verbose(0) << "g - globally manipulate atoms in molecule" << endl); 2563 DoLog(0) && (Log() << Verbose(0) << "M - Merge molecules" << endl); 2564 DoLog(0) && (Log() << Verbose(0) << "m - manipulate atoms" << endl); 2565 DoLog(0) && (Log() << Verbose(0) << "-----------------------------------------------" << endl); 2566 DoLog(0) && (Log() << Verbose(0) << "c - edit the current configuration" << endl); 2567 DoLog(0) && (Log() << Verbose(0) << "-----------------------------------------------" << endl); 2568 DoLog(0) && (Log() << Verbose(0) << "s - save current setup to config file" << endl); 2569 DoLog(0) && (Log() << Verbose(0) << "T - call the current test routine" << endl); 2570 DoLog(0) && (Log() << Verbose(0) << "q - quit" << endl); 2571 DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl); 2572 DoLog(0) && (Log() << Verbose(0) << "Input: "); 2283 2573 cin >> choice; 2284 2574 … … 2286 2576 case 'a': // (in)activate molecule 2287 2577 { 2288 Log() << Verbose(0) << "Enter index of molecule: ";2578 DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: "); 2289 2579 cin >> j; 2290 2580 for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) … … 2332 2622 // save element data base 2333 2623 if (periode->StorePeriodentafel(configuration->databasepath)) //ElementsFileName 2334 Log() << Verbose(0) << "Saving of elements.db successful." << endl;2624 DoLog(0) && (Log() << Verbose(0) << "Saving of elements.db successful." << endl); 2335 2625 else 2336 Log() << Verbose(0) << "Saving of elements.db failed." << endl;2626 DoLog(0) && (Log() << Verbose(0) << "Saving of elements.db failed." << endl); 2337 2627 2338 2628 delete(molecules); // also free's all molecules contained
Note:
See TracChangeset
for help on using the changeset viewer.
