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