Changeset 375b458 for src/builder.cpp
- Timestamp:
- Jul 23, 2009, 1:53:01 PM (16 years ago)
- Branches:
- Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
- Children:
- 51c910
- Parents:
- a5b2c3a
- File:
-
- 1 edited
-
src/builder.cpp (modified) (22 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/builder.cpp
ra5b2c3a r375b458 15 15 * \section about About the Program 16 16 * 17 * 18 * 19 * 17 * Molecuilder is a short program, written in C++, that enables the construction of a coordinate set for the 18 * atoms making up an molecule by the successive statement of binding angles and distances and referencing to 19 * already constructed atoms. 20 20 * 21 * 22 * 21 * A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello 22 * molecular dynamics implementation. 23 23 * 24 24 * \section install Installation 25 25 * 26 * 27 * 28 * 29 * 26 * Installation should without problems succeed as follows: 27 * -# ./configure (or: mkdir build;mkdir run;cd build; ../configure --bindir=../run) 28 * -# make 29 * -# make install 30 30 * 31 * 32 * 33 * 31 * Further useful commands are 32 * -# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n 33 * -# make doxygen-doc: Creates these html pages out of the documented source 34 34 * 35 35 * \section run Running 36 36 * 37 * 37 * The program can be executed by running: ./molecuilder 38 38 * 39 * 40 * 41 * 39 * Note, that it uses a database, called "elements.db", in the executable's directory. If the file is not found, 40 * it is created and any given data on elements of the periodic table will be stored therein and re-used on 41 * later re-execution. 42 42 * 43 43 * \section ref References 44 44 * 45 * 45 * For the special configuration file format, see the documentation of pcp. 46 46 * 47 47 */ … … 63 63 static void AddAtoms(periodentafel *periode, molecule *mol) 64 64 { 65 66 67 Vector x,y,z,n;// coordinates for absolute point in cell volume68 69 char choice;// menu choice char70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 65 atom *first, *second, *third, *fourth; 66 Vector **atoms; 67 Vector x,y,z,n; // coordinates for absolute point in cell volume 68 double a,b,c; 69 char choice; // menu choice char 70 bool valid; 71 72 cout << Verbose(0) << "===========ADD ATOM============================" << endl; 73 cout << Verbose(0) << " a - state absolute coordinates of atom" << endl; 74 cout << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl; 75 cout << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl; 76 cout << Verbose(0) << " d - state two atoms, two angles and a distance" << endl; 77 cout << Verbose(0) << " e - least square distance position to a set of atoms" << endl; 78 cout << Verbose(0) << "all else - go back" << endl; 79 cout << Verbose(0) << "===============================================" << endl; 80 cout << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl; 81 cout << Verbose(0) << "INPUT: "; 82 cin >> choice; 83 84 switch (choice) { 85 85 default: 86 86 cout << Verbose(0) << "Not a valid choice." << endl; 87 87 break; 88 88 case 'a': // absolute coordinates of atom 89 89 cout << Verbose(0) << "Enter absolute coordinates." << endl; 90 90 first = new atom; 91 91 first->x.AskPosition(mol->cell_size, false); 92 first->type = periode->AskElement(); 93 mol->AddAtom(first); 94 95 96 92 first->type = periode->AskElement(); // give type 93 mol->AddAtom(first); // add to molecule 94 break; 95 96 case 'b': // relative coordinates of atom wrt to reference point 97 97 first = new atom; 98 98 valid = true; … … 106 106 cout << Verbose(0) << "\n"; 107 107 } while (!(valid = mol->CheckBounds((const Vector *)&first->x))); 108 first->type = periode->AskElement(); 109 mol->AddAtom(first); 110 111 112 108 first->type = periode->AskElement(); // give type 109 mol->AddAtom(first); // add to molecule 110 break; 111 112 case 'c': // relative coordinates of atom wrt to already placed atom 113 113 first = new atom; 114 114 valid = true; … … 122 122 } 123 123 } while (!(valid = mol->CheckBounds((const Vector *)&first->x))); 124 first->type = periode->AskElement(); 125 mol->AddAtom(first); 124 first->type = periode->AskElement(); // give type 125 mol->AddAtom(first); // add to molecule 126 126 break; 127 127 … … 224 224 cout << Verbose(0) << endl; 225 225 } while (!(valid = mol->CheckBounds((const Vector *)&first->x))); 226 first->type = periode->AskElement(); 227 mol->AddAtom(first); 228 229 230 226 first->type = periode->AskElement(); // give type 227 mol->AddAtom(first); // add to molecule 228 break; 229 230 case 'e': // least square distance position to a set of atoms 231 231 first = new atom; 232 232 atoms = new (Vector*[128]); … … 248 248 249 249 first->x.Output((ofstream *)&cout); 250 first->type = periode->AskElement(); 251 mol->AddAtom(first); 250 first->type = periode->AskElement(); // give type 251 mol->AddAtom(first); // add to molecule 252 252 } else { 253 253 delete first; 254 254 cout << Verbose(0) << "Please enter at least two vectors!\n"; 255 255 } 256 257 256 break; 257 }; 258 258 }; 259 259 … … 263 263 static void CenterAtoms(molecule *mol) 264 264 { 265 266 char choice;// menu choice char267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 mol->CenterEdge((ofstream *)&cout, &x);// make every coordinate positive297 298 299 300 301 mol->SetBoxDimension(&helper);// update Box of atoms by boundary302 303 304 305 306 307 308 309 310 311 312 313 314 265 Vector x, y, helper; 266 char choice; // menu choice char 267 268 cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl; 269 cout << Verbose(0) << " a - on origin" << endl; 270 cout << Verbose(0) << " b - on center of gravity" << endl; 271 cout << Verbose(0) << " c - within box with additional boundary" << endl; 272 cout << Verbose(0) << " d - within given simulation box" << endl; 273 cout << Verbose(0) << "all else - go back" << endl; 274 cout << Verbose(0) << "===============================================" << endl; 275 cout << Verbose(0) << "INPUT: "; 276 cin >> choice; 277 278 switch (choice) { 279 default: 280 cout << Verbose(0) << "Not a valid choice." << endl; 281 break; 282 case 'a': 283 cout << Verbose(0) << "Centering atoms in config file on origin." << endl; 284 mol->CenterOrigin((ofstream *)&cout, &x); 285 break; 286 case 'b': 287 cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl; 288 mol->CenterGravity((ofstream *)&cout, &x); 289 break; 290 case 'c': 291 cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl; 292 for (int i=0;i<NDIM;i++) { 293 cout << Verbose(0) << "Enter axis " << i << " boundary: "; 294 cin >> y.x[i]; 295 } 296 mol->CenterEdge((ofstream *)&cout, &x); // make every coordinate positive 297 mol->Translate(&y); // translate by boundary 298 helper.CopyVector(&y); 299 helper.Scale(2.); 300 helper.AddVector(&x); 301 mol->SetBoxDimension(&helper); // update Box of atoms by boundary 302 break; 303 case 'd': 304 cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl; 305 for (int i=0;i<NDIM;i++) { 306 cout << Verbose(0) << "Enter axis " << i << " boundary: "; 307 cin >> x.x[i]; 308 } 309 // center 310 mol->CenterInBox((ofstream *)&cout, &x); 311 // update Box of atoms by boundary 312 mol->SetBoxDimension(&x); 313 break; 314 } 315 315 }; 316 316 … … 321 321 static void AlignAtoms(periodentafel *periode, molecule *mol) 322 322 { 323 324 325 char choice;// menu choice char326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 323 atom *first, *second, *third; 324 Vector x,n; 325 char choice; // menu choice char 326 327 cout << Verbose(0) << "===========ALIGN ATOMS=========================" << endl; 328 cout << Verbose(0) << " a - state three atoms defining align plane" << endl; 329 cout << Verbose(0) << " b - state alignment vector" << endl; 330 cout << Verbose(0) << " c - state two atoms in alignment direction" << endl; 331 cout << Verbose(0) << " d - align automatically by least square fit" << endl; 332 cout << Verbose(0) << "all else - go back" << endl; 333 cout << Verbose(0) << "===============================================" << endl; 334 cout << Verbose(0) << "INPUT: "; 335 cin >> choice; 336 337 switch (choice) { 338 default: 339 case 'a': // three atoms defining mirror plane 340 first = mol->AskAtom("Enter first atom: "); 341 second = mol->AskAtom("Enter second atom: "); 342 third = mol->AskAtom("Enter third atom: "); 343 344 n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x); 345 break; 346 case 'b': // normal vector of mirror plane 347 cout << Verbose(0) << "Enter normal vector of mirror plane." << endl; 348 n.AskPosition(mol->cell_size,0); 349 n.Normalize(); 350 break; 351 case 'c': // three atoms defining mirror plane 352 first = mol->AskAtom("Enter first atom: "); 353 second = mol->AskAtom("Enter second atom: "); 354 355 n.CopyVector((const Vector *)&first->x); 356 n.SubtractVector((const Vector *)&second->x); 357 n.Normalize(); 358 break; 359 case 'd': 360 char shorthand[4]; 361 Vector a; 362 struct lsq_params param; 363 do { 364 fprintf(stdout, "Enter the element of atoms to be chosen: "); 365 fscanf(stdin, "%3s", shorthand); 366 } while ((param.type = periode->FindElement(shorthand)) == NULL); 367 cout << Verbose(0) << "Element is " << param.type->name << endl; 368 mol->GetAlignvector(¶m); 369 for (int i=NDIM;i--;) { 370 x.x[i] = gsl_vector_get(param.x,i); 371 n.x[i] = gsl_vector_get(param.x,i+NDIM); 372 } 373 gsl_vector_free(param.x); 374 cout << Verbose(0) << "Offset vector: "; 375 x.Output((ofstream *)&cout); 376 cout << Verbose(0) << endl; 377 n.Normalize(); 378 break; 379 }; 380 cout << Verbose(0) << "Alignment vector: "; 381 n.Output((ofstream *)&cout); 382 cout << Verbose(0) << endl; 383 mol->Align(&n); 384 384 }; 385 385 … … 389 389 static void MirrorAtoms(molecule *mol) 390 390 { 391 392 393 char choice;// menu choice char394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 391 atom *first, *second, *third; 392 Vector n; 393 char choice; // menu choice char 394 395 cout << Verbose(0) << "===========MIRROR ATOMS=========================" << endl; 396 cout << Verbose(0) << " a - state three atoms defining mirror plane" << endl; 397 cout << Verbose(0) << " b - state normal vector of mirror plane" << endl; 398 cout << Verbose(0) << " c - state two atoms in normal direction" << endl; 399 cout << Verbose(0) << "all else - go back" << endl; 400 cout << Verbose(0) << "===============================================" << endl; 401 cout << Verbose(0) << "INPUT: "; 402 cin >> choice; 403 404 switch (choice) { 405 default: 406 case 'a': // three atoms defining mirror plane 407 first = mol->AskAtom("Enter first atom: "); 408 second = mol->AskAtom("Enter second atom: "); 409 third = mol->AskAtom("Enter third atom: "); 410 411 n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x); 412 break; 413 case 'b': // normal vector of mirror plane 414 cout << Verbose(0) << "Enter normal vector of mirror plane." << endl; 415 n.AskPosition(mol->cell_size,0); 416 n.Normalize(); 417 break; 418 case 'c': // three atoms defining mirror plane 419 first = mol->AskAtom("Enter first atom: "); 420 second = mol->AskAtom("Enter second atom: "); 421 422 n.CopyVector((const Vector *)&first->x); 423 n.SubtractVector((const Vector *)&second->x); 424 n.Normalize(); 425 break; 426 }; 427 cout << Verbose(0) << "Normal vector: "; 428 n.Output((ofstream *)&cout); 429 cout << Verbose(0) << endl; 430 mol->Mirror((const Vector *)&n); 431 431 }; 432 432 … … 436 436 static void RemoveAtoms(molecule *mol) 437 437 { 438 atom *first, *second, *third; 439 int axis; 440 double tmp1, tmp2; 441 char choice; // menu choice char 442 443 cout << Verbose(0) << "===========REMOVE ATOMS=========================" << endl; 444 cout << Verbose(0) << " a - state atom for removal by number" << endl; 445 cout << Verbose(0) << " b - keep only in radius around atom" << endl; 446 cout << Verbose(0) << " c - remove this with one axis greater value" << endl; 447 cout << Verbose(0) << "all else - go back" << endl; 448 cout << Verbose(0) << "===============================================" << endl; 449 cout << Verbose(0) << "INPUT: "; 450 cin >> choice; 451 452 switch (choice) { 453 default: 454 case 'a': 455 if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: "))) 456 cout << Verbose(1) << "Atom removed." << endl; 457 else 458 cout << Verbose(1) << "Atom not found." << endl; 459 break; 460 case 'b': 461 third = mol->AskAtom("Enter number of atom as reference point: "); 462 cout << Verbose(0) << "Enter radius: "; 463 cin >> tmp1; 464 first = mol->start; 465 second = first->next; 466 while(second != mol->end) { 467 first = second; 468 second = first->next; 469 if (first->x.DistanceSquared((const Vector *)&third->x) > tmp1*tmp1) // distance to first above radius ... 470 mol->RemoveAtom(first); 471 } 472 break; 473 case 'c': 474 cout << Verbose(0) << "Which axis is it: "; 475 cin >> axis; 476 cout << Verbose(0) << "Lower boundary: "; 477 cin >> tmp1; 478 cout << Verbose(0) << "Upper boundary: "; 479 cin >> tmp2; 480 first = mol->start; 438 atom *first, *second, *third; 439 int axis; 440 double tmp1, tmp2; 441 char choice; // menu choice char 442 443 cout << Verbose(0) << "===========REMOVE ATOMS=========================" << endl; 444 cout << Verbose(0) << " a - state atom for removal by number" << endl; 445 cout << Verbose(0) << " b - keep only in radius around atom" << endl; 446 cout << Verbose(0) << " c - remove this with one axis greater value" << endl; 447 cout << Verbose(0) << "all else - go back" << endl; 448 cout << Verbose(0) << "===============================================" << endl; 449 cout << Verbose(0) << "INPUT: "; 450 cin >> choice; 451 452 switch (choice) { 453 default: 454 case 'a': 455 if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: "))) 456 cout << Verbose(1) << "Atom removed." << endl; 457 else 458 cout << Verbose(1) << "Atom not found." << endl; 459 break; 460 case 'b': 461 third = mol->AskAtom("Enter number of atom as reference point: "); 462 cout << Verbose(0) << "Enter radius: "; 463 cin >> tmp1; 464 first = mol->start; 481 465 second = first->next; 482 466 while(second != mol->end) { 483 467 first = second; 484 468 second = first->next; 485 if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ... 486 //cout << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl; 487 mol->RemoveAtom(first); 488 } 489 } 490 break; 491 }; 492 //mol->Output((ofstream *)&cout); 493 choice = 'r'; 469 if (first->x.DistanceSquared((const Vector *)&third->x) > tmp1*tmp1) // distance to first above radius ... 470 mol->RemoveAtom(first); 471 } 472 break; 473 case 'c': 474 cout << Verbose(0) << "Which axis is it: "; 475 cin >> axis; 476 cout << Verbose(0) << "Lower boundary: "; 477 cin >> tmp1; 478 cout << Verbose(0) << "Upper boundary: "; 479 cin >> tmp2; 480 first = mol->start; 481 second = first->next; 482 while(second != mol->end) { 483 first = second; 484 second = first->next; 485 if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ... 486 //cout << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl; 487 mol->RemoveAtom(first); 488 } 489 } 490 break; 491 }; 492 //mol->Output((ofstream *)&cout); 493 choice = 'r'; 494 494 }; 495 495 … … 500 500 static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration) 501 501 { 502 503 504 505 506 char choice;// menu choice char507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 third= mol->AskAtom("Enter last atom: ");565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 502 atom *first, *second, *third; 503 Vector x,y; 504 double min[256], tmp1, tmp2, tmp3; 505 int Z; 506 char choice; // menu choice char 507 508 cout << Verbose(0) << "===========MEASURE ATOMS=========================" << endl; 509 cout << Verbose(0) << " a - calculate bond length between one atom and all others" << endl; 510 cout << Verbose(0) << " b - calculate bond length between two atoms" << endl; 511 cout << Verbose(0) << " c - calculate bond angle" << endl; 512 cout << Verbose(0) << " d - calculate principal axis of the system" << endl; 513 cout << Verbose(0) << " e - calculate volume of the convex envelope" << endl; 514 cout << Verbose(0) << " f - calculate temperature from current velocity" << endl; 515 cout << Verbose(0) << " g - output all temperatures per step from velocities" << endl; 516 cout << Verbose(0) << "all else - go back" << endl; 517 cout << Verbose(0) << "===============================================" << endl; 518 cout << Verbose(0) << "INPUT: "; 519 cin >> choice; 520 521 switch(choice) { 522 default: 523 cout << Verbose(1) << "Not a valid choice." << endl; 524 break; 525 case 'a': 526 first = mol->AskAtom("Enter first atom: "); 527 for (int i=MAX_ELEMENTS;i--;) 528 min[i] = 0.; 529 530 second = mol->start; 531 while ((second->next != mol->end)) { 532 second = second->next; // advance 533 Z = second->type->Z; 534 tmp1 = 0.; 535 if (first != second) { 536 x.CopyVector((const Vector *)&first->x); 537 x.SubtractVector((const Vector *)&second->x); 538 tmp1 = x.Norm(); 539 } 540 if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1; 541 //cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl; 542 } 543 for (int i=MAX_ELEMENTS;i--;) 544 if (min[i] != 0.) cout << Verbose(0) << "Minimum Bond length between " << first->type->name << " Atom " << first->nr << " and next Ion of type " << (periode->FindElement(i))->name << ": " << min[i] << " a.u." << endl; 545 break; 546 547 case 'b': 548 first = mol->AskAtom("Enter first atom: "); 549 second = mol->AskAtom("Enter second atom: "); 550 for (int i=NDIM;i--;) 551 min[i] = 0.; 552 x.CopyVector((const Vector *)&first->x); 553 x.SubtractVector((const Vector *)&second->x); 554 tmp1 = x.Norm(); 555 cout << Verbose(1) << "Distance vector is "; 556 x.Output((ofstream *)&cout); 557 cout << "." << endl << "Norm of distance is " << tmp1 << "." << endl; 558 break; 559 560 case 'c': 561 cout << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl; 562 first = mol->AskAtom("Enter first atom: "); 563 second = mol->AskAtom("Enter central atom: "); 564 third = mol->AskAtom("Enter last atom: "); 565 tmp1 = tmp2 = tmp3 = 0.; 566 x.CopyVector((const Vector *)&first->x); 567 x.SubtractVector((const Vector *)&second->x); 568 y.CopyVector((const Vector *)&third->x); 569 y.SubtractVector((const Vector *)&second->x); 570 cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": "; 571 cout << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl; 572 break; 573 case 'd': 574 cout << Verbose(0) << "Evaluating prinicipal axis." << endl; 575 cout << Verbose(0) << "Shall we rotate? [0/1]: "; 576 cin >> Z; 577 if ((Z >=0) && (Z <=1)) 578 mol->PrincipalAxisSystem((ofstream *)&cout, (bool)Z); 579 else 580 mol->PrincipalAxisSystem((ofstream *)&cout, false); 581 break; 582 case 'e': 583 cout << Verbose(0) << "Evaluating volume of the convex envelope."; 584 VolumeOfConvexEnvelope((ofstream *)&cout, NULL, configuration, NULL, mol); 585 break; 586 case 'f': 587 mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps, (ofstream *)&cout); 588 break; 589 case 'g': 590 { 591 char filename[255]; 592 cout << "Please enter filename: " << endl; 593 cin >> filename; 594 cout << Verbose(1) << "Storing temperatures in " << filename << "." << endl; 595 ofstream *output = new ofstream(filename, ios::trunc); 596 if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output)) 597 cout << Verbose(2) << "File could not be written." << endl; 598 else 599 cout << Verbose(2) << "File stored." << endl; 600 output->close(); 601 delete(output); 602 } 603 break; 604 } 605 605 }; 606 606 … … 611 611 static void FragmentAtoms(molecule *mol, config *configuration) 612 612 { 613 614 615 616 617 618 619 if (mol->first->next != mol->last) {// there are bonds620 621 622 623 624 625 613 int Order1; 614 clock_t start, end; 615 616 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl; 617 cout << Verbose(0) << "What's the desired bond order: "; 618 cin >> Order1; 619 if (mol->first->next != mol->last) { // there are bonds 620 start = clock(); 621 mol->FragmentMolecule((ofstream *)&cout, Order1, configuration); 622 end = clock(); 623 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl; 624 } else 625 cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl; 626 626 }; 627 627 … … 1053 1053 static void testroutine(MoleculeListClass *molecules) 1054 1054 { 1055 1056 1055 // the current test routine checks the functionality of the KeySet&Graph concept: 1056 // We want to have a multiindex (the KeySet) describing a unique subgraph 1057 1057 int i, comp, counter=0; 1058 1058 … … 1067 1067 atom *Walker = mol->start; 1068 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 while (A !=Subgraphs.end()) {1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1069 // generate some KeySets 1070 cout << "Generating KeySets." << endl; 1071 KeySet TestSets[mol->AtomCount+1]; 1072 i=1; 1073 while (Walker->next != mol->end) { 1074 Walker = Walker->next; 1075 for (int j=0;j<i;j++) { 1076 TestSets[j].insert(Walker->nr); 1077 } 1078 i++; 1079 } 1080 cout << "Testing insertion of already present item in KeySets." << endl; 1081 KeySetTestPair test; 1082 test = TestSets[mol->AtomCount-1].insert(Walker->nr); 1083 if (test.second) { 1084 cout << Verbose(1) << "Insertion worked?!" << endl; 1085 } else { 1086 cout << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl; 1087 } 1088 TestSets[mol->AtomCount].insert(mol->end->previous->nr); 1089 TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr); 1090 1091 // constructing Graph structure 1092 cout << "Generating Subgraph class." << endl; 1093 Graph Subgraphs; 1094 1095 // insert KeySets into Subgraphs 1096 cout << "Inserting KeySets into Subgraph class." << endl; 1097 for (int j=0;j<mol->AtomCount;j++) { 1098 Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.))); 1099 } 1100 cout << "Testing insertion of already present item in Subgraph." << endl; 1101 GraphTestPair test2; 1102 test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.))); 1103 if (test2.second) { 1104 cout << Verbose(1) << "Insertion worked?!" << endl; 1105 } else { 1106 cout << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl; 1107 } 1108 1109 // show graphs 1110 cout << "Showing Subgraph's contents, checking that it's sorted." << endl; 1111 Graph::iterator A = Subgraphs.begin(); 1112 while (A != Subgraphs.end()) { 1113 cout << (*A).second.first << ": "; 1114 KeySet::iterator key = (*A).first.begin(); 1115 comp = -1; 1116 while (key != (*A).first.end()) { 1117 if ((*key) > comp) 1118 cout << (*key) << " "; 1119 else 1120 cout << (*key) << "! "; 1121 comp = (*key); 1122 key++; 1123 } 1124 cout << endl; 1125 A++; 1126 } 1127 delete(mol); 1128 1128 }; 1129 1129 … … 1136 1136 static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, MoleculeListClass *molecules) 1137 1137 { 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1138 char filename[MAXSTRINGSIZE]; 1139 ofstream output; 1140 molecule *mol = new molecule(periode); 1141 1142 // merge all molecules in MoleculeListClass into this molecule 1143 int N = molecules->ListOfMolecules.size(); 1144 int *src = new int(N); 1145 N=0; 1146 for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) 1147 src[N++] = (*ListRunner)->IndexNr; 1148 molecules->SimpleMultiAdd(mol, src, N); 1149 1150 cout << Verbose(0) << "Storing configuration ... " << endl; 1151 // get correct valence orbitals 1152 mol->CalculateOrbitals(*configuration); 1153 configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble; 1154 strcpy(filename, ConfigFileName); 1155 if (ConfigFileName != NULL) { // test the file name 1156 output.open(ConfigFileName, ios::trunc); 1157 } else if (strlen(configuration->configname) != 0) { 1158 strcpy(filename, configuration->configname); 1159 output.open(configuration->configname, ios::trunc); 1160 } else { 1161 strcpy(filename, DEFAULTCONFIG); 1162 output.open(DEFAULTCONFIG, ios::trunc); 1163 } 1164 output.close(); 1165 output.clear(); 1166 cout << Verbose(0) << "Saving of config file "; 1167 if (configuration->Save(filename, periode, mol)) 1168 cout << "successful." << endl; 1169 else 1170 cout << "failed." << endl; 1171 1172 // and save to xyz file 1173 if (ConfigFileName != NULL) { 1174 strcpy(filename, ConfigFileName); 1175 strcat(filename, ".xyz"); 1176 output.open(filename, ios::trunc); 1177 } 1178 if (output == NULL) { 1179 strcpy(filename,"main_pcp_linux"); 1180 strcat(filename, ".xyz"); 1181 output.open(filename, ios::trunc); 1182 } 1183 cout << Verbose(0) << "Saving of XYZ file "; 1184 if (mol->MDSteps <= 1) { 1185 if (mol->OutputXYZ(&output)) 1186 cout << "successful." << endl; 1187 else 1188 cout << "failed." << endl; 1189 } else { 1190 if (mol->OutputTrajectoriesXYZ(&output)) 1191 cout << "successful." << endl; 1192 else 1193 cout << "failed." << endl; 1194 } 1195 output.close(); 1196 output.clear(); 1197 1198 // and save as MPQC configuration 1199 if (ConfigFileName != NULL) 1200 strcpy(filename, ConfigFileName); 1201 if (output == NULL) 1202 strcpy(filename,"main_pcp_linux"); 1203 cout << Verbose(0) << "Saving as mpqc input "; 1204 if (configuration->SaveMPQC(filename, mol)) 1205 cout << "done." << endl; 1206 else 1207 cout << "failed." << endl; 1208 1209 if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) { 1210 cerr << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl; 1211 } 1212 delete(mol); 1213 1213 }; 1214 1214 … … 1225 1225 static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode, config& configuration, char *&ConfigFileName, char *&PathToDatabases) 1226 1226 { 1227 Vector x,y,z,n;// coordinates for absolute point in cell volume1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 default:// no match? Step on1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 mol->AddAtom(first);// add to molecule1405 1406 1407 1408 1409 1410 1411 1412 default:// no match? Don't step on (this is done in next switch's default)1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 MoleculeLeafClass *Subgraphs = NULL;// list of subgraphs from DFS analysis1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 Subgraphs->next->FillBondStructureFromReference((ofstream *)&cout, mol, (FragmentCounter = 0), ListOfLocalAtoms, false);// we want to keep the created ListOfLocalAtoms1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 LinkedCell LCList(mol, atof(argv[argptr]));// \NOTE not twice the radius??1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1227 Vector x,y,z,n; // coordinates for absolute point in cell volume 1228 double *factor; // unit factor if desired 1229 ifstream test; 1230 ofstream output; 1231 string line; 1232 atom *first; 1233 bool SaveFlag = false; 1234 int ExitFlag = 0; 1235 int j; 1236 double volume = 0.; 1237 enum ConfigStatus config_present = absent; 1238 clock_t start,end; 1239 int argptr; 1240 PathToDatabases = LocalPath; 1241 1242 // simply create a new molecule, wherein the config file is loaded and the manipulation takes place 1243 molecule *mol = new molecule(periode); 1244 molecules->insert(mol); 1245 1246 if (argc > 1) { // config file specified as option 1247 // 1. : Parse options that just set variables or print help 1248 argptr = 1; 1249 do { 1250 if (argv[argptr][0] == '-') { 1251 cout << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n"; 1252 argptr++; 1253 switch(argv[argptr-1][1]) { 1254 case 'h': 1255 case 'H': 1256 case '?': 1257 cout << "MoleCuilder suite" << endl << "==================" << endl << endl; 1258 cout << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl; 1259 cout << "or simply " << argv[0] << " without arguments for interactive session." << endl; 1260 cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl; 1261 cout << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl; 1262 cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl; 1263 cout << "\t-B <basis>\tSetting basis to store to MPQC config files." << endl; 1264 cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl; 1265 cout << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl; 1266 cout << "\t-O\tCenter atoms in origin." << endl; 1267 cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl; 1268 cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl; 1269 cout << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl; 1270 cout << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl; 1271 cout << "\t-h/-H/-?\tGive this help screen." << endl; 1272 cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl; 1273 cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl; 1274 cout << "\t-N <radius> <file>\tGet non-convex-envelope." << endl; 1275 cout << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl; 1276 cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl; 1277 cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl; 1278 cout << "\t-r\t\tConvert file from an old pcp syntax." << endl; 1279 cout << "\t-R\t\tRemove all atoms out of sphere around a given one." << endl; 1280 cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl; 1281 cout << "\t-T <file> Store temperatures from the config file in <file>." << endl; 1282 cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl; 1283 cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl; 1284 cout << "\t-v/-V\t\tGives version information." << endl; 1285 cout << "Note: config files must not begin with '-' !" << endl; 1286 delete(mol); 1287 delete(periode); 1288 return (1); 1289 break; 1290 case 'v': 1291 case 'V': 1292 cout << argv[0] << " " << VERSIONSTRING << endl; 1293 cout << "Build your own molecule position set." << endl; 1294 delete(mol); 1295 delete(periode); 1296 return (1); 1297 break; 1298 case 'e': 1299 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1300 cerr << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl; 1301 } else { 1302 cout << "Using " << argv[argptr] << " as elements database." << endl; 1303 PathToDatabases = argv[argptr]; 1304 argptr+=1; 1305 } 1306 break; 1307 case 'n': 1308 cout << "I won't parse trajectories." << endl; 1309 configuration.FastParsing = true; 1310 break; 1311 default: // no match? Step on 1312 argptr++; 1313 break; 1314 } 1315 } else 1316 argptr++; 1317 } while (argptr < argc); 1318 1319 // 2. Parse the element database 1320 if (periode->LoadPeriodentafel(PathToDatabases)) { 1321 cout << Verbose(0) << "Element list loaded successfully." << endl; 1322 //periode->Output((ofstream *)&cout); 1323 } else { 1324 cout << Verbose(0) << "Element list loading failed." << endl; 1325 return 1; 1326 } 1327 // 3. Find config file name and parse if possible 1328 if (argv[1][0] != '-') { 1329 cout << Verbose(0) << "Config file given." << endl; 1330 test.open(argv[1], ios::in); 1331 if (test == NULL) { 1332 //return (1); 1333 output.open(argv[1], ios::out); 1334 if (output == NULL) { 1335 cout << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl; 1336 config_present = absent; 1337 } else { 1338 cout << "Empty configuration file." << endl; 1339 ConfigFileName = argv[1]; 1340 config_present = empty; 1341 output.close(); 1342 } 1343 } else { 1344 test.close(); 1345 ConfigFileName = argv[1]; 1346 cout << Verbose(1) << "Specified config file found, parsing ... "; 1347 switch (configuration.TestSyntax(ConfigFileName, periode, mol)) { 1348 case 1: 1349 cout << "new syntax." << endl; 1350 configuration.Load(ConfigFileName, periode, mol); 1351 config_present = present; 1352 break; 1353 case 0: 1354 cout << "old syntax." << endl; 1355 configuration.LoadOld(ConfigFileName, periode, mol); 1356 config_present = present; 1357 break; 1358 default: 1359 cout << "Unknown syntax or empty, yet present file." << endl; 1360 config_present = empty; 1361 } 1362 } 1363 } else 1364 config_present = absent; 1365 // 4. parse again through options, now for those depending on elements db and config presence 1366 argptr = 1; 1367 do { 1368 cout << "Current Command line argument: " << argv[argptr] << "." << endl; 1369 if (argv[argptr][0] == '-') { 1370 argptr++; 1371 if ((config_present == present) || (config_present == empty)) { 1372 switch(argv[argptr-1][1]) { 1373 case 'p': 1374 ExitFlag = 1; 1375 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1376 ExitFlag = 255; 1377 cerr << "Not enough arguments for parsing: -p <xyz file>" << endl; 1378 } else { 1379 SaveFlag = true; 1380 cout << Verbose(1) << "Parsing xyz file for new atoms." << endl; 1381 if (!mol->AddXYZFile(argv[argptr])) 1382 cout << Verbose(2) << "File not found." << endl; 1383 else { 1384 cout << Verbose(2) << "File found and parsed." << endl; 1385 config_present = present; 1386 } 1387 } 1388 break; 1389 case 'a': 1390 ExitFlag = 1; 1391 if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1]))) { 1392 ExitFlag = 255; 1393 cerr << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl; 1394 } else { 1395 SaveFlag = true; 1396 cout << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), "; 1397 first = new atom; 1398 first->type = periode->FindElement(atoi(argv[argptr])); 1399 if (first->type != NULL) 1400 cout << Verbose(2) << "found element " << first->type->name << endl; 1401 for (int i=NDIM;i--;) 1402 first->x.x[i] = atof(argv[argptr+1+i]); 1403 if (first->type != NULL) { 1404 mol->AddAtom(first); // add to molecule 1405 if ((config_present == empty) && (mol->AtomCount != 0)) 1406 config_present = present; 1407 } else 1408 cerr << Verbose(1) << "Could not find the specified element." << endl; 1409 argptr+=4; 1410 } 1411 break; 1412 default: // no match? Don't step on (this is done in next switch's default) 1413 break; 1414 } 1415 } 1416 if (config_present == present) { 1417 switch(argv[argptr-1][1]) { 1418 case 'B': 1419 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1420 ExitFlag = 255; 1421 cerr << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl; 1422 } else { 1423 configuration.basis = argv[argptr]; 1424 cout << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl; 1425 argptr+=1; 1426 } 1427 break; 1428 case 'D': 1429 ExitFlag = 1; 1430 { 1431 cout << Verbose(1) << "Depth-First-Search Analysis." << endl; 1432 MoleculeLeafClass *Subgraphs = NULL; // list of subgraphs from DFS analysis 1433 int *MinimumRingSize = new int[mol->AtomCount]; 1434 atom ***ListOfLocalAtoms = NULL; 1435 int FragmentCounter = 0; 1436 class StackClass<bond *> *BackEdgeStack = NULL; 1437 class StackClass<bond *> *LocalBackEdgeStack = NULL; 1438 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr]), configuration.GetIsAngstroem()); 1439 mol->CreateListOfBondsPerAtom((ofstream *)&cout); 1440 Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack); 1441 if (Subgraphs != NULL) { 1442 Subgraphs->next->FillBondStructureFromReference((ofstream *)&cout, mol, (FragmentCounter = 0), ListOfLocalAtoms, false); // we want to keep the created ListOfLocalAtoms 1443 while (Subgraphs->next != NULL) { 1444 Subgraphs = Subgraphs->next; 1445 LocalBackEdgeStack = new StackClass<bond *> (Subgraphs->Leaf->BondCount); 1446 Subgraphs->Leaf->PickLocalBackEdges((ofstream *)&cout, ListOfLocalAtoms[FragmentCounter++], BackEdgeStack, LocalBackEdgeStack); 1447 Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize); 1448 delete(LocalBackEdgeStack); 1449 delete(Subgraphs->previous); 1450 } 1451 delete(Subgraphs); 1452 for (int i=0;i<FragmentCounter;i++) 1453 Free((void **)&ListOfLocalAtoms[FragmentCounter], "ParseCommandLineOptions: **ListOfLocalAtoms[]"); 1454 Free((void **)&ListOfLocalAtoms, "ParseCommandLineOptions: ***ListOfLocalAtoms"); 1455 } 1456 delete(BackEdgeStack); 1457 delete[](MinimumRingSize); 1458 } 1459 //argptr+=1; 1460 break; 1461 case 'E': 1462 ExitFlag = 1; 1463 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) { 1464 ExitFlag = 255; 1465 cerr << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl; 1466 } else { 1467 SaveFlag = true; 1468 cout << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl; 1469 first = mol->FindAtom(atoi(argv[argptr])); 1470 first->type = periode->FindElement(atoi(argv[argptr+1])); 1471 argptr+=2; 1472 } 1473 break; 1474 case 'A': 1475 ExitFlag = 1; 1476 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1477 ExitFlag =255; 1478 cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl; 1479 } else { 1480 cout << "Parsing bonds from " << argv[argptr] << "." << endl; 1481 ifstream *input = new ifstream(argv[argptr]); 1482 mol->CreateAdjacencyList2((ofstream *)&cout, input); 1483 input->close(); 1484 argptr+=1; 1485 } 1486 break; 1487 case 'N': 1488 ExitFlag = 1; 1489 if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){ 1490 ExitFlag = 255; 1491 cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl; 1492 } else { 1493 class Tesselation T; 1494 int N = 15; 1495 int number = 100; 1496 string filename(argv[argptr+1]); 1497 filename.append(".csv"); 1498 cout << Verbose(0) << "Evaluating non-convex envelope."; 1499 cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl; 1500 LinkedCell LCList(mol, atof(argv[argptr])); // \NOTE not twice the radius?? 1501 Find_non_convex_border((ofstream *)&cout, mol, &T, &LCList, argv[argptr+1], atof(argv[argptr])); 1502 //FindDistributionOfEllipsoids((ofstream *)&cout, &T, &LCList, N, number, filename.c_str()); 1503 argptr+=2; 1504 } 1505 break; 1506 case 'T': 1507 ExitFlag = 1; 1508 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1509 ExitFlag = 255; 1510 cerr << "Not enough or invalid arguments given for storing tempature: -T <temperature file>" << endl; 1511 } else { 1512 cout << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl; 1513 ofstream *output = new ofstream(argv[argptr], ios::trunc); 1514 if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output)) 1515 cout << Verbose(2) << "File could not be written." << endl; 1516 else 1517 cout << Verbose(2) << "File stored." << endl; 1518 output->close(); 1519 delete(output); 1520 argptr+=1; 1521 } 1522 break; 1523 case 'P': 1524 ExitFlag = 1; 1525 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1526 ExitFlag = 255; 1527 cerr << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl; 1528 } else { 1529 SaveFlag = true; 1530 cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl; 1531 if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem())) 1532 cout << Verbose(2) << "File not found." << endl; 1533 else 1534 cout << Verbose(2) << "File found and parsed." << endl; 1535 argptr+=1; 1536 } 1537 break; 1538 1538 case 'R': 1539 1539 ExitFlag = 1; … … 1561 1561 } 1562 1562 break; 1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596 1597 1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616 1617 1618 1619 1620 1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677 1678 1679 1680 1681 1682 1683 1684 1685 1686 1687 1688 1689 1690 1691 1692 1693 1694 1695 1696 1697 1698 1699 1700 1701 1702 1703 1704 1705 1706 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720 1721 1722 1723 1724 1725 1726 1727 1728 1729 1730 1731 1732 // 1733 // 1734 // 1735 // 1736 // 1737 // 1738 PrepareClustersinWater((ofstream *)&cout, &configuration, mol, volume, density);// if volume == 0, will calculate from ConvexEnvelope1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 1750 1751 1752 1753 1754 1755 1756 1757 mol->CountAtoms((ofstream *)&cout);// recount atoms1758 if (mol->AtomCount != 0) {// if there is more than none1759 count = mol->AtomCount;// is changed becausing of adding, thus has to be stored away beforehand1760 1761 1762 1763 1764 while (first->next != mol->end) {// make a list of all atoms with coordinates and element1765 1766 1767 1768 1769 1770 1771 1772 1773 1774 1775 for (int i=1;i<faktor;i++) {// then add this list with respective translation factor times1776 1777 1778 1779 first->x.CopyVector(vectors[k]);// use coordinate of original atom1780 first->x.AddVector(&x);// translate the coordinates1781 first->type = Elements[k];// insert original element1782 mol->AddAtom(first);// and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)1783 1784 1785 1786 1787 1788 1789 1790 1791 1792 1793 1794 1795 1796 1797 1798 1799 1800 default:// no match? Step on1801 1802 1803 1804 1805 1806 1807 1808 1809 1810 1811 1812 1813 1814 1815 } else {// no arguments, hence scan the elements db1816 1817 1818 1819 1820 1821 1822 1563 case 't': 1564 ExitFlag = 1; 1565 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1566 ExitFlag = 255; 1567 cerr << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl; 1568 } else { 1569 ExitFlag = 1; 1570 SaveFlag = true; 1571 cout << Verbose(1) << "Translating all ions to new origin." << endl; 1572 for (int i=NDIM;i--;) 1573 x.x[i] = atof(argv[argptr+i]); 1574 mol->Translate((const Vector *)&x); 1575 argptr+=3; 1576 } 1577 break; 1578 case 's': 1579 ExitFlag = 1; 1580 if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) { 1581 ExitFlag = 255; 1582 cerr << "Not enough or invalid arguments given for scaling: -s <factor/[factor_x]> [factor_y] [factor_z]" << endl; 1583 } else { 1584 SaveFlag = true; 1585 j = -1; 1586 cout << Verbose(1) << "Scaling all ion positions by factor." << endl; 1587 factor = new double[NDIM]; 1588 factor[0] = atof(argv[argptr]); 1589 if ((argptr < argc) && (IsValidNumber(argv[argptr]))) 1590 argptr++; 1591 factor[1] = atof(argv[argptr]); 1592 if ((argptr < argc) && (IsValidNumber(argv[argptr]))) 1593 argptr++; 1594 factor[2] = atof(argv[argptr]); 1595 mol->Scale(&factor); 1596 for (int i=0;i<NDIM;i++) { 1597 j += i+1; 1598 x.x[i] = atof(argv[NDIM+i]); 1599 mol->cell_size[j]*=factor[i]; 1600 } 1601 delete[](factor); 1602 argptr+=1; 1603 } 1604 break; 1605 case 'b': 1606 ExitFlag = 1; 1607 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1608 ExitFlag = 255; 1609 cerr << "Not enough or invalid arguments given for centering in box: -b <length_x> <length_y> <length_z>" << endl; 1610 } else { 1611 SaveFlag = true; 1612 j = -1; 1613 cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl; 1614 j=-1; 1615 for (int i=0;i<NDIM;i++) { 1616 j += i+1; 1617 x.x[i] = atof(argv[argptr++]); 1618 mol->cell_size[j] += x.x[i]*2.; 1619 } 1620 // center 1621 mol->CenterInBox((ofstream *)&cout, &x); 1622 // update Box of atoms by boundary 1623 mol->SetBoxDimension(&x); 1624 } 1625 break; 1626 case 'c': 1627 ExitFlag = 1; 1628 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1629 ExitFlag = 255; 1630 cerr << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl; 1631 } else { 1632 SaveFlag = true; 1633 j = -1; 1634 cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl; 1635 // make every coordinate positive 1636 mol->CenterEdge((ofstream *)&cout, &x); 1637 // update Box of atoms by boundary 1638 mol->SetBoxDimension(&x); 1639 // translate each coordinate by boundary 1640 j=-1; 1641 for (int i=0;i<NDIM;i++) { 1642 j += i+1; 1643 x.x[i] = atof(argv[argptr++]); 1644 mol->cell_size[j] += x.x[i]*2.; 1645 } 1646 mol->Translate((const Vector *)&x); 1647 } 1648 break; 1649 case 'O': 1650 ExitFlag = 1; 1651 SaveFlag = true; 1652 cout << Verbose(1) << "Centering atoms in origin." << endl; 1653 mol->CenterOrigin((ofstream *)&cout, &x); 1654 mol->SetBoxDimension(&x); 1655 break; 1656 case 'r': 1657 ExitFlag = 1; 1658 SaveFlag = true; 1659 cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl; 1660 break; 1661 case 'F': 1662 case 'f': 1663 ExitFlag = 1; 1664 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) { 1665 ExitFlag = 255; 1666 cerr << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl; 1667 } else { 1668 cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl; 1669 cout << Verbose(0) << "Creating connection matrix..." << endl; 1670 start = clock(); 1671 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr++]), configuration.GetIsAngstroem()); 1672 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl; 1673 if (mol->first->next != mol->last) { 1674 ExitFlag = mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration); 1675 } 1676 end = clock(); 1677 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl; 1678 argptr+=2; 1679 } 1680 break; 1681 case 'm': 1682 ExitFlag = 1; 1683 j = atoi(argv[argptr++]); 1684 if ((j<0) || (j>1)) { 1685 cerr << Verbose(1) << "ERROR: Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl; 1686 j = 0; 1687 } 1688 if (j) { 1689 SaveFlag = true; 1690 cout << Verbose(0) << "Converting to prinicipal axis system." << endl; 1691 } else 1692 cout << Verbose(0) << "Evaluating prinicipal axis." << endl; 1693 mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j); 1694 break; 1695 case 'o': 1696 ExitFlag = 1; 1697 if ((argptr >= argc) || (argv[argptr][0] == '-')){ 1698 ExitFlag = 255; 1699 cerr << "Not enough or invalid arguments given for convex envelope: -o <tecplot output file>" << endl; 1700 } else { 1701 cout << Verbose(0) << "Evaluating volume of the convex envelope."; 1702 cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl; 1703 VolumeOfConvexEnvelope((ofstream *)&cout, argv[argptr], &configuration, NULL, mol); 1704 argptr+=1; 1705 } 1706 break; 1707 case 'U': 1708 ExitFlag = 1; 1709 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) { 1710 ExitFlag = 255; 1711 cerr << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl; 1712 volume = -1; // for case 'u': don't print error again 1713 } else { 1714 volume = atof(argv[argptr++]); 1715 cout << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl; 1716 } 1717 case 'u': 1718 ExitFlag = 1; 1719 if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) { 1720 if (volume != -1) 1721 ExitFlag = 255; 1722 cerr << "Not enough arguments given for suspension: -u <density>" << endl; 1723 } else { 1724 double density; 1725 SaveFlag = true; 1726 cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water."; 1727 density = atof(argv[argptr++]); 1728 if (density < 1.0) { 1729 cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl; 1730 density = 1.3; 1731 } 1732 // for(int i=0;i<NDIM;i++) { 1733 // repetition[i] = atoi(argv[argptr++]); 1734 // if (repetition[i] < 1) 1735 // cerr << Verbose(0) << "ERROR: repetition value must be greater 1!" << endl; 1736 // repetition[i] = 1; 1737 // } 1738 PrepareClustersinWater((ofstream *)&cout, &configuration, mol, volume, density); // if volume == 0, will calculate from ConvexEnvelope 1739 } 1740 break; 1741 case 'd': 1742 ExitFlag = 1; 1743 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1744 ExitFlag = 255; 1745 cerr << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl; 1746 } else { 1747 SaveFlag = true; 1748 for (int axis = 1; axis <= NDIM; axis++) { 1749 int faktor = atoi(argv[argptr++]); 1750 int count; 1751 element ** Elements; 1752 Vector ** vectors; 1753 if (faktor < 1) { 1754 cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl; 1755 faktor = 1; 1756 } 1757 mol->CountAtoms((ofstream *)&cout); // recount atoms 1758 if (mol->AtomCount != 0) { // if there is more than none 1759 count = mol->AtomCount; // is changed becausing of adding, thus has to be stored away beforehand 1760 Elements = new element *[count]; 1761 vectors = new Vector *[count]; 1762 j = 0; 1763 first = mol->start; 1764 while (first->next != mol->end) { // make a list of all atoms with coordinates and element 1765 first = first->next; 1766 Elements[j] = first->type; 1767 vectors[j] = &first->x; 1768 j++; 1769 } 1770 if (count != j) 1771 cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl; 1772 x.Zero(); 1773 y.Zero(); 1774 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 1775 for (int i=1;i<faktor;i++) { // then add this list with respective translation factor times 1776 x.AddVector(&y); // per factor one cell width further 1777 for (int k=count;k--;) { // go through every atom of the original cell 1778 first = new atom(); // create a new body 1779 first->x.CopyVector(vectors[k]); // use coordinate of original atom 1780 first->x.AddVector(&x); // translate the coordinates 1781 first->type = Elements[k]; // insert original element 1782 mol->AddAtom(first); // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...) 1783 } 1784 } 1785 // free memory 1786 delete[](Elements); 1787 delete[](vectors); 1788 // correct cell size 1789 if (axis < 0) { // if sign was negative, we have to translate everything 1790 x.Zero(); 1791 x.AddVector(&y); 1792 x.Scale(-(faktor-1)); 1793 mol->Translate(&x); 1794 } 1795 mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor; 1796 } 1797 } 1798 } 1799 break; 1800 default: // no match? Step on 1801 if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step! 1802 argptr++; 1803 break; 1804 } 1805 } 1806 } else argptr++; 1807 } while (argptr < argc); 1808 if (SaveFlag) 1809 SaveConfig(ConfigFileName, &configuration, periode, molecules); 1810 if ((ExitFlag >= 1)) { 1811 delete(mol); 1812 delete(periode); 1813 return (ExitFlag); 1814 } 1815 } else { // no arguments, hence scan the elements db 1816 if (periode->LoadPeriodentafel(PathToDatabases)) 1817 cout << Verbose(0) << "Element list loaded successfully." << endl; 1818 else 1819 cout << Verbose(0) << "Element list loading failed." << endl; 1820 configuration.RetrieveConfigPathAndName("main_pcp_linux"); 1821 } 1822 return(0); 1823 1823 }; 1824 1824 … … 1827 1827 int main(int argc, char **argv) 1828 1828 { 1829 1830 1831 1832 1833 1834 1835 char choice;// menu choice char1836 Vector x,y,z,n;// coordinates for absolute point in cell volume1837 1838 1839 1840 1841 1842 1843 1844 1845 1846 1847 1848 1849 if (j) return j;// something went wrong1850 1851 1852 1829 periodentafel *periode = new periodentafel; // and a period table of all elements 1830 MoleculeListClass *molecules = new MoleculeListClass; // list of all molecules 1831 molecule *mol = NULL; 1832 config configuration; 1833 double tmp1; 1834 atom *first, *second; 1835 char choice; // menu choice char 1836 Vector x,y,z,n; // coordinates for absolute point in cell volume 1837 bool valid; // flag if input was valid or not 1838 ifstream test; 1839 ofstream output; 1840 string line; 1841 char *ConfigFileName = NULL; 1842 char *ElementsFileName = NULL; 1843 int Z; 1844 int j, axis, count, faktor; 1845 1846 // =========================== PARSE COMMAND LINE OPTIONS ==================================== 1847 j = ParseCommandLineOptions(argc, argv, molecules, periode, configuration, ConfigFileName, ElementsFileName); 1848 if (j == 1) return 0; // just for -v and -h options 1849 if (j) return j; // something went wrong 1850 1851 // General stuff 1852 if (molecules->ListOfMolecules.size() == 0) { 1853 1853 mol = new molecule(periode); 1854 1854 if (mol->cell_size[0] == 0.) { … … 1860 1860 } 1861 1861 molecules->insert(mol); 1862 1863 1864 1865 1866 1867 1868 1869 1870 1871 1872 1862 } 1863 1864 // =========================== START INTERACTIVE SESSION ==================================== 1865 1866 // now the main construction loop 1867 cout << Verbose(0) << endl << "Now comes the real construction..." << endl; 1868 do { 1869 cout << Verbose(0) << endl << endl; 1870 cout << Verbose(0) << "============Molecule list=======================" << endl; 1871 molecules->Enumerate((ofstream *)&cout); 1872 cout << Verbose(0) << "============Menu===============================" << endl; 1873 1873 cout << Verbose(0) << "a - set molecule (in)active" << endl; 1874 1874 cout << Verbose(0) << "e - edit new molecules" << endl; … … 1886 1886 cin >> choice; 1887 1887 1888 1889 1888 switch (choice) { 1889 case 'a': // (in)activate molecule 1890 1890 { 1891 1891 cout << "Enter index of molecule: "; … … 1897 1897 (*ListRunner)->ActiveFlag = !(*ListRunner)->ActiveFlag; 1898 1898 } 1899 1900 1901 1902 1903 1899 break; 1900 1901 case 'c': // edit each field of the configuration 1902 configuration.Edit(); 1903 break; 1904 1904 1905 1905 case 'e': // create molecule … … 1919 1919 break; 1920 1920 1921 1922 1923 1924 1925 1926 1927 1928 1929 1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 1945 1921 case 'q': // quit 1922 break; 1923 1924 case 's': // save to config file 1925 SaveConfig(ConfigFileName, &configuration, periode, molecules); 1926 break; 1927 1928 case 'T': 1929 testroutine(molecules); 1930 break; 1931 1932 default: 1933 break; 1934 }; 1935 } while (choice != 'q'); 1936 1937 // save element data base 1938 if (periode->StorePeriodentafel(ElementsFileName)) //ElementsFileName 1939 cout << Verbose(0) << "Saving of elements.db successful." << endl; 1940 else 1941 cout << Verbose(0) << "Saving of elements.db failed." << endl; 1942 1943 delete(molecules); 1944 delete(periode); 1945 return (0); 1946 1946 } 1947 1947
Note:
See TracChangeset
for help on using the changeset viewer.
