Changes in src/builder.cpp [481601:a67d19]
- File:
-
- 1 edited
-
src/builder.cpp (modified) (135 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/builder.cpp
r481601 ra67d19 52 52 #include <cstring> 53 53 54 #include "analysis_bonds.hpp" 54 55 #include "analysis_correlation.hpp" 55 56 #include "atom.hpp" … … 68 69 #include "periodentafel.hpp" 69 70 #include "version.h" 71 #include "World.hpp" 70 72 71 73 /********************************************* Subsubmenu routine ************************************/ … … 84 86 bool valid; 85 87 86 Log()<< Verbose(0) << "===========ADD ATOM============================" << endl;87 Log()<< Verbose(0) << " a - state absolute coordinates of atom" << endl;88 Log()<< Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;89 Log()<< Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;90 Log()<< Verbose(0) << " d - state two atoms, two angles and a distance" << endl;91 Log()<< Verbose(0) << " e - least square distance position to a set of atoms" << endl;92 Log()<< Verbose(0) << "all else - go back" << endl;93 Log()<< Verbose(0) << "===============================================" << endl;94 Log()<< Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;95 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: "; 96 98 cin >> choice; 97 99 98 100 switch (choice) { 99 101 default: 100 eLog() << Verbose(2) << "Not a valid choice." << endl;102 DoeLog(2) && (eLog()<< Verbose(2) << "Not a valid choice." << endl); 101 103 break; 102 104 case 'a': // absolute coordinates of atom 103 Log()<< Verbose(0) << "Enter absolute coordinates." << endl;105 cout << Verbose(0) << "Enter absolute coordinates." << endl; 104 106 first = new atom; 105 first->x.AskPosition( mol->cell_size, false);107 first->x.AskPosition(World::get()->cell_size, false); 106 108 first->type = periode->AskElement(); // give type 107 109 mol->AddAtom(first); // add to molecule … … 112 114 valid = true; 113 115 do { 114 if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;115 Log()<< Verbose(0) << "Enter reference coordinates." << endl;116 x.AskPosition( mol->cell_size, true);117 Log()<< Verbose(0) << "Enter relative coordinates." << endl;118 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); 119 121 first->x.AddVector((const Vector *)&x); 120 Log()<< Verbose(0) << "\n";122 cout << Verbose(0) << "\n"; 121 123 } while (!(valid = mol->CheckBounds((const Vector *)&first->x))); 122 124 first->type = periode->AskElement(); // give type … … 128 130 valid = true; 129 131 do { 130 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); 131 133 second = mol->AskAtom("Enter atom number: "); 132 Log() << Verbose(0) << "Enter relative coordinates." << endl;133 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); 134 136 for (int i=NDIM;i--;) { 135 137 first->x.x[i] += second->x.x[i]; … … 145 147 do { 146 148 if (!valid) { 147 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); 148 150 } 149 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; 150 152 second = mol->AskAtom("Enter central atom: "); 151 153 third = mol->AskAtom("Enter second atom (specifying the axis for first angle): "); … … 158 160 c *= M_PI/180.; 159 161 bound(&c, -M_PI, M_PI); 160 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; 161 163 /* 162 164 second->Output(1,1,(ofstream *)&cout); … … 170 172 171 173 if (!z.SolveSystem(&x,&y,&n, b, c, a)) { 172 Log() << Verbose(0) << "Failure solving self-dependent linear system!" << endl;174 coutg() << Verbose(0) << "Failure solving self-dependent linear system!" << endl; 173 175 continue; 174 176 } 175 Log() << Verbose(0) << "resulting relative coordinates: ";177 DoLog(0) && (Log() << Verbose(0) << "resulting relative coordinates: "); 176 178 z.Output(); 177 Log() << Verbose(0) << endl;179 DoLog(0) && (Log() << Verbose(0) << endl); 178 180 */ 179 181 // calc axis vector … … 183 185 Log() << Verbose(0) << "x: ", 184 186 x.Output(); 185 Log() << Verbose(0) << endl;187 DoLog(0) && (Log() << Verbose(0) << endl); 186 188 z.MakeNormalVector(&second->x,&third->x,&fourth->x); 187 189 Log() << Verbose(0) << "z: ", 188 190 z.Output(); 189 Log() << Verbose(0) << endl;191 DoLog(0) && (Log() << Verbose(0) << endl); 190 192 y.MakeNormalVector(&x,&z); 191 193 Log() << Verbose(0) << "y: ", 192 194 y.Output(); 193 Log() << Verbose(0) << endl;195 DoLog(0) && (Log() << Verbose(0) << endl); 194 196 195 197 // rotate vector around first angle … … 198 200 Log() << Verbose(0) << "Rotated vector: ", 199 201 first->x.Output(); 200 Log() << Verbose(0) << endl;202 DoLog(0) && (Log() << Verbose(0) << endl); 201 203 // remove the projection onto the rotation plane of the second angle 202 204 n.CopyVector(&y); … … 204 206 Log() << Verbose(0) << "N1: ", 205 207 n.Output(); 206 Log() << Verbose(0) << endl;208 DoLog(0) && (Log() << Verbose(0) << endl); 207 209 first->x.SubtractVector(&n); 208 210 Log() << Verbose(0) << "Subtracted vector: ", 209 211 first->x.Output(); 210 Log() << Verbose(0) << endl;212 DoLog(0) && (Log() << Verbose(0) << endl); 211 213 n.CopyVector(&z); 212 214 n.Scale(first->x.ScalarProduct(&z)); 213 215 Log() << Verbose(0) << "N2: ", 214 216 n.Output(); 215 Log() << Verbose(0) << endl;217 DoLog(0) && (Log() << Verbose(0) << endl); 216 218 first->x.SubtractVector(&n); 217 219 Log() << Verbose(0) << "2nd subtracted vector: ", 218 220 first->x.Output(); 219 Log() << Verbose(0) << endl;221 DoLog(0) && (Log() << Verbose(0) << endl); 220 222 221 223 // rotate another vector around second angle … … 224 226 Log() << Verbose(0) << "2nd Rotated vector: ", 225 227 n.Output(); 226 Log() << Verbose(0) << endl;228 DoLog(0) && (Log() << Verbose(0) << endl); 227 229 228 230 // add the two linear independent vectors … … 232 234 first->x.AddVector(&second->x); 233 235 234 Log() << Verbose(0) << "resulting coordinates: ";236 DoLog(0) && (Log() << Verbose(0) << "resulting coordinates: "); 235 237 first->x.Output(); 236 Log() << Verbose(0) << endl;238 DoLog(0) && (Log() << Verbose(0) << endl); 237 239 } while (!(valid = mol->CheckBounds((const Vector *)&first->x))); 238 240 first->type = periode->AskElement(); // give type … … 247 249 atoms[i] = NULL; 248 250 int i=0, j=0; 249 Log()<< Verbose(0) << "Now we need at least three molecules.\n";251 cout << Verbose(0) << "Now we need at least three molecules.\n"; 250 252 do { 251 Log()<< Verbose(0) << "Enter " << i+1 << "th atom: ";253 cout << Verbose(0) << "Enter " << i+1 << "th atom: "; 252 254 cin >> j; 253 255 if (j != -1) { … … 264 266 } else { 265 267 delete first; 266 Log()<< Verbose(0) << "Please enter at least two vectors!\n";268 cout << Verbose(0) << "Please enter at least two vectors!\n"; 267 269 } 268 270 break; … … 278 280 char choice; // menu choice char 279 281 280 Log()<< Verbose(0) << "===========CENTER ATOMS=========================" << endl;281 Log()<< Verbose(0) << " a - on origin" << endl;282 Log()<< Verbose(0) << " b - on center of gravity" << endl;283 Log()<< Verbose(0) << " c - within box with additional boundary" << endl;284 Log()<< Verbose(0) << " d - within given simulation box" << endl;285 Log()<< Verbose(0) << "all else - go back" << endl;286 Log()<< Verbose(0) << "===============================================" << endl;287 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: "; 288 290 cin >> choice; 289 291 290 292 switch (choice) { 291 293 default: 292 Log()<< Verbose(0) << "Not a valid choice." << endl;294 cout << Verbose(0) << "Not a valid choice." << endl; 293 295 break; 294 296 case 'a': 295 Log()<< Verbose(0) << "Centering atoms in config file on origin." << endl;297 cout << Verbose(0) << "Centering atoms in config file on origin." << endl; 296 298 mol->CenterOrigin(); 297 299 break; 298 300 case 'b': 299 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; 300 302 mol->CenterPeriodic(); 301 303 break; 302 304 case 'c': 303 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; 304 306 for (int i=0;i<NDIM;i++) { 305 Log()<< Verbose(0) << "Enter axis " << i << " boundary: ";307 cout << Verbose(0) << "Enter axis " << i << " boundary: "; 306 308 cin >> y.x[i]; 307 309 } … … 314 316 break; 315 317 case 'd': 316 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; 317 319 for (int i=0;i<NDIM;i++) { 318 Log()<< Verbose(0) << "Enter axis " << i << " boundary: ";320 cout << Verbose(0) << "Enter axis " << i << " boundary: "; 319 321 cin >> x.x[i]; 320 322 } … … 337 339 char choice; // menu choice char 338 340 339 Log()<< Verbose(0) << "===========ALIGN ATOMS=========================" << endl;340 Log()<< Verbose(0) << " a - state three atoms defining align plane" << endl;341 Log()<< Verbose(0) << " b - state alignment vector" << endl;342 Log()<< Verbose(0) << " c - state two atoms in alignment direction" << endl;343 Log()<< Verbose(0) << " d - align automatically by least square fit" << endl;344 Log()<< Verbose(0) << "all else - go back" << endl;345 Log()<< Verbose(0) << "===============================================" << endl;346 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: "; 347 349 cin >> choice; 348 350 … … 357 359 break; 358 360 case 'b': // normal vector of mirror plane 359 Log()<< Verbose(0) << "Enter normal vector of mirror plane." << endl;360 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); 361 363 n.Normalize(); 362 364 break; … … 377 379 fscanf(stdin, "%3s", shorthand); 378 380 } while ((param.type = periode->FindElement(shorthand)) == NULL); 379 Log()<< Verbose(0) << "Element is " << param.type->name << endl;381 cout << Verbose(0) << "Element is " << param.type->name << endl; 380 382 mol->GetAlignvector(¶m); 381 383 for (int i=NDIM;i--;) { … … 384 386 } 385 387 gsl_vector_free(param.x); 386 Log()<< Verbose(0) << "Offset vector: ";388 cout << Verbose(0) << "Offset vector: "; 387 389 x.Output(); 388 Log() << Verbose(0) << endl;390 DoLog(0) && (Log() << Verbose(0) << endl); 389 391 n.Normalize(); 390 392 break; 391 393 }; 392 Log() << Verbose(0) << "Alignment vector: ";394 DoLog(0) && (Log() << Verbose(0) << "Alignment vector: "); 393 395 n.Output(); 394 Log() << Verbose(0) << endl;396 DoLog(0) && (Log() << Verbose(0) << endl); 395 397 mol->Align(&n); 396 398 }; … … 405 407 char choice; // menu choice char 406 408 407 Log() << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;408 Log() << Verbose(0) << " a - state three atoms defining mirror plane" << endl;409 Log() << Verbose(0) << " b - state normal vector of mirror plane" << endl;410 Log() << Verbose(0) << " c - state two atoms in normal direction" << endl;411 Log() << Verbose(0) << "all else - go back" << endl;412 Log() << Verbose(0) << "===============================================" << endl;413 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: "); 414 416 cin >> choice; 415 417 … … 424 426 break; 425 427 case 'b': // normal vector of mirror plane 426 Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;427 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); 428 430 n.Normalize(); 429 431 break; … … 437 439 break; 438 440 }; 439 Log() << Verbose(0) << "Normal vector: ";441 DoLog(0) && (Log() << Verbose(0) << "Normal vector: "); 440 442 n.Output(); 441 Log() << Verbose(0) << endl;443 DoLog(0) && (Log() << Verbose(0) << endl); 442 444 mol->Mirror((const Vector *)&n); 443 445 }; … … 453 455 char choice; // menu choice char 454 456 455 Log() << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;456 Log() << Verbose(0) << " a - state atom for removal by number" << endl;457 Log() << Verbose(0) << " b - keep only in radius around atom" << endl;458 Log() << Verbose(0) << " c - remove this with one axis greater value" << endl;459 Log() << Verbose(0) << "all else - go back" << endl;460 Log() << Verbose(0) << "===============================================" << endl;461 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: "); 462 464 cin >> choice; 463 465 … … 466 468 case 'a': 467 469 if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: "))) 468 Log() << Verbose(1) << "Atom removed." << endl;470 DoLog(1) && (Log() << Verbose(1) << "Atom removed." << endl); 469 471 else 470 Log() << Verbose(1) << "Atom not found." << endl;472 DoLog(1) && (Log() << Verbose(1) << "Atom not found." << endl); 471 473 break; 472 474 case 'b': 473 475 second = mol->AskAtom("Enter number of atom as reference point: "); 474 Log() << Verbose(0) << "Enter radius: ";476 DoLog(0) && (Log() << Verbose(0) << "Enter radius: "); 475 477 cin >> tmp1; 476 478 first = mol->start; … … 484 486 break; 485 487 case 'c': 486 Log() << Verbose(0) << "Which axis is it: ";488 DoLog(0) && (Log() << Verbose(0) << "Which axis is it: "); 487 489 cin >> axis; 488 Log() << Verbose(0) << "Lower boundary: ";490 DoLog(0) && (Log() << Verbose(0) << "Lower boundary: "); 489 491 cin >> tmp1; 490 Log() << Verbose(0) << "Upper boundary: ";492 DoLog(0) && (Log() << Verbose(0) << "Upper boundary: "); 491 493 cin >> tmp2; 492 494 first = mol->start; … … 518 520 char choice; // menu choice char 519 521 520 Log() << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;521 Log() << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;522 Log() << Verbose(0) << " b - calculate bond length between two atoms" << endl;523 Log() << Verbose(0) << " c - calculate bond angle" << endl;524 Log() << Verbose(0) << " d - calculate principal axis of the system" << endl;525 Log() << Verbose(0) << " e - calculate volume of the convex envelope" << endl;526 Log() << Verbose(0) << " f - calculate temperature from current velocity" << endl;527 Log() << Verbose(0) << " g - output all temperatures per step from velocities" << endl;528 Log() << Verbose(0) << "all else - go back" << endl;529 Log() << Verbose(0) << "===============================================" << endl;530 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: "); 531 533 cin >> choice; 532 534 533 535 switch(choice) { 534 536 default: 535 Log() << Verbose(1) << "Not a valid choice." << endl;537 DoLog(1) && (Log() << Verbose(1) << "Not a valid choice." << endl); 536 538 break; 537 539 case 'a': … … 565 567 x.SubtractVector((const Vector *)&second->x); 566 568 tmp1 = x.Norm(); 567 Log() << Verbose(1) << "Distance vector is ";569 DoLog(1) && (Log() << Verbose(1) << "Distance vector is "); 568 570 x.Output(); 569 Log() << Verbose(0) << "." << endl << "Norm of distance is " << tmp1 << "." << endl;571 DoLog(0) && (Log() << Verbose(0) << "." << endl << "Norm of distance is " << tmp1 << "." << endl); 570 572 break; 571 573 572 574 case 'c': 573 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); 574 576 first = mol->AskAtom("Enter first atom: "); 575 577 second = mol->AskAtom("Enter central atom: "); … … 580 582 y.CopyVector((const Vector *)&third->x); 581 583 y.SubtractVector((const Vector *)&second->x); 582 Log() << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";583 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); 584 586 break; 585 587 case 'd': 586 Log() << Verbose(0) << "Evaluating prinicipal axis." << endl;587 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]: "); 588 590 cin >> Z; 589 591 if ((Z >=0) && (Z <=1)) … … 594 596 case 'e': 595 597 { 596 Log() << Verbose(0) << "Evaluating volume of the convex envelope.";598 DoLog(0) && (Log() << Verbose(0) << "Evaluating volume of the convex envelope."); 597 599 class Tesselation *TesselStruct = NULL; 598 600 const LinkedCell *LCList = NULL; … … 600 602 FindConvexBorder(mol, TesselStruct, LCList, NULL); 601 603 double clustervolume = VolumeOfConvexEnvelope(TesselStruct, configuration); 602 Log() << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl;\604 DoLog(0) && (Log() << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl);\ 603 605 delete(LCList); 604 606 delete(TesselStruct); … … 611 613 { 612 614 char filename[255]; 613 Log() << Verbose(0) << "Please enter filename: " << endl;615 DoLog(0) && (Log() << Verbose(0) << "Please enter filename: " << endl); 614 616 cin >> filename; 615 Log() << Verbose(1) << "Storing temperatures in " << filename << "." << endl;617 DoLog(1) && (Log() << Verbose(1) << "Storing temperatures in " << filename << "." << endl); 616 618 ofstream *output = new ofstream(filename, ios::trunc); 617 619 if (!mol->OutputTemperatureFromTrajectories(output, 0, mol->MDSteps)) 618 Log() << Verbose(2) << "File could not be written." << endl;620 DoLog(2) && (Log() << Verbose(2) << "File could not be written." << endl); 619 621 else 620 Log() << Verbose(2) << "File stored." << endl;622 DoLog(2) && (Log() << Verbose(2) << "File stored." << endl); 621 623 output->close(); 622 624 delete(output); … … 635 637 clock_t start, end; 636 638 637 Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;638 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: "); 639 641 cin >> Order1; 640 642 if (mol->first->next != mol->last) { // there are bonds … … 642 644 mol->FragmentMolecule(Order1, configuration); 643 645 end = clock(); 644 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); 645 647 } else 646 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); 647 649 }; 648 650 … … 655 657 static void ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration) 656 658 { 657 atom *first, *second ;659 atom *first, *second, *third; 658 660 molecule *mol = NULL; 659 661 Vector x,y,z,n; // coordinates for absolute point in cell volume … … 663 665 bool valid; 664 666 665 Log() << Verbose(0) << "=========MANIPULATE ATOMS======================" << endl; 666 Log() << Verbose(0) << "a - add an atom" << endl; 667 Log() << Verbose(0) << "r - remove an atom" << endl; 668 Log() << Verbose(0) << "b - scale a bond between atoms" << endl; 669 Log() << Verbose(0) << "u - change an atoms element" << endl; 670 Log() << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl; 671 Log() << Verbose(0) << "all else - go back" << endl; 672 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); 673 676 if (molecules->NumberOfActiveMolecules() > 1) 674 eLog() << Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl;675 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: "); 676 679 cin >> choice; 677 680 678 681 switch (choice) { 679 682 default: 680 Log() << Verbose(0) << "Not a valid choice." << endl;683 DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl); 681 684 break; 682 685 … … 685 688 if ((*ListRunner)->ActiveFlag) { 686 689 mol = *ListRunner; 687 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); 688 691 AddAtoms(periode, mol); 689 692 } … … 694 697 if ((*ListRunner)->ActiveFlag) { 695 698 mol = *ListRunner; 696 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;697 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); 698 701 first = mol->AskAtom("Enter first (fixed) atom: "); 699 702 second = mol->AskAtom("Enter second (shifting) atom: "); … … 702 705 minBond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]); 703 706 minBond = sqrt(minBond); 704 Log() << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << minBond << " a.u." << endl;705 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.]: "); 706 709 cin >> bond; 707 710 for (int i=NDIM;i--;) { … … 717 720 if ((*ListRunner)->ActiveFlag) { 718 721 mol = *ListRunner; 719 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;720 Log() << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;721 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: "); 722 725 factor = new double[NDIM]; 723 726 cin >> factor[0]; … … 734 737 if ((*ListRunner)->ActiveFlag) { 735 738 mol = *ListRunner; 736 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); 737 740 MeasureAtoms(periode, mol, configuration); 738 741 } … … 743 746 if ((*ListRunner)->ActiveFlag) { 744 747 mol = *ListRunner; 745 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); 746 749 RemoveAtoms(mol); 747 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 } 748 784 break; 749 785 … … 753 789 int Z; 754 790 mol = *ListRunner; 755 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); 756 792 first = NULL; 757 793 do { 758 Log() << Verbose(0) << "Change the element of which atom: ";794 DoLog(0) && (Log() << Verbose(0) << "Change the element of which atom: "); 759 795 cin >> Z; 760 796 } while ((first = mol->FindAtom(Z)) == NULL); 761 Log() << Verbose(0) << "New element by atomic number Z: ";797 DoLog(0) && (Log() << Verbose(0) << "New element by atomic number Z: "); 762 798 cin >> Z; 763 799 first->type = periode->FindElement(Z); 764 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); 765 801 } 766 802 break; … … 783 819 MoleculeLeafClass *Subgraphs = NULL; 784 820 785 Log() << Verbose(0) << "=========MANIPULATE GLOBALLY===================" << endl;786 Log() << Verbose(0) << "c - scale by unit transformation" << endl;787 Log() << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;788 Log() << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;789 Log() << Verbose(0) << "g - center atoms in box" << endl;790 Log() << Verbose(0) << "i - realign molecule" << endl;791 Log() << Verbose(0) << "m - mirror all molecules" << endl;792 Log() << Verbose(0) << "o - create connection matrix" << endl;793 Log() << Verbose(0) << "t - translate molecule by vector" << endl;794 Log() << Verbose(0) << "all else - go back" << endl;795 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); 796 832 if (molecules->NumberOfActiveMolecules() > 1) 797 eLog() << Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl;798 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: "); 799 835 cin >> choice; 800 836 801 837 switch (choice) { 802 838 default: 803 Log() << Verbose(0) << "Not a valid choice." << endl;839 DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl); 804 840 break; 805 841 … … 808 844 if ((*ListRunner)->ActiveFlag) { 809 845 mol = *ListRunner; 810 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;811 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]: "); 812 848 cin >> axis; 813 Log() << Verbose(0) << "State the factor: ";849 DoLog(0) && (Log() << Verbose(0) << "State the factor: "); 814 850 cin >> faktor; 815 851 … … 828 864 } 829 865 if (count != j) 830 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); 831 867 x.Zero(); 832 868 y.Zero(); 833 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 834 870 for (int i=1;i<faktor;i++) { // then add this list with respective translation factor times 835 871 x.AddVector(&y); // per factor one cell width further … … 854 890 mol->Translate(&x); 855 891 } 856 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; 857 893 } 858 894 } … … 867 903 if ((*ListRunner)->ActiveFlag) { 868 904 mol = *ListRunner; 869 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); 870 906 CenterAtoms(mol); 871 907 } … … 876 912 if ((*ListRunner)->ActiveFlag) { 877 913 mol = *ListRunner; 878 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); 879 915 AlignAtoms(periode, mol); 880 916 } … … 885 921 if ((*ListRunner)->ActiveFlag) { 886 922 mol = *ListRunner; 887 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); 888 924 MirrorAtoms(mol); 889 925 } … … 896 932 double bonddistance; 897 933 clock_t start,end; 898 Log() << Verbose(0) << "What's the maximum bond distance: ";934 DoLog(0) && (Log() << Verbose(0) << "What's the maximum bond distance: "); 899 935 cin >> bonddistance; 900 936 start = clock(); 901 937 mol->CreateAdjacencyList(bonddistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL); 902 938 end = clock(); 903 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); 904 940 } 905 941 break; … … 909 945 if ((*ListRunner)->ActiveFlag) { 910 946 mol = *ListRunner; 911 Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;912 Log() << Verbose(0) << "Enter translation vector." << endl;913 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); 914 950 mol->Center.AddVector((const Vector *)&x); 915 951 } … … 938 974 molecule *mol = NULL; 939 975 940 Log() << Verbose(0) << "==========EDIT MOLECULES=====================" << endl;941 Log() << Verbose(0) << "c - create new molecule" << endl;942 Log() << Verbose(0) << "l - load molecule from xyz file" << endl;943 Log() << Verbose(0) << "n - change molecule's name" << endl;944 Log() << Verbose(0) << "N - give molecules filename" << endl;945 Log() << Verbose(0) << "p - parse atoms in xyz file into molecule" << endl;946 Log() << Verbose(0) << "r - remove a molecule" << endl;947 Log() << Verbose(0) << "all else - go back" << endl;948 Log() << Verbose(0) << "===============================================" << endl;949 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: "); 950 986 cin >> choice; 951 987 952 988 switch (choice) { 953 989 default: 954 Log() << Verbose(0) << "Not a valid choice." << endl;990 DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl); 955 991 break; 956 992 case 'c': … … 962 998 { 963 999 char filename[MAXSTRINGSIZE]; 964 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); 965 1001 mol = new molecule(periode); 966 1002 do { 967 Log() << Verbose(0) << "Enter file name: ";1003 DoLog(0) && (Log() << Verbose(0) << "Enter file name: "); 968 1004 cin >> filename; 969 1005 } while (!mol->AddXYZFile(filename)); … … 971 1007 // center at set box dimensions 972 1008 mol->CenterEdge(¢er); 973 mol->cell_size[0] = center.x[0]; 974 mol->cell_size[1] = 0; 975 mol->cell_size[2] = center.x[1]; 976 mol->cell_size[3] = 0; 977 mol->cell_size[4] = 0; 978 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]; 979 1016 molecules->insert(mol); 980 1017 } … … 985 1022 char filename[MAXSTRINGSIZE]; 986 1023 do { 987 Log() << Verbose(0) << "Enter index of molecule: ";1024 DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: "); 988 1025 cin >> nr; 989 1026 mol = molecules->ReturnIndex(nr); 990 1027 } while (mol == NULL); 991 Log() << Verbose(0) << "Enter name: ";1028 DoLog(0) && (Log() << Verbose(0) << "Enter name: "); 992 1029 cin >> filename; 993 1030 strcpy(mol->name, filename); … … 999 1036 char filename[MAXSTRINGSIZE]; 1000 1037 do { 1001 Log() << Verbose(0) << "Enter index of molecule: ";1038 DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: "); 1002 1039 cin >> nr; 1003 1040 mol = molecules->ReturnIndex(nr); 1004 1041 } while (mol == NULL); 1005 Log() << Verbose(0) << "Enter name: ";1042 DoLog(0) && (Log() << Verbose(0) << "Enter name: "); 1006 1043 cin >> filename; 1007 1044 mol->SetNameFromFilename(filename); … … 1014 1051 mol = NULL; 1015 1052 do { 1016 Log() << Verbose(0) << "Enter index of molecule: ";1053 DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: "); 1017 1054 cin >> nr; 1018 1055 mol = molecules->ReturnIndex(nr); 1019 1056 } while (mol == NULL); 1020 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); 1021 1058 do { 1022 Log() << Verbose(0) << "Enter file name: ";1059 DoLog(0) && (Log() << Verbose(0) << "Enter file name: "); 1023 1060 cin >> filename; 1024 1061 } while (!mol->AddXYZFile(filename)); … … 1028 1065 1029 1066 case 'r': 1030 Log() << Verbose(0) << "Enter index of molecule: ";1067 DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: "); 1031 1068 cin >> nr; 1032 1069 count = 1; … … 1051 1088 char choice; // menu choice char 1052 1089 1053 Log() << Verbose(0) << "===========MERGE MOLECULES=====================" << endl; 1054 Log() << Verbose(0) << "a - simple add of one molecule to another" << endl; 1055 Log() << Verbose(0) << "e - embedding merge of two molecules" << endl; 1056 Log() << Verbose(0) << "m - multi-merge of all molecules" << endl; 1057 Log() << Verbose(0) << "s - scatter merge of two molecules" << endl; 1058 Log() << Verbose(0) << "t - simple merge of two molecules" << endl; 1059 Log() << Verbose(0) << "all else - go back" << endl; 1060 Log() << Verbose(0) << "===============================================" << endl; 1061 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: "); 1062 1103 cin >> choice; 1063 1104 1064 1105 switch (choice) { 1065 1106 default: 1066 Log() << Verbose(0) << "Not a valid choice." << endl;1107 DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl); 1067 1108 break; 1068 1109 … … 1073 1114 { 1074 1115 do { 1075 Log() << Verbose(0) << "Enter index of destination molecule: ";1116 DoLog(0) && (Log() << Verbose(0) << "Enter index of destination molecule: "); 1076 1117 cin >> dest; 1077 1118 destmol = molecules->ReturnIndex(dest); 1078 1119 } while ((destmol == NULL) && (dest != -1)); 1079 1120 do { 1080 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: "); 1081 1122 cin >> src; 1082 1123 srcmol = molecules->ReturnIndex(src); … … 1088 1129 break; 1089 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 1090 1183 case 'e': 1091 1184 { … … 1093 1186 molecule *srcmol = NULL, *destmol = NULL; 1094 1187 do { 1095 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): "); 1096 1189 cin >> src; 1097 1190 srcmol = molecules->ReturnIndex(src); 1098 1191 } while ((srcmol == NULL) && (src != -1)); 1099 1192 do { 1100 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): "); 1101 1194 cin >> dest; 1102 1195 destmol = molecules->ReturnIndex(dest); … … 1107 1200 break; 1108 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 1109 1212 case 'm': 1110 1213 { … … 1112 1215 molecule *mol = NULL; 1113 1216 do { 1114 Log() << Verbose(0) << "Enter index of molecule to merge into: ";1217 DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule to merge into: "); 1115 1218 cin >> nr; 1116 1219 mol = molecules->ReturnIndex(nr); … … 1129 1232 1130 1233 case 's': 1131 Log() << Verbose(0) << "Not implemented yet." << endl;1234 DoLog(0) && (Log() << Verbose(0) << "Not implemented yet." << endl); 1132 1235 break; 1133 1236 … … 1138 1241 { 1139 1242 do { 1140 Log() << Verbose(0) << "Enter index of destination molecule: ";1243 DoLog(0) && (Log() << Verbose(0) << "Enter index of destination molecule: "); 1141 1244 cin >> dest; 1142 1245 destmol = molecules->ReturnIndex(dest); 1143 1246 } while ((destmol == NULL) && (dest != -1)); 1144 1247 do { 1145 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: "); 1146 1249 cin >> src; 1147 1250 srcmol = molecules->ReturnIndex(src); … … 1172 1275 mol = (molecules->ListOfMolecules.front())->CopyMolecule(); 1173 1276 else { 1174 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 ... "); 1175 1278 performCriticalExit(); 1176 1279 return; … … 1179 1282 1180 1283 // generate some KeySets 1181 Log() << Verbose(0) << "Generating KeySets." << endl;1284 DoLog(0) && (Log() << Verbose(0) << "Generating KeySets." << endl); 1182 1285 KeySet TestSets[mol->AtomCount+1]; 1183 1286 i=1; … … 1189 1292 i++; 1190 1293 } 1191 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); 1192 1295 KeySetTestPair test; 1193 1296 test = TestSets[mol->AtomCount-1].insert(Walker->nr); 1194 1297 if (test.second) { 1195 Log() << Verbose(1) << "Insertion worked?!" << endl;1298 DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl); 1196 1299 } else { 1197 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); 1198 1301 } 1199 1302 TestSets[mol->AtomCount].insert(mol->end->previous->nr); … … 1201 1304 1202 1305 // constructing Graph structure 1203 Log() << Verbose(0) << "Generating Subgraph class." << endl;1306 DoLog(0) && (Log() << Verbose(0) << "Generating Subgraph class." << endl); 1204 1307 Graph Subgraphs; 1205 1308 1206 1309 // insert KeySets into Subgraphs 1207 Log() << Verbose(0) << "Inserting KeySets into Subgraph class." << endl;1310 DoLog(0) && (Log() << Verbose(0) << "Inserting KeySets into Subgraph class." << endl); 1208 1311 for (int j=0;j<mol->AtomCount;j++) { 1209 1312 Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.))); 1210 1313 } 1211 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); 1212 1315 GraphTestPair test2; 1213 1316 test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.))); 1214 1317 if (test2.second) { 1215 Log() << Verbose(1) << "Insertion worked?!" << endl;1318 DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl); 1216 1319 } else { 1217 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); 1218 1321 } 1219 1322 1220 1323 // show graphs 1221 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); 1222 1325 Graph::iterator A = Subgraphs.begin(); 1223 1326 while (A != Subgraphs.end()) { 1224 Log() << Verbose(0) << (*A).second.first << ": ";1327 DoLog(0) && (Log() << Verbose(0) << (*A).second.first << ": "); 1225 1328 KeySet::iterator key = (*A).first.begin(); 1226 1329 comp = -1; 1227 1330 while (key != (*A).first.end()) { 1228 1331 if ((*key) > comp) 1229 Log() << Verbose(0) << (*key) << " ";1332 DoLog(0) && (Log() << Verbose(0) << (*key) << " "); 1230 1333 else 1231 Log() << Verbose(0) << (*key) << "! ";1334 DoLog(0) && (Log() << Verbose(0) << (*key) << "! "); 1232 1335 comp = (*key); 1233 1336 key++; 1234 1337 } 1235 Log() << Verbose(0) << endl;1338 DoLog(0) && (Log() << Verbose(0) << endl); 1236 1339 A++; 1237 1340 } … … 1253 1356 1254 1357 if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) { 1255 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); 1256 1359 } 1257 1360 … … 1262 1365 if (output == NULL) 1263 1366 strcpy(filename,"main_pcp_linux"); 1264 Log() << Verbose(0) << "Saving as pdb input ";1367 DoLog(0) && (Log() << Verbose(0) << "Saving as pdb input "); 1265 1368 if (configuration->SavePDB(filename, molecules)) 1266 Log() << Verbose(0) << "done." << endl;1369 DoLog(0) && (Log() << Verbose(0) << "done." << endl); 1267 1370 else 1268 Log() << Verbose(0) << "failed." << endl;1371 DoLog(0) && (Log() << Verbose(0) << "failed." << endl); 1269 1372 1270 1373 // then save as tremolo data file … … 1273 1376 if (output == NULL) 1274 1377 strcpy(filename,"main_pcp_linux"); 1275 Log() << Verbose(0) << "Saving as tremolo data input ";1378 DoLog(0) && (Log() << Verbose(0) << "Saving as tremolo data input "); 1276 1379 if (configuration->SaveTREMOLO(filename, molecules)) 1277 Log() << Verbose(0) << "done." << endl;1380 DoLog(0) && (Log() << Verbose(0) << "done." << endl); 1278 1381 else 1279 Log() << Verbose(0) << "failed." << endl;1382 DoLog(0) && (Log() << Verbose(0) << "failed." << endl); 1280 1383 1281 1384 // translate each to its center and merge all molecules in MoleculeListClass into this molecule … … 1297 1400 } 1298 1401 1299 Log() << Verbose(0) << "Storing configuration ... " << endl;1402 DoLog(0) && (Log() << Verbose(0) << "Storing configuration ... " << endl); 1300 1403 // get correct valence orbitals 1301 1404 mol->CalculateOrbitals(*configuration); … … 1313 1416 output.close(); 1314 1417 output.clear(); 1315 Log() << Verbose(0) << "Saving of config file ";1418 DoLog(0) && (Log() << Verbose(0) << "Saving of config file "); 1316 1419 if (configuration->Save(filename, periode, mol)) 1317 Log() << Verbose(0) << "successful." << endl;1420 DoLog(0) && (Log() << Verbose(0) << "successful." << endl); 1318 1421 else 1319 Log() << Verbose(0) << "failed." << endl;1422 DoLog(0) && (Log() << Verbose(0) << "failed." << endl); 1320 1423 1321 1424 // and save to xyz file … … 1330 1433 output.open(filename, ios::trunc); 1331 1434 } 1332 Log() << Verbose(0) << "Saving of XYZ file ";1435 DoLog(0) && (Log() << Verbose(0) << "Saving of XYZ file "); 1333 1436 if (mol->MDSteps <= 1) { 1334 1437 if (mol->OutputXYZ(&output)) 1335 Log() << Verbose(0) << "successful." << endl;1438 DoLog(0) && (Log() << Verbose(0) << "successful." << endl); 1336 1439 else 1337 Log() << Verbose(0) << "failed." << endl;1440 DoLog(0) && (Log() << Verbose(0) << "failed." << endl); 1338 1441 } else { 1339 1442 if (mol->OutputTrajectoriesXYZ(&output)) 1340 Log() << Verbose(0) << "successful." << endl;1443 DoLog(0) && (Log() << Verbose(0) << "successful." << endl); 1341 1444 else 1342 Log() << Verbose(0) << "failed." << endl;1445 DoLog(0) && (Log() << Verbose(0) << "failed." << endl); 1343 1446 } 1344 1447 output.close(); … … 1350 1453 if (output == NULL) 1351 1454 strcpy(filename,"main_pcp_linux"); 1352 Log() << Verbose(0) << "Saving as mpqc input ";1455 DoLog(0) && (Log() << Verbose(0) << "Saving as mpqc input "); 1353 1456 if (configuration->SaveMPQC(filename, mol)) 1354 Log() << Verbose(0) << "done." << endl;1457 DoLog(0) && (Log() << Verbose(0) << "done." << endl); 1355 1458 else 1356 Log() << Verbose(0) << "failed." << endl;1459 DoLog(0) && (Log() << Verbose(0) << "failed." << endl); 1357 1460 1358 1461 if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) { 1359 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); 1360 1463 } 1361 1464 … … 1387 1490 enum ConfigStatus configPresent = absent; 1388 1491 clock_t start,end; 1492 double MaxDistance = -1; 1389 1493 int argptr; 1390 1494 molecule *mol = NULL; … … 1398 1502 do { 1399 1503 if (argv[argptr][0] == '-') { 1400 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"); 1401 1505 argptr++; 1402 1506 switch(argv[argptr-1][1]) { … … 1404 1508 case 'H': 1405 1509 case '?': 1406 Log() << Verbose(0) << "MoleCuilder suite" << endl << "==================" << endl << endl; 1407 Log() << Verbose(0) << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl; 1408 Log() << Verbose(0) << "or simply " << argv[0] << " without arguments for interactive session." << endl; 1409 Log() << Verbose(0) << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl; 1410 Log() << Verbose(0) << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl; 1411 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; 1412 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; 1413 Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl; 1414 Log() << Verbose(0) << "\t-C <Z> <output> <bin output>\tPair Correlation analysis." << endl; 1415 Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl; 1416 Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl; 1417 Log() << Verbose(0) << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl; 1418 Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl; 1419 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; 1420 Log() << Verbose(0) << "\t-F <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl; 1421 Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl; 1422 Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl; 1423 Log() << Verbose(0) << "\t-I\t Dissect current system of molecules into a set of disconnected (subgraphs of) molecules." << endl; 1424 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; 1425 Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl; 1426 Log() << Verbose(0) << "\t-M <basis>\tSetting basis to store to MPQC config files." << endl; 1427 Log() << Verbose(0) << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl; 1428 Log() << Verbose(0) << "\t-N <radius> <file>\tGet non-convex-envelope." << endl; 1429 Log() << Verbose(0) << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl; 1430 Log() << Verbose(0) << "\t-O\tCenter atoms in origin." << endl; 1431 Log() << Verbose(0) << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl; 1432 Log() << Verbose(0) << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl; 1433 Log() << Verbose(0) << "\t-r <id>\t\tRemove an atom with given id." << endl; 1434 Log() << Verbose(0) << "\t-R <id> <radius>\t\tRemove all atoms out of sphere around a given one." << endl; 1435 Log() << Verbose(0) << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl; 1436 Log() << Verbose(0) << "\t-S <file> Store temperatures from the config file in <file>." << endl; 1437 Log() << Verbose(0) << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl; 1438 Log() << Verbose(0) << "\t-T x1 x2 x3\tTranslate periodically all atoms by this vector (x1,x2,x3)." << endl; 1439 Log() << Verbose(0) << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl; 1440 Log() << Verbose(0) << "\t-v\t\tsets verbosity (more is more)." << endl; 1441 Log() << Verbose(0) << "\t-V\t\tGives version information." << endl; 1442 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); 1443 1551 return (1); 1444 1552 break; … … 1448 1556 } 1449 1557 setVerbosity(verbosity); 1450 Log() << Verbose(0) << "Setting verbosity to " << verbosity << "." << endl;1558 DoLog(0) && (Log() << Verbose(0) << "Setting verbosity to " << verbosity << "." << endl); 1451 1559 break; 1452 1560 case 'V': 1453 Log() << Verbose(0) << argv[0] << " " << VERSIONSTRING << endl;1454 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); 1455 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 } 1456 1581 break; 1457 1582 case 'e': 1458 1583 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1459 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); 1460 1585 performCriticalExit(); 1461 1586 } else { 1462 Log() << Verbose(0) << "Using " << argv[argptr] << " as elements database." << endl;1587 DoLog(0) && (Log() << Verbose(0) << "Using " << argv[argptr] << " as elements database." << endl); 1463 1588 strncpy (configuration.databasepath, argv[argptr], MAXSTRINGSIZE-1); 1464 1589 argptr+=1; … … 1467 1592 case 'g': 1468 1593 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1469 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); 1470 1595 performCriticalExit(); 1471 1596 } else { 1472 1597 BondGraphFileName = argv[argptr]; 1473 Log() << Verbose(0) << "Using " << BondGraphFileName << " as bond length table." << endl;1598 DoLog(0) && (Log() << Verbose(0) << "Using " << BondGraphFileName << " as bond length table." << endl); 1474 1599 argptr+=1; 1475 1600 } 1476 1601 break; 1477 1602 case 'n': 1478 Log() << Verbose(0) << "I won't parse trajectories." << endl;1603 DoLog(0) && (Log() << Verbose(0) << "I won't parse trajectories." << endl); 1479 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 } 1480 1615 break; 1481 1616 default: // no match? Step on … … 1489 1624 // 3a. Parse the element database 1490 1625 if (periode->LoadPeriodentafel(configuration.databasepath)) { 1491 Log() << Verbose(0) << "Element list loaded successfully." << endl;1626 DoLog(0) && (Log() << Verbose(0) << "Element list loaded successfully." << endl); 1492 1627 //periode->Output(); 1493 1628 } else { 1494 Log() << Verbose(0) << "Element list loading failed." << endl;1629 DoLog(0) && (Log() << Verbose(0) << "Element list loading failed." << endl); 1495 1630 return 1; 1496 1631 } … … 1498 1633 if (argv[1][0] != '-') { 1499 1634 // simply create a new molecule, wherein the config file is loaded and the manipulation takes place 1500 Log() << Verbose(0) << "Config file given." << endl;1635 DoLog(0) && (Log() << Verbose(0) << "Config file given." << endl); 1501 1636 test.open(argv[1], ios::in); 1502 1637 if (test == NULL) { … … 1504 1639 output.open(argv[1], ios::out); 1505 1640 if (output == NULL) { 1506 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); 1507 1642 configPresent = absent; 1508 1643 } else { 1509 Log() << Verbose(0) << "Empty configuration file." << endl;1644 DoLog(0) && (Log() << Verbose(0) << "Empty configuration file." << endl); 1510 1645 ConfigFileName = argv[1]; 1511 1646 configPresent = empty; … … 1515 1650 test.close(); 1516 1651 ConfigFileName = argv[1]; 1517 Log() << Verbose(1) << "Specified config file found, parsing ... ";1652 DoLog(1) && (Log() << Verbose(1) << "Specified config file found, parsing ... "); 1518 1653 switch (configuration.TestSyntax(ConfigFileName, periode)) { 1519 1654 case 1: 1520 Log() << Verbose(0) << "new syntax." << endl;1655 DoLog(0) && (Log() << Verbose(0) << "new syntax." << endl); 1521 1656 configuration.Load(ConfigFileName, BondGraphFileName, periode, molecules); 1522 1657 configPresent = present; 1523 1658 break; 1524 1659 case 0: 1525 Log() << Verbose(0) << "old syntax." << endl;1660 DoLog(0) && (Log() << Verbose(0) << "old syntax." << endl); 1526 1661 configuration.LoadOld(ConfigFileName, BondGraphFileName, periode, molecules); 1527 1662 configPresent = present; 1528 1663 break; 1529 1664 default: 1530 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); 1531 1666 configPresent = empty; 1532 1667 } … … 1552 1687 configuration.BG = new BondGraph(configuration.GetIsAngstroem()); 1553 1688 if ((!BondGraphFileName.empty()) && (configuration.BG->LoadBondLengthTable(BondGraphFileName))) { 1554 Log() << Verbose(0) << "Bond length table loaded successfully." << endl;1689 DoLog(0) && (Log() << Verbose(0) << "Bond length table loaded successfully." << endl); 1555 1690 } else { 1556 eLog() << Verbose(1) << "Bond length table loading failed." << endl;1691 DoeLog(1) && (eLog()<< Verbose(1) << "Bond length table loading failed." << endl); 1557 1692 } 1558 1693 } … … 1561 1696 argptr = 1; 1562 1697 do { 1563 Log() << Verbose(0) << "Current Command line argument: " << argv[argptr] << "." << endl;1698 DoLog(0) && (Log() << Verbose(0) << "Current Command line argument: " << argv[argptr] << "." << endl); 1564 1699 if (argv[argptr][0] == '-') { 1565 1700 argptr++; … … 1570 1705 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1571 1706 ExitFlag = 255; 1572 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); 1573 1708 performCriticalExit(); 1574 1709 } else { 1575 1710 SaveFlag = true; 1576 Log() << Verbose(1) << "Parsing xyz file for new atoms." << endl;1711 DoLog(1) && (Log() << Verbose(1) << "Parsing xyz file for new atoms." << endl); 1577 1712 if (!mol->AddXYZFile(argv[argptr])) 1578 Log() << Verbose(2) << "File not found." << endl;1713 DoLog(2) && (Log() << Verbose(2) << "File not found." << endl); 1579 1714 else { 1580 Log() << Verbose(2) << "File found and parsed." << endl;1715 DoLog(2) && (Log() << Verbose(2) << "File found and parsed." << endl); 1581 1716 configPresent = present; 1582 1717 } … … 1587 1722 if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3]))) { 1588 1723 ExitFlag = 255; 1589 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); 1590 1725 performCriticalExit(); 1591 1726 } else { 1592 1727 SaveFlag = true; 1593 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] << "), "); 1594 1729 first = new atom; 1595 1730 first->type = periode->FindElement(atoi(argv[argptr])); 1596 1731 if (first->type != NULL) 1597 Log() << Verbose(2) << "found element " << first->type->name << endl;1732 DoLog(2) && (Log() << Verbose(2) << "found element " << first->type->name << endl); 1598 1733 for (int i=NDIM;i--;) 1599 1734 first->x.x[i] = atof(argv[argptr+1+i]); … … 1603 1738 configPresent = present; 1604 1739 } else 1605 eLog() << Verbose(1) << "Could not find the specified element." << endl;1740 DoeLog(1) && (eLog()<< Verbose(1) << "Could not find the specified element." << endl); 1606 1741 argptr+=4; 1607 1742 } … … 1616 1751 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1617 1752 ExitFlag = 255; 1618 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); 1619 1754 performCriticalExit(); 1620 1755 } else { 1621 1756 configuration.basis = argv[argptr]; 1622 Log() << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl;1757 DoLog(1) && (Log() << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl); 1623 1758 argptr+=1; 1624 1759 } … … 1627 1762 if (ExitFlag == 0) ExitFlag = 1; 1628 1763 { 1629 Log() << Verbose(1) << "Depth-First-Search Analysis." << endl;1764 DoLog(1) && (Log() << Verbose(1) << "Depth-First-Search Analysis." << endl); 1630 1765 MoleculeLeafClass *Subgraphs = NULL; // list of subgraphs from DFS analysis 1631 1766 int *MinimumRingSize = new int[mol->AtomCount]; … … 1658 1793 break; 1659 1794 case 'I': 1660 Log() << Verbose(1) << "Dissecting molecular system into a set of disconnected subgraphs ... " << endl;1795 DoLog(1) && (Log() << Verbose(1) << "Dissecting molecular system into a set of disconnected subgraphs ... " << endl); 1661 1796 // @TODO rather do the dissection afterwards 1662 1797 molecules->DissectMoleculeIntoConnectedSubgraphs(periode, &configuration); … … 1669 1804 } 1670 1805 } 1671 if ( mol == NULL) {1806 if ((mol == NULL) && (!molecules->ListOfMolecules.empty())) { 1672 1807 mol = *(molecules->ListOfMolecules.begin()); 1673 mol->ActiveFlag = true; 1808 if (mol != NULL) 1809 mol->ActiveFlag = true; 1674 1810 } 1675 1811 break; 1676 1812 case 'C': 1677 if (ExitFlag == 0) ExitFlag = 1; 1678 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr][0] == '-') || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-')) { 1679 ExitFlag = 255; 1680 eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C <Z> <output> <bin output>" << endl; 1681 performCriticalExit(); 1682 } else { 1683 ofstream output(argv[argptr+1]); 1684 ofstream binoutput(argv[argptr+2]); 1685 const double radius = 5.; 1686 1687 // get the boundary 1688 class molecule *Boundary = NULL; 1689 class Tesselation *TesselStruct = NULL; 1690 const LinkedCell *LCList = NULL; 1691 // find biggest molecule 1692 int counter = 0; 1693 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) { 1694 if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) { 1695 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; 1696 1962 } 1697 counter++;1698 1963 } 1699 bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives"); 1700 counter = 0; 1701 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) { 1702 Actives[counter++] = (*BigFinder)->ActiveFlag; 1703 (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true; 1704 } 1705 LCList = new LinkedCell(Boundary, 2.*radius); 1706 element *elemental = periode->FindElement((const int) atoi(argv[argptr])); 1707 FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL); 1708 int ranges[NDIM] = {1,1,1}; 1709 CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges ); 1710 OutputCorrelationToSurface(&output, surfacemap); 1711 BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 0. ); 1712 OutputCorrelation ( &binoutput, binmap ); 1713 output.close(); 1714 binoutput.close(); 1715 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) 1716 (*BigFinder)->ActiveFlag = Actives[counter++]; 1717 Free(&Actives); 1718 delete(LCList); 1719 delete(TesselStruct); 1720 argptr+=3; 1721 } 1722 break; 1964 break; 1965 } 1723 1966 case 'E': 1724 1967 if (ExitFlag == 0) ExitFlag = 1; 1725 1968 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) { 1726 1969 ExitFlag = 255; 1727 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); 1728 1971 performCriticalExit(); 1729 1972 } else { 1730 1973 SaveFlag = true; 1731 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); 1732 1975 first = mol->FindAtom(atoi(argv[argptr])); 1733 1976 first->type = periode->FindElement(atoi(argv[argptr+1])); … … 1737 1980 case 'F': 1738 1981 if (ExitFlag == 0) ExitFlag = 1; 1739 if (argptr+6 >=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]))) { 1740 1989 ExitFlag = 255; 1741 eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <boundary> <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); 1742 1991 performCriticalExit(); 1743 1992 } else { 1744 1993 SaveFlag = true; 1745 Log() << Verbose(1) << "Filling Box with water molecules." << endl;1994 DoLog(1) && (Log() << Verbose(1) << "Filling Box with water molecules." << endl); 1746 1995 // construct water molecule 1747 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); 1748 2002 molecule *Filling = NULL; 1749 atom *second = NULL, *third = NULL;1750 // first = new atom();1751 // first->type = periode->FindElement(5);1752 // first->x.Zero();1753 // filler->AddAtom(first);1754 first = new atom();1755 first->type = periode->FindElement(1);1756 first->x.Init(0.441, -0.143, 0.);1757 filler->AddAtom(first);1758 second = new atom();1759 second->type = periode->FindElement(1);1760 second->x.Init(-0.464, 1.137, 0.0);1761 filler->AddAtom(second);1762 third = new atom();1763 third->type = periode->FindElement(8);1764 third->x.Init(-0.464, 0.177, 0.);1765 filler->AddAtom(third);1766 filler->AddBond(first, third, 1);1767 filler->AddBond(second, third, 1);1768 2003 // call routine 1769 2004 double distance[NDIM]; 1770 2005 for (int i=0;i<NDIM;i++) 1771 distance[i] = atof(argv[argptr+i ]);1772 Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atof(argv[argptr+5]), atoi(argv[argptr+6]));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])); 1773 2008 if (Filling != NULL) { 1774 2009 Filling->ActiveFlag = false; … … 1783 2018 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1784 2019 ExitFlag =255; 1785 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); 1786 2021 performCriticalExit(); 1787 2022 } else { 1788 Log() << Verbose(0) << "Parsing bonds from " << argv[argptr] << "." << endl;2023 DoLog(0) && (Log() << Verbose(0) << "Parsing bonds from " << argv[argptr] << "." << endl); 1789 2024 ifstream *input = new ifstream(argv[argptr]); 1790 2025 mol->CreateAdjacencyListFromDbondFile(input); … … 1793 2028 } 1794 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 1795 2059 case 'N': 1796 2060 if (ExitFlag == 0) ExitFlag = 1; 1797 2061 if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){ 1798 2062 ExitFlag = 255; 1799 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); 1800 2064 performCriticalExit(); 1801 2065 } else { … … 1805 2069 //string filename(argv[argptr+1]); 1806 2070 //filename.append(".csv"); 1807 Log() << Verbose(0) << "Evaluating non-convex envelope of biggest molecule.";1808 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); 1809 2073 // find biggest molecule 1810 2074 int counter = 0; … … 1816 2080 counter++; 1817 2081 } 1818 Log() << Verbose(1) << "Biggest molecule has " << Boundary->AtomCount << " atoms." << endl;2082 DoLog(1) && (Log() << Verbose(1) << "Biggest molecule has " << Boundary->AtomCount << " atoms." << endl); 1819 2083 start = clock(); 1820 2084 LCList = new LinkedCell(Boundary, atof(argv[argptr])*2.); … … 1823 2087 //FindDistributionOfEllipsoids(T, &LCList, N, number, filename.c_str()); 1824 2088 end = clock(); 1825 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); 1826 2090 delete(LCList); 1827 2091 delete(T); … … 1833 2097 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1834 2098 ExitFlag = 255; 1835 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); 1836 2100 performCriticalExit(); 1837 2101 } else { 1838 Log() << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;2102 DoLog(1) && (Log() << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl); 1839 2103 ofstream *output = new ofstream(argv[argptr], ios::trunc); 1840 2104 if (!mol->OutputTemperatureFromTrajectories(output, 0, mol->MDSteps)) 1841 Log() << Verbose(2) << "File could not be written." << endl;2105 DoLog(2) && (Log() << Verbose(2) << "File could not be written." << endl); 1842 2106 else 1843 Log() << Verbose(2) << "File stored." << endl;2107 DoLog(2) && (Log() << Verbose(2) << "File stored." << endl); 1844 2108 output->close(); 1845 2109 delete(output); … … 1851 2115 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1852 2116 ExitFlag = 255; 1853 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); 1854 2118 performCriticalExit(); 1855 2119 } else { 1856 2120 SaveFlag = true; 1857 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); 1858 2122 if (atoi(argv[argptr+3]) == 1) 1859 Log() << Verbose(1) << "Using Identity for the permutation map." << endl;2123 DoLog(1) && (Log() << Verbose(1) << "Using Identity for the permutation map." << endl); 1860 2124 if (!mol->LinearInterpolationBetweenConfiguration(atoi(argv[argptr]), atoi(argv[argptr+1]), argv[argptr+2], configuration, atoi(argv[argptr+3])) == 1 ? true : false) 1861 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); 1862 2126 else 1863 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); 1864 2128 argptr+=4; 1865 2129 } … … 1869 2133 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1870 2134 ExitFlag = 255; 1871 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); 1872 2136 performCriticalExit(); 1873 2137 } else { 1874 2138 SaveFlag = true; 1875 Log() << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;2139 DoLog(1) && (Log() << Verbose(1) << "Parsing forces file and Verlet integrating." << endl); 1876 2140 if (!mol->VerletForceIntegration(argv[argptr], configuration)) 1877 Log() << Verbose(2) << "File not found." << endl;2141 DoLog(2) && (Log() << Verbose(2) << "File not found." << endl); 1878 2142 else 1879 Log() << Verbose(2) << "File found and parsed." << endl;2143 DoLog(2) && (Log() << Verbose(2) << "File found and parsed." << endl); 1880 2144 argptr+=1; 1881 2145 } … … 1885 2149 if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) { 1886 2150 ExitFlag = 255; 1887 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); 1888 2152 performCriticalExit(); 1889 2153 } else { 1890 2154 SaveFlag = true; 1891 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); 1892 2156 double tmp1 = atof(argv[argptr+1]); 1893 2157 atom *third = mol->FindAtom(atoi(argv[argptr])); … … 1902 2166 } 1903 2167 } else { 1904 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); 1905 2169 } 1906 2170 argptr+=2; … … 1911 2175 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1912 2176 ExitFlag = 255; 1913 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); 1914 2178 performCriticalExit(); 1915 2179 } else { 1916 2180 if (ExitFlag == 0) ExitFlag = 1; 1917 2181 SaveFlag = true; 1918 Log() << Verbose(1) << "Translating all ions by given vector." << endl;2182 DoLog(1) && (Log() << Verbose(1) << "Translating all ions by given vector." << endl); 1919 2183 for (int i=NDIM;i--;) 1920 2184 x.x[i] = atof(argv[argptr+i]); … … 1927 2191 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1928 2192 ExitFlag = 255; 1929 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); 1930 2194 performCriticalExit(); 1931 2195 } else { 1932 2196 if (ExitFlag == 0) ExitFlag = 1; 1933 2197 SaveFlag = true; 1934 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); 1935 2199 for (int i=NDIM;i--;) 1936 2200 x.x[i] = atof(argv[argptr+i]); … … 1943 2207 if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1944 2208 ExitFlag = 255; 1945 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); 1946 2210 performCriticalExit(); 1947 2211 } else { 1948 2212 SaveFlag = true; 1949 2213 j = -1; 1950 Log() << Verbose(1) << "Scaling all ion positions by factor." << endl;2214 DoLog(1) && (Log() << Verbose(1) << "Scaling all ion positions by factor." << endl); 1951 2215 factor = new double[NDIM]; 1952 2216 factor[0] = atof(argv[argptr]); … … 1954 2218 factor[2] = atof(argv[argptr+2]); 1955 2219 mol->Scale((const double ** const)&factor); 2220 double * const cell_size = World::get()->cell_size; 1956 2221 for (int i=0;i<NDIM;i++) { 1957 2222 j += i+1; 1958 2223 x.x[i] = atof(argv[NDIM+i]); 1959 mol->cell_size[j]*=factor[i];2224 cell_size[j]*=factor[i]; 1960 2225 } 1961 2226 delete[](factor); … … 1967 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])) ) { 1968 2233 ExitFlag = 255; 1969 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); 1970 2235 performCriticalExit(); 1971 2236 } else { 1972 2237 SaveFlag = true; 1973 2238 j = -1; 1974 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; 1975 2241 for (int i=0;i<6;i++) { 1976 mol->cell_size[i] = atof(argv[argptr+i]);2242 cell_size[i] = atof(argv[argptr+i]); 1977 2243 } 1978 2244 // center … … 1985 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])) ) { 1986 2252 ExitFlag = 255; 1987 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); 1988 2254 performCriticalExit(); 1989 2255 } else { 1990 2256 SaveFlag = true; 1991 2257 j = -1; 1992 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; 1993 2260 for (int i=0;i<6;i++) { 1994 mol->cell_size[i] = atof(argv[argptr+i]);2261 cell_size[i] = atof(argv[argptr+i]); 1995 2262 } 1996 2263 // center … … 2003 2270 if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 2004 2271 ExitFlag = 255; 2005 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); 2006 2273 performCriticalExit(); 2007 2274 } else { 2008 2275 SaveFlag = true; 2009 2276 j = -1; 2010 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); 2011 2278 // make every coordinate positive 2012 2279 mol->CenterEdge(&x); … … 2014 2281 mol->SetBoxDimension(&x); 2015 2282 // translate each coordinate by boundary 2283 double * const cell_size = World::get()->cell_size; 2016 2284 j=-1; 2017 2285 for (int i=0;i<NDIM;i++) { 2018 2286 j += i+1; 2019 2287 x.x[i] = atof(argv[argptr+i]); 2020 mol->cell_size[j] += x.x[i]*2.;2288 cell_size[j] += x.x[i]*2.; 2021 2289 } 2022 2290 mol->Translate((const Vector *)&x); … … 2027 2295 if (ExitFlag == 0) ExitFlag = 1; 2028 2296 SaveFlag = true; 2029 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); 2030 2298 x.Zero(); 2031 2299 mol->CenterEdge(&x); … … 2037 2305 if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr]))) { 2038 2306 ExitFlag = 255; 2039 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); 2040 2308 performCriticalExit(); 2041 2309 } else { 2042 2310 SaveFlag = true; 2043 Log() << Verbose(1) << "Removing atom " << argv[argptr] << "." << endl;2311 DoLog(1) && (Log() << Verbose(1) << "Removing atom " << argv[argptr] << "." << endl); 2044 2312 atom *first = mol->FindAtom(atoi(argv[argptr])); 2045 2313 mol->RemoveAtom(first); … … 2051 2319 if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) { 2052 2320 ExitFlag = 255; 2053 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); 2054 2322 performCriticalExit(); 2055 2323 } else { 2056 Log() << Verbose(0) << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;2057 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); 2058 2326 start = clock(); 2059 2327 mol->CreateAdjacencyList(atof(argv[argptr++]), configuration.GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL); 2060 Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;2328 DoLog(0) && (Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl); 2061 2329 if (mol->first->next != mol->last) { 2062 2330 ExitFlag = mol->FragmentMolecule(atoi(argv[argptr]), &configuration); 2063 2331 } 2064 2332 end = clock(); 2065 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); 2066 2334 argptr+=2; 2067 2335 } … … 2071 2339 j = atoi(argv[argptr++]); 2072 2340 if ((j<0) || (j>1)) { 2073 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); 2074 2342 j = 0; 2075 2343 } 2076 2344 if (j) { 2077 2345 SaveFlag = true; 2078 Log() << Verbose(0) << "Converting to prinicipal axis system." << endl;2346 DoLog(0) && (Log() << Verbose(0) << "Converting to prinicipal axis system." << endl); 2079 2347 } else 2080 Log() << Verbose(0) << "Evaluating prinicipal axis." << endl;2348 DoLog(0) && (Log() << Verbose(0) << "Evaluating prinicipal axis." << endl); 2081 2349 mol->PrincipalAxisSystem((bool)j); 2082 2350 break; … … 2085 2353 if ((argptr+1 >= argc) || (argv[argptr][0] == '-')){ 2086 2354 ExitFlag = 255; 2087 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); 2088 2356 performCriticalExit(); 2089 2357 } else { 2090 2358 class Tesselation *TesselStruct = NULL; 2091 2359 const LinkedCell *LCList = NULL; 2092 Log() << Verbose(0) << "Evaluating volume of the convex envelope.";2093 Log() << Verbose(1) << "Storing tecplot convex data in " << argv[argptr] << "." << endl;2094 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); 2095 2363 LCList = new LinkedCell(mol, 10.); 2096 2364 //FindConvexBorder(mol, LCList, argv[argptr]); … … 2099 2367 double volumedifference = ConvexizeNonconvexEnvelope(TesselStruct, mol, argv[argptr]); 2100 2368 double clustervolume = VolumeOfConvexEnvelope(TesselStruct, &configuration); 2101 Log() << Verbose(0) << "The tesselated volume area is " << clustervolume << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;2102 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); 2103 2371 delete(TesselStruct); 2104 2372 delete(LCList); … … 2110 2378 if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) { 2111 2379 ExitFlag = 255; 2112 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); 2113 2381 performCriticalExit(); 2114 2382 } else { 2115 2383 volume = atof(argv[argptr++]); 2116 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); 2117 2385 } 2118 2386 case 'u': … … 2121 2389 if (volume != -1) 2122 2390 ExitFlag = 255; 2123 eLog() << Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl;2391 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl); 2124 2392 performCriticalExit(); 2125 2393 } else { 2126 2394 double density; 2127 2395 SaveFlag = true; 2128 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."); 2129 2397 density = atof(argv[argptr++]); 2130 2398 if (density < 1.0) { 2131 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); 2132 2400 density = 1.3; 2133 2401 } … … 2135 2403 // repetition[i] = atoi(argv[argptr++]); 2136 2404 // if (repetition[i] < 1) 2137 // eLog() << Verbose(1) << "repetition value must be greater 1!" << endl;2405 // DoeLog(1) && (eLog()<< Verbose(1) << "repetition value must be greater 1!" << endl); 2138 2406 // repetition[i] = 1; 2139 2407 // } … … 2145 2413 if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 2146 2414 ExitFlag = 255; 2147 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); 2148 2416 performCriticalExit(); 2149 2417 } else { 2150 2418 SaveFlag = true; 2419 double * const cell_size = World::get()->cell_size; 2151 2420 for (int axis = 1; axis <= NDIM; axis++) { 2152 2421 int faktor = atoi(argv[argptr++]); … … 2155 2424 Vector ** vectors; 2156 2425 if (faktor < 1) { 2157 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); 2158 2427 faktor = 1; 2159 2428 } … … 2172 2441 } 2173 2442 if (count != j) 2174 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); 2175 2444 x.Zero(); 2176 2445 y.Zero(); 2177 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 2178 2447 for (int i=1;i<faktor;i++) { // then add this list with respective translation factor times 2179 2448 x.AddVector(&y); // per factor one cell width further … … 2196 2465 mol->Translate(&x); 2197 2466 } 2198 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; 2199 2468 } 2200 2469 } … … 2213 2482 } else { // no arguments, hence scan the elements db 2214 2483 if (periode->LoadPeriodentafel(configuration.databasepath)) 2215 Log() << Verbose(0) << "Element list loaded successfully." << endl;2484 DoLog(0) && (Log() << Verbose(0) << "Element list loaded successfully." << endl); 2216 2485 else 2217 Log() << Verbose(0) << "Element list loading failed." << endl;2486 DoLog(0) && (Log() << Verbose(0) << "Element list loading failed." << endl); 2218 2487 configuration.RetrieveConfigPathAndName("main_pcp_linux"); 2219 2488 } … … 2238 2507 2239 2508 cout << ESPACKVersion << endl; 2509 2510 DoLog(1) && (Log() << Verbose(1) << "test" << endl); 2511 DoLog(3) && (Log() << Verbose(1) << "test"); 2240 2512 2241 2513 setVerbosity(0); … … 2265 2537 if (molecules->ListOfMolecules.size() == 0) { 2266 2538 mol = new molecule(periode); 2267 if (mol->cell_size[0] == 0.) { 2268 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); 2269 2542 for (int i=0;i<6;i++) { 2270 Log() << Verbose(1) << "Cell size" << i << ": ";2271 cin >> mol->cell_size[i];2543 DoLog(1) && (Log() << Verbose(1) << "Cell size" << i << ": "); 2544 cin >> cell_size[i]; 2272 2545 } 2273 2546 } … … 2279 2552 2280 2553 // now the main construction loop 2281 Log() << Verbose(0) << endl << "Now comes the real construction..." << endl;2554 DoLog(0) && (Log() << Verbose(0) << endl << "Now comes the real construction..." << endl); 2282 2555 do { 2283 Log() << Verbose(0) << endl << endl;2284 Log() << Verbose(0) << "============Molecule list=======================" << endl;2556 DoLog(0) && (Log() << Verbose(0) << endl << endl); 2557 DoLog(0) && (Log() << Verbose(0) << "============Molecule list=======================" << endl); 2285 2558 molecules->Enumerate((ofstream *)&cout); 2286 Log() << Verbose(0) << "============Menu===============================" << endl;2287 Log() << Verbose(0) << "a - set molecule (in)active" << endl;2288 Log() << Verbose(0) << "e - edit molecules (load, parse, save)" << endl;2289 Log() << Verbose(0) << "g - globally manipulate atoms in molecule" << endl;2290 Log() << Verbose(0) << "M - Merge molecules" << endl;2291 Log() << Verbose(0) << "m - manipulate atoms" << endl;2292 Log() << Verbose(0) << "-----------------------------------------------" << endl;2293 Log() << Verbose(0) << "c - edit the current configuration" << endl;2294 Log() << Verbose(0) << "-----------------------------------------------" << endl;2295 Log() << Verbose(0) << "s - save current setup to config file" << endl;2296 Log() << Verbose(0) << "T - call the current test routine" << endl;2297 Log() << Verbose(0) << "q - quit" << endl;2298 Log() << Verbose(0) << "===============================================" << endl;2299 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: "); 2300 2573 cin >> choice; 2301 2574 … … 2303 2576 case 'a': // (in)activate molecule 2304 2577 { 2305 Log() << Verbose(0) << "Enter index of molecule: ";2578 DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: "); 2306 2579 cin >> j; 2307 2580 for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) … … 2349 2622 // save element data base 2350 2623 if (periode->StorePeriodentafel(configuration->databasepath)) //ElementsFileName 2351 Log() << Verbose(0) << "Saving of elements.db successful." << endl;2624 DoLog(0) && (Log() << Verbose(0) << "Saving of elements.db successful." << endl); 2352 2625 else 2353 Log() << Verbose(0) << "Saving of elements.db failed." << endl;2626 DoLog(0) && (Log() << Verbose(0) << "Saving of elements.db failed." << endl); 2354 2627 2355 2628 delete(molecules); // also free's all molecules contained
Note:
See TracChangeset
for help on using the changeset viewer.
