Changeset 042f82 for src/builder.cpp
- Timestamp:
- Jul 23, 2009, 2:21:57 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, Candidate_v1.7.0, 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:
- 36ec71
- Parents:
- 205ccd
- File:
-
- 1 edited
-
src/builder.cpp (modified) (22 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/builder.cpp
r205ccd r042f82 15 15 * \section about About the Program 16 16 * 17 * Molecuilder is a short program, written in C++, that enables the construction of a coordinate set for the18 * atoms making up an molecule by the successive statement of binding angles and distances and referencing to19 * already constructed atoms.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 * A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello22 * molecular dynamics implementation.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 * Installation should without problems succeed as follows:27 * -# ./configure (or: mkdir build;mkdir run;cd build; ../configure --bindir=../run)28 * -# make29 * -# make install26 * 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 * Further useful commands are32 * -# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n33 * -# make doxygen-doc: Creates these html pages out of the documented source31 * 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 * The program can be executed by running: ./molecuilder37 * The program can be executed by running: ./molecuilder 38 38 * 39 * Note, that it uses a database, called "elements.db", in the executable's directory. If the file is not found,40 * it is created and any given data on elements of the periodic table will be stored therein and re-used on41 * later re-execution.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 * For the special configuration file format, see the documentation of pcp.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 atom *first, *second, *third, *fourth;66 Vector **atoms;67 Vector x,y,z,n;// coordinates for absolute point in cell volume68 double a,b,c;69 char choice;// menu choice char70 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) {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 case 'a': // absolute coordinates of atom88 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(); // give type93 mol->AddAtom(first); // add to molecule94 break;95 96 case 'b': // relative coordinates of atom wrt to reference point92 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(); // give type109 mol->AddAtom(first); // add to molecule110 break;111 112 case 'c': // relative coordinates of atom wrt to already placed atom108 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(); // give type125 mol->AddAtom(first); // add to molecule124 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(); // give type227 mol->AddAtom(first); // add to molecule228 break;229 230 case 'e': // least square distance position to a set of atoms226 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(); // give type251 mol->AddAtom(first); // add to molecule250 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 break;257 };256 break; 257 }; 258 258 }; 259 259 … … 263 263 static void CenterAtoms(molecule *mol) 264 264 { 265 Vector x, y, helper;266 char choice;// menu choice char267 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 positive297 mol->Translate(&y); // translate by boundary298 helper.CopyVector(&y);299 helper.Scale(2.);300 helper.AddVector(&x);301 mol->SetBoxDimension(&helper);// update Box of atoms by boundary302 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 // center310 mol->CenterInBox((ofstream *)&cout);311 // update Box of atoms by boundary312 mol->SetBoxDimension(&x);313 break;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); 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 atom *first, *second, *third;324 Vector x,n;325 char choice;// menu choice char326 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 plane340 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 plane347 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 plane352 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);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 atom *first, *second, *third;392 Vector n;393 char choice;// menu choice char394 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 plane407 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 plane414 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 plane419 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);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;439 int axis;440 double tmp1, tmp2;441 char choice;// menu choice char442 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 else458 cout << Verbose(1) << "Atom not found." << endl;459 break;460 case 'b':461 second = mol->AskAtom("Enter number of atom as reference point: ");462 cout << Verbose(0) << "Enter radius: ";463 cin >> tmp1;464 first = mol->start;465 while(first->next != mol->end) {466 first = first->next;467 if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...468 mol->RemoveAtom(first);469 }470 break;471 case 'c':472 cout << Verbose(0) << "Which axis is it: ";473 cin >> axis;474 cout << Verbose(0) << "Left inward boundary: ";475 cin >> tmp1;476 cout << Verbose(0) << "Right inward boundary: ";477 cin >> tmp2;478 first = mol->start;479 while(first->next != mol->end) {480 first = first->next;481 if ((first->x.x[axis] > tmp2) || (first->x.x[axis] < tmp1)) // out of boundary ...482 mol->RemoveAtom(first);483 }484 break;485 };486 //mol->Output((ofstream *)&cout);487 choice = 'r';438 atom *first, *second; 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 second = mol->AskAtom("Enter number of atom as reference point: "); 462 cout << Verbose(0) << "Enter radius: "; 463 cin >> tmp1; 464 first = mol->start; 465 while(first->next != mol->end) { 466 first = first->next; 467 if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ... 468 mol->RemoveAtom(first); 469 } 470 break; 471 case 'c': 472 cout << Verbose(0) << "Which axis is it: "; 473 cin >> axis; 474 cout << Verbose(0) << "Left inward boundary: "; 475 cin >> tmp1; 476 cout << Verbose(0) << "Right inward boundary: "; 477 cin >> tmp2; 478 first = mol->start; 479 while(first->next != mol->end) { 480 first = first->next; 481 if ((first->x.x[axis] > tmp2) || (first->x.x[axis] < tmp1)) // out of boundary ... 482 mol->RemoveAtom(first); 483 } 484 break; 485 }; 486 //mol->Output((ofstream *)&cout); 487 choice = 'r'; 488 488 }; 489 489 … … 494 494 static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration) 495 495 { 496 atom *first, *second, *third;497 Vector x,y;498 double min[256], tmp1, tmp2, tmp3;499 int Z;500 char choice;// menu choice char501 502 cout << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;503 cout << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;504 cout << Verbose(0) << " b - calculate bond length between two atoms" << endl;505 cout << Verbose(0) << " c - calculate bond angle" << endl;506 cout << Verbose(0) << " d - calculate principal axis of the system" << endl;507 cout << Verbose(0) << " e - calculate volume of the convex envelope" << endl;508 cout << Verbose(0) << " f - calculate temperature from current velocity" << endl;509 cout << Verbose(0) << " g - output all temperatures per step from velocities" << endl;510 cout << Verbose(0) << "all else - go back" << endl;511 cout << Verbose(0) << "===============================================" << endl;512 cout << Verbose(0) << "INPUT: ";513 cin >> choice;514 515 switch(choice) {516 default:517 cout << Verbose(1) << "Not a valid choice." << endl;518 break;519 case 'a':520 first = mol->AskAtom("Enter first atom: ");521 for (int i=MAX_ELEMENTS;i--;)522 min[i] = 0.;523 524 second = mol->start;525 while ((second->next != mol->end)) {526 second = second->next; // advance527 Z = second->type->Z;528 tmp1 = 0.;529 if (first != second) {530 x.CopyVector((const Vector *)&first->x);531 x.SubtractVector((const Vector *)&second->x);532 tmp1 = x.Norm();533 }534 if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;535 //cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;536 }537 for (int i=MAX_ELEMENTS;i--;)538 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;539 break;540 541 case 'b':542 first = mol->AskAtom("Enter first atom: ");543 second = mol->AskAtom("Enter second atom: ");544 for (int i=NDIM;i--;)545 min[i] = 0.;546 x.CopyVector((const Vector *)&first->x);547 x.SubtractVector((const Vector *)&second->x);548 tmp1 = x.Norm();549 cout << Verbose(1) << "Distance vector is ";550 x.Output((ofstream *)&cout);551 cout << "." << endl << "Norm of distance is " << tmp1 << "." << endl;552 break;553 554 case 'c':555 cout << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl;556 first = mol->AskAtom("Enter first atom: ");557 second = mol->AskAtom("Enter central atom: ");558 third= mol->AskAtom("Enter last atom: ");559 tmp1 = tmp2 = tmp3 = 0.;560 x.CopyVector((const Vector *)&first->x);561 x.SubtractVector((const Vector *)&second->x);562 y.CopyVector((const Vector *)&third->x);563 y.SubtractVector((const Vector *)&second->x);564 cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";565 cout << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;566 break;567 case 'd':568 cout << Verbose(0) << "Evaluating prinicipal axis." << endl;569 cout << Verbose(0) << "Shall we rotate? [0/1]: ";570 cin >> Z;571 if ((Z >=0) && (Z <=1))572 mol->PrincipalAxisSystem((ofstream *)&cout, (bool)Z);573 else574 mol->PrincipalAxisSystem((ofstream *)&cout, false);575 break;576 case 'e':577 cout << Verbose(0) << "Evaluating volume of the convex envelope.";578 VolumeOfConvexEnvelope((ofstream *)&cout, NULL, configuration, NULL, mol);579 break;580 case 'f':581 mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps, (ofstream *)&cout);582 break;583 case 'g':584 {585 char filename[255];586 cout << "Please enter filename: " << endl;587 cin >> filename;588 cout << Verbose(1) << "Storing temperatures in " << filename << "." << endl;589 ofstream *output = new ofstream(filename, ios::trunc);590 if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output))591 cout << Verbose(2) << "File could not be written." << endl;592 else593 cout << Verbose(2) << "File stored." << endl;594 output->close();595 delete(output);596 }597 break;598 }496 atom *first, *second, *third; 497 Vector x,y; 498 double min[256], tmp1, tmp2, tmp3; 499 int Z; 500 char choice; // menu choice char 501 502 cout << Verbose(0) << "===========MEASURE ATOMS=========================" << endl; 503 cout << Verbose(0) << " a - calculate bond length between one atom and all others" << endl; 504 cout << Verbose(0) << " b - calculate bond length between two atoms" << endl; 505 cout << Verbose(0) << " c - calculate bond angle" << endl; 506 cout << Verbose(0) << " d - calculate principal axis of the system" << endl; 507 cout << Verbose(0) << " e - calculate volume of the convex envelope" << endl; 508 cout << Verbose(0) << " f - calculate temperature from current velocity" << endl; 509 cout << Verbose(0) << " g - output all temperatures per step from velocities" << endl; 510 cout << Verbose(0) << "all else - go back" << endl; 511 cout << Verbose(0) << "===============================================" << endl; 512 cout << Verbose(0) << "INPUT: "; 513 cin >> choice; 514 515 switch(choice) { 516 default: 517 cout << Verbose(1) << "Not a valid choice." << endl; 518 break; 519 case 'a': 520 first = mol->AskAtom("Enter first atom: "); 521 for (int i=MAX_ELEMENTS;i--;) 522 min[i] = 0.; 523 524 second = mol->start; 525 while ((second->next != mol->end)) { 526 second = second->next; // advance 527 Z = second->type->Z; 528 tmp1 = 0.; 529 if (first != second) { 530 x.CopyVector((const Vector *)&first->x); 531 x.SubtractVector((const Vector *)&second->x); 532 tmp1 = x.Norm(); 533 } 534 if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1; 535 //cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl; 536 } 537 for (int i=MAX_ELEMENTS;i--;) 538 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; 539 break; 540 541 case 'b': 542 first = mol->AskAtom("Enter first atom: "); 543 second = mol->AskAtom("Enter second atom: "); 544 for (int i=NDIM;i--;) 545 min[i] = 0.; 546 x.CopyVector((const Vector *)&first->x); 547 x.SubtractVector((const Vector *)&second->x); 548 tmp1 = x.Norm(); 549 cout << Verbose(1) << "Distance vector is "; 550 x.Output((ofstream *)&cout); 551 cout << "." << endl << "Norm of distance is " << tmp1 << "." << endl; 552 break; 553 554 case 'c': 555 cout << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl; 556 first = mol->AskAtom("Enter first atom: "); 557 second = mol->AskAtom("Enter central atom: "); 558 third = mol->AskAtom("Enter last atom: "); 559 tmp1 = tmp2 = tmp3 = 0.; 560 x.CopyVector((const Vector *)&first->x); 561 x.SubtractVector((const Vector *)&second->x); 562 y.CopyVector((const Vector *)&third->x); 563 y.SubtractVector((const Vector *)&second->x); 564 cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": "; 565 cout << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl; 566 break; 567 case 'd': 568 cout << Verbose(0) << "Evaluating prinicipal axis." << endl; 569 cout << Verbose(0) << "Shall we rotate? [0/1]: "; 570 cin >> Z; 571 if ((Z >=0) && (Z <=1)) 572 mol->PrincipalAxisSystem((ofstream *)&cout, (bool)Z); 573 else 574 mol->PrincipalAxisSystem((ofstream *)&cout, false); 575 break; 576 case 'e': 577 cout << Verbose(0) << "Evaluating volume of the convex envelope."; 578 VolumeOfConvexEnvelope((ofstream *)&cout, NULL, configuration, NULL, mol); 579 break; 580 case 'f': 581 mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps, (ofstream *)&cout); 582 break; 583 case 'g': 584 { 585 char filename[255]; 586 cout << "Please enter filename: " << endl; 587 cin >> filename; 588 cout << Verbose(1) << "Storing temperatures in " << filename << "." << endl; 589 ofstream *output = new ofstream(filename, ios::trunc); 590 if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output)) 591 cout << Verbose(2) << "File could not be written." << endl; 592 else 593 cout << Verbose(2) << "File stored." << endl; 594 output->close(); 595 delete(output); 596 } 597 break; 598 } 599 599 }; 600 600 … … 605 605 static void FragmentAtoms(molecule *mol, config *configuration) 606 606 { 607 int Order1;608 clock_t start, end;609 610 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;611 cout << Verbose(0) << "What's the desired bond order: ";612 cin >> Order1;613 if (mol->first->next != mol->last) {// there are bonds614 start = clock();615 mol->FragmentMolecule((ofstream *)&cout, Order1, configuration);616 end = clock();617 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;618 } else619 cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl;607 int Order1; 608 clock_t start, end; 609 610 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl; 611 cout << Verbose(0) << "What's the desired bond order: "; 612 cin >> Order1; 613 if (mol->first->next != mol->last) { // there are bonds 614 start = clock(); 615 mol->FragmentMolecule((ofstream *)&cout, Order1, configuration); 616 end = clock(); 617 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl; 618 } else 619 cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl; 620 620 }; 621 621 … … 1038 1038 static void testroutine(MoleculeListClass *molecules) 1039 1039 { 1040 // the current test routine checks the functionality of the KeySet&Graph concept:1041 // We want to have a multiindex (the KeySet) describing a unique subgraph1040 // the current test routine checks the functionality of the KeySet&Graph concept: 1041 // We want to have a multiindex (the KeySet) describing a unique subgraph 1042 1042 int i, comp, counter=0; 1043 1043 … … 1052 1052 atom *Walker = mol->start; 1053 1053 1054 // generate some KeySets1055 cout << "Generating KeySets." << endl;1056 KeySet TestSets[mol->AtomCount+1];1057 i=1;1058 while (Walker->next != mol->end) {1059 Walker = Walker->next;1060 for (int j=0;j<i;j++) {1061 TestSets[j].insert(Walker->nr);1062 }1063 i++;1064 }1065 cout << "Testing insertion of already present item in KeySets." << endl;1066 KeySetTestPair test;1067 test = TestSets[mol->AtomCount-1].insert(Walker->nr);1068 if (test.second) {1069 cout << Verbose(1) << "Insertion worked?!" << endl;1070 } else {1071 cout << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;1072 }1073 TestSets[mol->AtomCount].insert(mol->end->previous->nr);1074 TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);1075 1076 // constructing Graph structure1077 cout << "Generating Subgraph class." << endl;1078 Graph Subgraphs;1079 1080 // insert KeySets into Subgraphs1081 cout << "Inserting KeySets into Subgraph class." << endl;1082 for (int j=0;j<mol->AtomCount;j++) {1083 Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));1084 }1085 cout << "Testing insertion of already present item in Subgraph." << endl;1086 GraphTestPair test2;1087 test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.)));1088 if (test2.second) {1089 cout << Verbose(1) << "Insertion worked?!" << endl;1090 } else {1091 cout << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;1092 }1093 1094 // show graphs1095 cout << "Showing Subgraph's contents, checking that it's sorted." << endl;1096 Graph::iterator A = Subgraphs.begin();1097 while (A !=Subgraphs.end()) {1098 cout << (*A).second.first << ": ";1099 KeySet::iterator key = (*A).first.begin();1100 comp = -1;1101 while (key != (*A).first.end()) {1102 if ((*key) > comp)1103 cout << (*key) << " ";1104 else1105 cout << (*key) << "! ";1106 comp = (*key);1107 key++;1108 }1109 cout << endl;1110 A++;1111 }1112 delete(mol);1054 // generate some KeySets 1055 cout << "Generating KeySets." << endl; 1056 KeySet TestSets[mol->AtomCount+1]; 1057 i=1; 1058 while (Walker->next != mol->end) { 1059 Walker = Walker->next; 1060 for (int j=0;j<i;j++) { 1061 TestSets[j].insert(Walker->nr); 1062 } 1063 i++; 1064 } 1065 cout << "Testing insertion of already present item in KeySets." << endl; 1066 KeySetTestPair test; 1067 test = TestSets[mol->AtomCount-1].insert(Walker->nr); 1068 if (test.second) { 1069 cout << Verbose(1) << "Insertion worked?!" << endl; 1070 } else { 1071 cout << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl; 1072 } 1073 TestSets[mol->AtomCount].insert(mol->end->previous->nr); 1074 TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr); 1075 1076 // constructing Graph structure 1077 cout << "Generating Subgraph class." << endl; 1078 Graph Subgraphs; 1079 1080 // insert KeySets into Subgraphs 1081 cout << "Inserting KeySets into Subgraph class." << endl; 1082 for (int j=0;j<mol->AtomCount;j++) { 1083 Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.))); 1084 } 1085 cout << "Testing insertion of already present item in Subgraph." << endl; 1086 GraphTestPair test2; 1087 test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.))); 1088 if (test2.second) { 1089 cout << Verbose(1) << "Insertion worked?!" << endl; 1090 } else { 1091 cout << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl; 1092 } 1093 1094 // show graphs 1095 cout << "Showing Subgraph's contents, checking that it's sorted." << endl; 1096 Graph::iterator A = Subgraphs.begin(); 1097 while (A != Subgraphs.end()) { 1098 cout << (*A).second.first << ": "; 1099 KeySet::iterator key = (*A).first.begin(); 1100 comp = -1; 1101 while (key != (*A).first.end()) { 1102 if ((*key) > comp) 1103 cout << (*key) << " "; 1104 else 1105 cout << (*key) << "! "; 1106 comp = (*key); 1107 key++; 1108 } 1109 cout << endl; 1110 A++; 1111 } 1112 delete(mol); 1113 1113 }; 1114 1114 … … 1121 1121 static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, MoleculeListClass *molecules) 1122 1122 { 1123 char filename[MAXSTRINGSIZE];1124 ofstream output;1125 molecule *mol = new molecule(periode);1126 1127 // merge all molecules in MoleculeListClass into this molecule1128 int N = molecules->ListOfMolecules.size();1129 int *src = new int(N);1130 N=0;1131 for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)1132 src[N++] = (*ListRunner)->IndexNr;1133 molecules->SimpleMultiAdd(mol, src, N);1134 1135 cout << Verbose(0) << "Storing configuration ... " << endl;1136 // get correct valence orbitals1137 mol->CalculateOrbitals(*configuration);1138 configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble;1139 strcpy(filename, ConfigFileName);1140 if (ConfigFileName != NULL) { // test the file name1141 output.open(ConfigFileName, ios::trunc);1142 } else if (strlen(configuration->configname) != 0) {1143 strcpy(filename, configuration->configname);1144 output.open(configuration->configname, ios::trunc);1145 } else {1146 strcpy(filename, DEFAULTCONFIG);1147 output.open(DEFAULTCONFIG, ios::trunc);1148 }1149 output.close();1150 output.clear();1151 cout << Verbose(0) << "Saving of config file ";1152 if (configuration->Save(filename, periode, mol))1153 cout << "successful." << endl;1154 else1155 cout << "failed." << endl;1156 1157 // and save to xyz file1158 if (ConfigFileName != NULL) {1159 strcpy(filename, ConfigFileName);1160 strcat(filename, ".xyz");1161 output.open(filename, ios::trunc);1162 }1163 if (output == NULL) {1164 strcpy(filename,"main_pcp_linux");1165 strcat(filename, ".xyz");1166 output.open(filename, ios::trunc);1167 }1168 cout << Verbose(0) << "Saving of XYZ file ";1169 if (mol->MDSteps <= 1) {1170 if (mol->OutputXYZ(&output))1171 cout << "successful." << endl;1172 else1173 cout << "failed." << endl;1174 } else {1175 if (mol->OutputTrajectoriesXYZ(&output))1176 cout << "successful." << endl;1177 else1178 cout << "failed." << endl;1179 }1180 output.close();1181 output.clear();1182 1183 // and save as MPQC configuration1184 if (ConfigFileName != NULL)1185 strcpy(filename, ConfigFileName);1186 if (output == NULL)1187 strcpy(filename,"main_pcp_linux");1188 cout << Verbose(0) << "Saving as mpqc input ";1189 if (configuration->SaveMPQC(filename, mol))1190 cout << "done." << endl;1191 else1192 cout << "failed." << endl;1193 1194 if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {1195 cerr << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl;1196 }1197 delete(mol);1123 char filename[MAXSTRINGSIZE]; 1124 ofstream output; 1125 molecule *mol = new molecule(periode); 1126 1127 // merge all molecules in MoleculeListClass into this molecule 1128 int N = molecules->ListOfMolecules.size(); 1129 int *src = new int(N); 1130 N=0; 1131 for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) 1132 src[N++] = (*ListRunner)->IndexNr; 1133 molecules->SimpleMultiAdd(mol, src, N); 1134 1135 cout << Verbose(0) << "Storing configuration ... " << endl; 1136 // get correct valence orbitals 1137 mol->CalculateOrbitals(*configuration); 1138 configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble; 1139 strcpy(filename, ConfigFileName); 1140 if (ConfigFileName != NULL) { // test the file name 1141 output.open(ConfigFileName, ios::trunc); 1142 } else if (strlen(configuration->configname) != 0) { 1143 strcpy(filename, configuration->configname); 1144 output.open(configuration->configname, ios::trunc); 1145 } else { 1146 strcpy(filename, DEFAULTCONFIG); 1147 output.open(DEFAULTCONFIG, ios::trunc); 1148 } 1149 output.close(); 1150 output.clear(); 1151 cout << Verbose(0) << "Saving of config file "; 1152 if (configuration->Save(filename, periode, mol)) 1153 cout << "successful." << endl; 1154 else 1155 cout << "failed." << endl; 1156 1157 // and save to xyz file 1158 if (ConfigFileName != NULL) { 1159 strcpy(filename, ConfigFileName); 1160 strcat(filename, ".xyz"); 1161 output.open(filename, ios::trunc); 1162 } 1163 if (output == NULL) { 1164 strcpy(filename,"main_pcp_linux"); 1165 strcat(filename, ".xyz"); 1166 output.open(filename, ios::trunc); 1167 } 1168 cout << Verbose(0) << "Saving of XYZ file "; 1169 if (mol->MDSteps <= 1) { 1170 if (mol->OutputXYZ(&output)) 1171 cout << "successful." << endl; 1172 else 1173 cout << "failed." << endl; 1174 } else { 1175 if (mol->OutputTrajectoriesXYZ(&output)) 1176 cout << "successful." << endl; 1177 else 1178 cout << "failed." << endl; 1179 } 1180 output.close(); 1181 output.clear(); 1182 1183 // and save as MPQC configuration 1184 if (ConfigFileName != NULL) 1185 strcpy(filename, ConfigFileName); 1186 if (output == NULL) 1187 strcpy(filename,"main_pcp_linux"); 1188 cout << Verbose(0) << "Saving as mpqc input "; 1189 if (configuration->SaveMPQC(filename, mol)) 1190 cout << "done." << endl; 1191 else 1192 cout << "failed." << endl; 1193 1194 if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) { 1195 cerr << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl; 1196 } 1197 delete(mol); 1198 1198 }; 1199 1199 … … 1210 1210 static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode, config& configuration, char *&ConfigFileName) 1211 1211 { 1212 Vector x,y,z,n;// coordinates for absolute point in cell volume1213 double *factor; // unit factor if desired1214 ifstream test;1215 ofstream output;1216 string line;1217 atom *first;1218 bool SaveFlag = false;1219 int ExitFlag = 0;1220 int j;1221 double volume = 0.;1222 enum ConfigStatus config_present = absent;1223 clock_t start,end;1224 int argptr;1212 Vector x,y,z,n; // coordinates for absolute point in cell volume 1213 double *factor; // unit factor if desired 1214 ifstream test; 1215 ofstream output; 1216 string line; 1217 atom *first; 1218 bool SaveFlag = false; 1219 int ExitFlag = 0; 1220 int j; 1221 double volume = 0.; 1222 enum ConfigStatus config_present = absent; 1223 clock_t start,end; 1224 int argptr; 1225 1225 strncpy(configuration.databasepath, LocalPath, MAXSTRINGSIZE-1); 1226 1226 1227 // simply create a new molecule, wherein the config file is loaded and the manipulation takes place1228 molecule *mol = new molecule(periode);1229 molecules->insert(mol);1230 1231 if (argc > 1) { // config file specified as option1232 // 1. : Parse options that just set variables or print help1233 argptr = 1;1234 do {1235 if (argv[argptr][0] == '-') {1236 cout << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n";1237 argptr++;1238 switch(argv[argptr-1][1]) {1239 case 'h':1240 case 'H':1241 case '?':1242 cout << "MoleCuilder suite" << endl << "==================" << endl << endl;1243 cout << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl;1244 cout << "or simply " << argv[0] << " without arguments for interactive session." << endl;1245 cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;1246 cout << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl;1247 cout << "\t-b xx xy xz yy yz zz\tCenter atoms in domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;1248 cout << "\t-B <basis>\tSetting basis to store to MPQC config files." << endl;1249 cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;1250 cout << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;1251 cout << "\t-O\tCenter atoms in origin." << endl;1252 cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;1253 cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;1254 cout << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;1255 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;1256 cout << "\t-h/-H/-?\tGive this help screen." << endl;1257 cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;1258 cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;1259 cout << "\t-N\tGet non-convex-envelope." << endl;1260 cout << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;1261 cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;1262 cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;1263 cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;1264 cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;1227 // simply create a new molecule, wherein the config file is loaded and the manipulation takes place 1228 molecule *mol = new molecule(periode); 1229 molecules->insert(mol); 1230 1231 if (argc > 1) { // config file specified as option 1232 // 1. : Parse options that just set variables or print help 1233 argptr = 1; 1234 do { 1235 if (argv[argptr][0] == '-') { 1236 cout << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n"; 1237 argptr++; 1238 switch(argv[argptr-1][1]) { 1239 case 'h': 1240 case 'H': 1241 case '?': 1242 cout << "MoleCuilder suite" << endl << "==================" << endl << endl; 1243 cout << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl; 1244 cout << "or simply " << argv[0] << " without arguments for interactive session." << endl; 1245 cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl; 1246 cout << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl; 1247 cout << "\t-b xx xy xz yy yz zz\tCenter atoms in domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl; 1248 cout << "\t-B <basis>\tSetting basis to store to MPQC config files." << endl; 1249 cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl; 1250 cout << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl; 1251 cout << "\t-O\tCenter atoms in origin." << endl; 1252 cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl; 1253 cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl; 1254 cout << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl; 1255 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; 1256 cout << "\t-h/-H/-?\tGive this help screen." << endl; 1257 cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl; 1258 cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl; 1259 cout << "\t-N\tGet non-convex-envelope." << endl; 1260 cout << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl; 1261 cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl; 1262 cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl; 1263 cout << "\t-r\t\tConvert file from an old pcp syntax." << endl; 1264 cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl; 1265 1265 cout << "\t-T x1 x2 x3\tTranslate periodically all atoms by this vector (x1,x2,x3)." << endl; 1266 cout << "\t-S <file> Store temperatures from the config file in <file>." << endl;1267 cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;1268 cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl;1269 cout << "\t-v/-V\t\tGives version information." << endl;1270 cout << "Note: config files must not begin with '-' !" << endl;1271 delete(mol);1272 delete(periode);1273 return (1);1274 break;1275 case 'v':1276 case 'V':1277 cout << argv[0] << " " << VERSIONSTRING << endl;1278 cout << "Build your own molecule position set." << endl;1279 delete(mol);1280 delete(periode);1281 return (1);1282 break;1283 case 'e':1284 if ((argptr >= argc) || (argv[argptr][0] == '-')) {1285 cerr << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl;1286 } else {1287 cout << "Using " << argv[argptr] << " as elements database." << endl;1288 strncpy (configuration.databasepath, argv[argptr], MAXSTRINGSIZE-1);1289 argptr+=1;1290 }1291 break;1292 case 'n':1293 cout << "I won't parse trajectories." << endl;1294 configuration.FastParsing = true;1295 break;1296 default:// no match? Step on1297 argptr++;1298 break;1299 }1300 } else1301 argptr++;1302 } while (argptr < argc);1303 1304 // 2. Parse the element database1305 if (periode->LoadPeriodentafel(configuration.databasepath)) {1306 cout << Verbose(0) << "Element list loaded successfully." << endl;1307 //periode->Output((ofstream *)&cout);1308 } else {1309 cout << Verbose(0) << "Element list loading failed." << endl;1310 return 1;1311 }1312 // 3. Find config file name and parse if possible1313 if (argv[1][0] != '-') {1314 cout << Verbose(0) << "Config file given." << endl;1315 test.open(argv[1], ios::in);1316 if (test == NULL) {1317 //return (1);1318 output.open(argv[1], ios::out);1319 if (output == NULL) {1320 cout << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl;1321 config_present = absent;1322 } else {1323 cout << "Empty configuration file." << endl;1324 ConfigFileName = argv[1];1325 config_present = empty;1326 output.close();1327 }1328 } else {1329 test.close();1330 ConfigFileName = argv[1];1331 cout << Verbose(1) << "Specified config file found, parsing ... ";1332 switch (configuration.TestSyntax(ConfigFileName, periode, mol)) {1333 case 1:1334 cout << "new syntax." << endl;1335 configuration.Load(ConfigFileName, periode, mol);1336 config_present = present;1337 break;1338 case 0:1339 cout << "old syntax." << endl;1340 configuration.LoadOld(ConfigFileName, periode, mol);1341 config_present = present;1342 break;1343 default:1344 cout << "Unknown syntax or empty, yet present file." << endl;1345 config_present = empty;1346 }1347 }1348 } else1349 config_present = absent;1350 // 4. parse again through options, now for those depending on elements db and config presence1351 argptr = 1;1352 do {1353 cout << "Current Command line argument: " << argv[argptr] << "." << endl;1354 if (argv[argptr][0] == '-') {1355 argptr++;1356 if ((config_present == present) || (config_present == empty)) {1357 switch(argv[argptr-1][1]) {1358 case 'p':1359 ExitFlag = 1;1360 if ((argptr >= argc) || (argv[argptr][0] == '-')) {1361 ExitFlag = 255;1362 cerr << "Not enough arguments for parsing: -p <xyz file>" << endl;1363 } else {1364 SaveFlag = true;1365 cout << Verbose(1) << "Parsing xyz file for new atoms." << endl;1366 if (!mol->AddXYZFile(argv[argptr]))1367 cout << Verbose(2) << "File not found." << endl;1368 else {1369 cout << Verbose(2) << "File found and parsed." << endl;1370 config_present = present;1371 }1372 }1373 break;1374 case 'a':1375 ExitFlag = 1;1376 if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1]))) {1377 ExitFlag = 255;1378 cerr << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl;1379 } else {1380 SaveFlag = true;1381 cout << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";1382 first = new atom;1383 first->type = periode->FindElement(atoi(argv[argptr]));1384 if (first->type != NULL)1385 cout << Verbose(2) << "found element " << first->type->name << endl;1386 for (int i=NDIM;i--;)1387 first->x.x[i] = atof(argv[argptr+1+i]);1388 if (first->type != NULL) {1389 mol->AddAtom(first);// add to molecule1390 if ((config_present == empty) && (mol->AtomCount != 0))1391 config_present = present;1392 } else1393 cerr << Verbose(1) << "Could not find the specified element." << endl;1394 argptr+=4;1395 }1396 break;1397 default:// no match? Don't step on (this is done in next switch's default)1398 break;1399 }1400 }1401 if (config_present == present) {1402 switch(argv[argptr-1][1]) {1403 case 'B':1404 if ((argptr >= argc) || (argv[argptr][0] == '-')) {1405 ExitFlag = 255;1406 cerr << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl;1407 } else {1408 configuration.basis = argv[argptr];1409 cout << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl;1410 argptr+=1;1411 }1412 break;1413 case 'D':1414 ExitFlag = 1;1415 {1416 cout << Verbose(1) << "Depth-First-Search Analysis." << endl;1417 MoleculeLeafClass *Subgraphs = NULL;// list of subgraphs from DFS analysis1418 int *MinimumRingSize = new int[mol->AtomCount];1419 atom ***ListOfLocalAtoms = NULL;1420 int FragmentCounter = 0;1421 class StackClass<bond *> *BackEdgeStack = NULL;1422 class StackClass<bond *> *LocalBackEdgeStack = NULL;1423 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr]), configuration.GetIsAngstroem());1424 mol->CreateListOfBondsPerAtom((ofstream *)&cout);1425 Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);1426 if (Subgraphs != NULL) {1427 Subgraphs->next->FillBondStructureFromReference((ofstream *)&cout, mol, (FragmentCounter = 0), ListOfLocalAtoms, false);// we want to keep the created ListOfLocalAtoms1428 while (Subgraphs->next != NULL) {1429 Subgraphs = Subgraphs->next;1430 LocalBackEdgeStack = new StackClass<bond *> (Subgraphs->Leaf->BondCount);1431 Subgraphs->Leaf->PickLocalBackEdges((ofstream *)&cout, ListOfLocalAtoms[FragmentCounter++], BackEdgeStack, LocalBackEdgeStack);1432 Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize);1433 delete(LocalBackEdgeStack);1434 delete(Subgraphs->previous);1435 }1436 delete(Subgraphs);1437 for (int i=0;i<FragmentCounter;i++)1438 Free((void **)&ListOfLocalAtoms[FragmentCounter], "ParseCommandLineOptions: **ListOfLocalAtoms[]");1439 Free((void **)&ListOfLocalAtoms, "ParseCommandLineOptions: ***ListOfLocalAtoms");1440 }1441 delete(BackEdgeStack);1442 delete[](MinimumRingSize);1443 }1444 //argptr+=1;1445 break;1446 case 'E':1447 ExitFlag = 1;1448 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) {1449 ExitFlag = 255;1450 cerr << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl;1451 } else {1452 SaveFlag = true;1453 cout << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl;1454 first = mol->FindAtom(atoi(argv[argptr]));1455 first->type = periode->FindElement(atoi(argv[argptr+1]));1456 argptr+=2;1457 }1458 break;1459 case 'A':1460 ExitFlag = 1;1461 if ((argptr >= argc) || (argv[argptr][0] == '-')) {1462 ExitFlag =255;1463 cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;1464 } else {1465 cout << "Parsing bonds from " << argv[argptr] << "." << endl;1466 ifstream *input = new ifstream(argv[argptr]);1467 mol->CreateAdjacencyList2((ofstream *)&cout, input);1468 input->close();1469 argptr+=1;1470 }1471 break;1472 case 'N':1473 ExitFlag = 1;1474 if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){1475 ExitFlag = 255;1476 cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;1477 } else {1478 class Tesselation T;1479 string filename(argv[argptr+1]);1480 filename.append(".csv");1481 cout << Verbose(0) << "Evaluating non-convex envelope.";1482 cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;1483 LinkedCell LCList(mol, atof(argv[argptr])*2.);1484 Find_non_convex_border((ofstream *)&cout, mol, &T, &LCList, argv[argptr+1], atof(argv[argptr]));1485 //FindDistributionOfEllipsoids((ofstream *)&cout, &T, &LCList, N, number, filename.c_str());1486 argptr+=2;1487 }1488 break;1489 case 'S':1490 ExitFlag = 1;1491 if ((argptr >= argc) || (argv[argptr][0] == '-')) {1492 ExitFlag = 255;1493 cerr << "Not enough or invalid arguments given for storing tempature: -S <temperature file>" << endl;1494 } else {1495 cout << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;1496 ofstream *output = new ofstream(argv[argptr], ios::trunc);1497 if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output))1498 cout << Verbose(2) << "File could not be written." << endl;1499 else1500 cout << Verbose(2) << "File stored." << endl;1501 output->close();1502 delete(output);1503 argptr+=1;1504 }1505 break;1506 case 'P':1507 ExitFlag = 1;1508 if ((argptr >= argc) || (argv[argptr][0] == '-')) {1509 ExitFlag = 255;1510 cerr << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl;1511 } else {1512 SaveFlag = true;1513 cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;1514 if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem()))1515 cout << Verbose(2) << "File not found." << endl;1516 else1517 cout << Verbose(2) << "File found and parsed." << endl;1518 argptr+=1;1519 }1520 break;1521 case 't':1522 ExitFlag = 1;1523 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {1524 ExitFlag = 255;1525 cerr << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl;1526 } else {1527 ExitFlag = 1;1528 SaveFlag = true;1529 cout << Verbose(1) << "Translating all ions to new origin." << endl;1530 for (int i=NDIM;i--;)1531 x.x[i] = atof(argv[argptr+i]);1532 mol->Translate((const Vector *)&x);1533 argptr+=3;1534 }1535 break;1266 cout << "\t-S <file> Store temperatures from the config file in <file>." << endl; 1267 cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl; 1268 cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl; 1269 cout << "\t-v/-V\t\tGives version information." << endl; 1270 cout << "Note: config files must not begin with '-' !" << endl; 1271 delete(mol); 1272 delete(periode); 1273 return (1); 1274 break; 1275 case 'v': 1276 case 'V': 1277 cout << argv[0] << " " << VERSIONSTRING << endl; 1278 cout << "Build your own molecule position set." << endl; 1279 delete(mol); 1280 delete(periode); 1281 return (1); 1282 break; 1283 case 'e': 1284 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1285 cerr << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl; 1286 } else { 1287 cout << "Using " << argv[argptr] << " as elements database." << endl; 1288 strncpy (configuration.databasepath, argv[argptr], MAXSTRINGSIZE-1); 1289 argptr+=1; 1290 } 1291 break; 1292 case 'n': 1293 cout << "I won't parse trajectories." << endl; 1294 configuration.FastParsing = true; 1295 break; 1296 default: // no match? Step on 1297 argptr++; 1298 break; 1299 } 1300 } else 1301 argptr++; 1302 } while (argptr < argc); 1303 1304 // 2. Parse the element database 1305 if (periode->LoadPeriodentafel(configuration.databasepath)) { 1306 cout << Verbose(0) << "Element list loaded successfully." << endl; 1307 //periode->Output((ofstream *)&cout); 1308 } else { 1309 cout << Verbose(0) << "Element list loading failed." << endl; 1310 return 1; 1311 } 1312 // 3. Find config file name and parse if possible 1313 if (argv[1][0] != '-') { 1314 cout << Verbose(0) << "Config file given." << endl; 1315 test.open(argv[1], ios::in); 1316 if (test == NULL) { 1317 //return (1); 1318 output.open(argv[1], ios::out); 1319 if (output == NULL) { 1320 cout << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl; 1321 config_present = absent; 1322 } else { 1323 cout << "Empty configuration file." << endl; 1324 ConfigFileName = argv[1]; 1325 config_present = empty; 1326 output.close(); 1327 } 1328 } else { 1329 test.close(); 1330 ConfigFileName = argv[1]; 1331 cout << Verbose(1) << "Specified config file found, parsing ... "; 1332 switch (configuration.TestSyntax(ConfigFileName, periode, mol)) { 1333 case 1: 1334 cout << "new syntax." << endl; 1335 configuration.Load(ConfigFileName, periode, mol); 1336 config_present = present; 1337 break; 1338 case 0: 1339 cout << "old syntax." << endl; 1340 configuration.LoadOld(ConfigFileName, periode, mol); 1341 config_present = present; 1342 break; 1343 default: 1344 cout << "Unknown syntax or empty, yet present file." << endl; 1345 config_present = empty; 1346 } 1347 } 1348 } else 1349 config_present = absent; 1350 // 4. parse again through options, now for those depending on elements db and config presence 1351 argptr = 1; 1352 do { 1353 cout << "Current Command line argument: " << argv[argptr] << "." << endl; 1354 if (argv[argptr][0] == '-') { 1355 argptr++; 1356 if ((config_present == present) || (config_present == empty)) { 1357 switch(argv[argptr-1][1]) { 1358 case 'p': 1359 ExitFlag = 1; 1360 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1361 ExitFlag = 255; 1362 cerr << "Not enough arguments for parsing: -p <xyz file>" << endl; 1363 } else { 1364 SaveFlag = true; 1365 cout << Verbose(1) << "Parsing xyz file for new atoms." << endl; 1366 if (!mol->AddXYZFile(argv[argptr])) 1367 cout << Verbose(2) << "File not found." << endl; 1368 else { 1369 cout << Verbose(2) << "File found and parsed." << endl; 1370 config_present = present; 1371 } 1372 } 1373 break; 1374 case 'a': 1375 ExitFlag = 1; 1376 if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1]))) { 1377 ExitFlag = 255; 1378 cerr << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl; 1379 } else { 1380 SaveFlag = true; 1381 cout << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), "; 1382 first = new atom; 1383 first->type = periode->FindElement(atoi(argv[argptr])); 1384 if (first->type != NULL) 1385 cout << Verbose(2) << "found element " << first->type->name << endl; 1386 for (int i=NDIM;i--;) 1387 first->x.x[i] = atof(argv[argptr+1+i]); 1388 if (first->type != NULL) { 1389 mol->AddAtom(first); // add to molecule 1390 if ((config_present == empty) && (mol->AtomCount != 0)) 1391 config_present = present; 1392 } else 1393 cerr << Verbose(1) << "Could not find the specified element." << endl; 1394 argptr+=4; 1395 } 1396 break; 1397 default: // no match? Don't step on (this is done in next switch's default) 1398 break; 1399 } 1400 } 1401 if (config_present == present) { 1402 switch(argv[argptr-1][1]) { 1403 case 'B': 1404 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1405 ExitFlag = 255; 1406 cerr << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl; 1407 } else { 1408 configuration.basis = argv[argptr]; 1409 cout << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl; 1410 argptr+=1; 1411 } 1412 break; 1413 case 'D': 1414 ExitFlag = 1; 1415 { 1416 cout << Verbose(1) << "Depth-First-Search Analysis." << endl; 1417 MoleculeLeafClass *Subgraphs = NULL; // list of subgraphs from DFS analysis 1418 int *MinimumRingSize = new int[mol->AtomCount]; 1419 atom ***ListOfLocalAtoms = NULL; 1420 int FragmentCounter = 0; 1421 class StackClass<bond *> *BackEdgeStack = NULL; 1422 class StackClass<bond *> *LocalBackEdgeStack = NULL; 1423 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr]), configuration.GetIsAngstroem()); 1424 mol->CreateListOfBondsPerAtom((ofstream *)&cout); 1425 Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack); 1426 if (Subgraphs != NULL) { 1427 Subgraphs->next->FillBondStructureFromReference((ofstream *)&cout, mol, (FragmentCounter = 0), ListOfLocalAtoms, false); // we want to keep the created ListOfLocalAtoms 1428 while (Subgraphs->next != NULL) { 1429 Subgraphs = Subgraphs->next; 1430 LocalBackEdgeStack = new StackClass<bond *> (Subgraphs->Leaf->BondCount); 1431 Subgraphs->Leaf->PickLocalBackEdges((ofstream *)&cout, ListOfLocalAtoms[FragmentCounter++], BackEdgeStack, LocalBackEdgeStack); 1432 Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize); 1433 delete(LocalBackEdgeStack); 1434 delete(Subgraphs->previous); 1435 } 1436 delete(Subgraphs); 1437 for (int i=0;i<FragmentCounter;i++) 1438 Free((void **)&ListOfLocalAtoms[FragmentCounter], "ParseCommandLineOptions: **ListOfLocalAtoms[]"); 1439 Free((void **)&ListOfLocalAtoms, "ParseCommandLineOptions: ***ListOfLocalAtoms"); 1440 } 1441 delete(BackEdgeStack); 1442 delete[](MinimumRingSize); 1443 } 1444 //argptr+=1; 1445 break; 1446 case 'E': 1447 ExitFlag = 1; 1448 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) { 1449 ExitFlag = 255; 1450 cerr << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl; 1451 } else { 1452 SaveFlag = true; 1453 cout << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl; 1454 first = mol->FindAtom(atoi(argv[argptr])); 1455 first->type = periode->FindElement(atoi(argv[argptr+1])); 1456 argptr+=2; 1457 } 1458 break; 1459 case 'A': 1460 ExitFlag = 1; 1461 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1462 ExitFlag =255; 1463 cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl; 1464 } else { 1465 cout << "Parsing bonds from " << argv[argptr] << "." << endl; 1466 ifstream *input = new ifstream(argv[argptr]); 1467 mol->CreateAdjacencyList2((ofstream *)&cout, input); 1468 input->close(); 1469 argptr+=1; 1470 } 1471 break; 1472 case 'N': 1473 ExitFlag = 1; 1474 if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){ 1475 ExitFlag = 255; 1476 cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl; 1477 } else { 1478 class Tesselation T; 1479 string filename(argv[argptr+1]); 1480 filename.append(".csv"); 1481 cout << Verbose(0) << "Evaluating non-convex envelope."; 1482 cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl; 1483 LinkedCell LCList(mol, atof(argv[argptr])*2.); 1484 Find_non_convex_border((ofstream *)&cout, mol, &T, &LCList, argv[argptr+1], atof(argv[argptr])); 1485 //FindDistributionOfEllipsoids((ofstream *)&cout, &T, &LCList, N, number, filename.c_str()); 1486 argptr+=2; 1487 } 1488 break; 1489 case 'S': 1490 ExitFlag = 1; 1491 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1492 ExitFlag = 255; 1493 cerr << "Not enough or invalid arguments given for storing tempature: -S <temperature file>" << endl; 1494 } else { 1495 cout << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl; 1496 ofstream *output = new ofstream(argv[argptr], ios::trunc); 1497 if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output)) 1498 cout << Verbose(2) << "File could not be written." << endl; 1499 else 1500 cout << Verbose(2) << "File stored." << endl; 1501 output->close(); 1502 delete(output); 1503 argptr+=1; 1504 } 1505 break; 1506 case 'P': 1507 ExitFlag = 1; 1508 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1509 ExitFlag = 255; 1510 cerr << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl; 1511 } else { 1512 SaveFlag = true; 1513 cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl; 1514 if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem())) 1515 cout << Verbose(2) << "File not found." << endl; 1516 else 1517 cout << Verbose(2) << "File found and parsed." << endl; 1518 argptr+=1; 1519 } 1520 break; 1521 case 't': 1522 ExitFlag = 1; 1523 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1524 ExitFlag = 255; 1525 cerr << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl; 1526 } else { 1527 ExitFlag = 1; 1528 SaveFlag = true; 1529 cout << Verbose(1) << "Translating all ions to new origin." << endl; 1530 for (int i=NDIM;i--;) 1531 x.x[i] = atof(argv[argptr+i]); 1532 mol->Translate((const Vector *)&x); 1533 argptr+=3; 1534 } 1535 break; 1536 1536 case 'T': 1537 1537 ExitFlag = 1; … … 1549 1549 } 1550 1550 break; 1551 case 's':1552 ExitFlag = 1;1553 if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {1554 ExitFlag = 255;1555 cerr << "Not enough or invalid arguments given for scaling: -s <factor/[factor_x]> [factor_y] [factor_z]" << endl;1556 } else {1557 SaveFlag = true;1558 j = -1;1559 cout << Verbose(1) << "Scaling all ion positions by factor." << endl;1560 factor = new double[NDIM];1561 factor[0] = atof(argv[argptr]);1562 if ((argptr < argc) && (IsValidNumber(argv[argptr])))1563 argptr++;1564 factor[1] = atof(argv[argptr]);1565 if ((argptr < argc) && (IsValidNumber(argv[argptr])))1566 argptr++;1567 factor[2] = atof(argv[argptr]);1568 mol->Scale(&factor);1569 for (int i=0;i<NDIM;i++) {1570 j += i+1;1571 x.x[i] = atof(argv[NDIM+i]);1572 mol->cell_size[j]*=factor[i];1573 }1574 delete[](factor);1575 argptr+=1;1576 }1577 break;1578 case 'b':1579 ExitFlag = 1;1580 if ((argptr+5 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {1581 ExitFlag = 255;1582 cerr << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl;1583 } else {1584 SaveFlag = true;1585 cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;1586 for (int i=0;i<6;i++) {1587 mol->cell_size[i] = atof(argv[argptr+i]);1588 }1589 // center1590 mol->CenterInBox((ofstream *)&cout);1551 case 's': 1552 ExitFlag = 1; 1553 if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) { 1554 ExitFlag = 255; 1555 cerr << "Not enough or invalid arguments given for scaling: -s <factor/[factor_x]> [factor_y] [factor_z]" << endl; 1556 } else { 1557 SaveFlag = true; 1558 j = -1; 1559 cout << Verbose(1) << "Scaling all ion positions by factor." << endl; 1560 factor = new double[NDIM]; 1561 factor[0] = atof(argv[argptr]); 1562 if ((argptr < argc) && (IsValidNumber(argv[argptr]))) 1563 argptr++; 1564 factor[1] = atof(argv[argptr]); 1565 if ((argptr < argc) && (IsValidNumber(argv[argptr]))) 1566 argptr++; 1567 factor[2] = atof(argv[argptr]); 1568 mol->Scale(&factor); 1569 for (int i=0;i<NDIM;i++) { 1570 j += i+1; 1571 x.x[i] = atof(argv[NDIM+i]); 1572 mol->cell_size[j]*=factor[i]; 1573 } 1574 delete[](factor); 1575 argptr+=1; 1576 } 1577 break; 1578 case 'b': 1579 ExitFlag = 1; 1580 if ((argptr+5 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) { 1581 ExitFlag = 255; 1582 cerr << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl; 1583 } else { 1584 SaveFlag = true; 1585 cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl; 1586 for (int i=0;i<6;i++) { 1587 mol->cell_size[i] = atof(argv[argptr+i]); 1588 } 1589 // center 1590 mol->CenterInBox((ofstream *)&cout); 1591 1591 argptr+=6; 1592 }1593 break;1594 case 'c':1595 ExitFlag = 1;1596 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {1597 ExitFlag = 255;1598 cerr << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl;1599 } else {1600 SaveFlag = true;1601 j = -1;1602 cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;1603 // make every coordinate positive1604 mol->CenterEdge((ofstream *)&cout, &x);1605 // update Box of atoms by boundary1606 mol->SetBoxDimension(&x);1607 // translate each coordinate by boundary1608 j=-1;1609 for (int i=0;i<NDIM;i++) {1610 j += i+1;1611 x.x[i] = atof(argv[argptr++]);1612 mol->cell_size[j] += x.x[i]*2.;1613 }1614 mol->Translate((const Vector *)&x);1592 } 1593 break; 1594 case 'c': 1595 ExitFlag = 1; 1596 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1597 ExitFlag = 255; 1598 cerr << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl; 1599 } else { 1600 SaveFlag = true; 1601 j = -1; 1602 cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl; 1603 // make every coordinate positive 1604 mol->CenterEdge((ofstream *)&cout, &x); 1605 // update Box of atoms by boundary 1606 mol->SetBoxDimension(&x); 1607 // translate each coordinate by boundary 1608 j=-1; 1609 for (int i=0;i<NDIM;i++) { 1610 j += i+1; 1611 x.x[i] = atof(argv[argptr++]); 1612 mol->cell_size[j] += x.x[i]*2.; 1613 } 1614 mol->Translate((const Vector *)&x); 1615 1615 argptr+=3; 1616 }1617 break;1618 case 'O':1619 ExitFlag = 1;1620 SaveFlag = true;1621 cout << Verbose(1) << "Centering atoms in origin." << endl;1622 mol->CenterOrigin((ofstream *)&cout, &x);1623 mol->SetBoxDimension(&x);1616 } 1617 break; 1618 case 'O': 1619 ExitFlag = 1; 1620 SaveFlag = true; 1621 cout << Verbose(1) << "Centering atoms in origin." << endl; 1622 mol->CenterOrigin((ofstream *)&cout, &x); 1623 mol->SetBoxDimension(&x); 1624 1624 argptr+=0; 1625 break;1626 case 'r':1627 ExitFlag = 1;1628 SaveFlag = true;1629 cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;1630 break;1631 case 'F':1632 case 'f':1633 ExitFlag = 1;1634 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {1635 ExitFlag = 255;1636 cerr << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl;1637 } else {1638 cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;1639 cout << Verbose(0) << "Creating connection matrix..." << endl;1640 start = clock();1641 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr++]), configuration.GetIsAngstroem());1642 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;1643 if (mol->first->next != mol->last) {1644 ExitFlag = mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);1645 }1646 end = clock();1647 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;1648 argptr+=2;1649 }1650 break;1651 case 'm':1652 ExitFlag = 1;1653 j = atoi(argv[argptr++]);1654 if ((j<0) || (j>1)) {1655 cerr << Verbose(1) << "ERROR: Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;1656 j = 0;1657 }1658 if (j) {1659 SaveFlag = true;1660 cout << Verbose(0) << "Converting to prinicipal axis system." << endl;1661 } else1662 cout << Verbose(0) << "Evaluating prinicipal axis." << endl;1663 mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j);1664 break;1665 case 'o':1666 ExitFlag = 1;1667 if ((argptr >= argc) || (argv[argptr][0] == '-')){1668 ExitFlag = 255;1669 cerr << "Not enough or invalid arguments given for convex envelope: -o <tecplot output file>" << endl;1670 } else {1671 cout << Verbose(0) << "Evaluating volume of the convex envelope.";1672 cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;1673 VolumeOfConvexEnvelope((ofstream *)&cout, argv[argptr], &configuration, NULL, mol);1674 argptr+=1;1675 }1676 break;1677 case 'U':1678 ExitFlag = 1;1679 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {1680 ExitFlag = 255;1681 cerr << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl;1682 volume = -1; // for case 'u': don't print error again1683 } else {1684 volume = atof(argv[argptr++]);1685 cout << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;1686 }1687 case 'u':1688 ExitFlag = 1;1689 if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {1690 if (volume != -1)1691 ExitFlag = 255;1692 cerr << "Not enough arguments given for suspension: -u <density>" << endl;1693 } else {1694 double density;1695 SaveFlag = true;1696 cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";1697 density = atof(argv[argptr++]);1698 if (density < 1.0) {1699 cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;1700 density = 1.3;1701 }1702 // for(int i=0;i<NDIM;i++) {1703 // repetition[i] = atoi(argv[argptr++]);1704 // if (repetition[i] < 1)1705 // cerr << Verbose(0) << "ERROR: repetition value must be greater 1!" << endl;1706 // repetition[i] = 1;1707 // }1708 PrepareClustersinWater((ofstream *)&cout, &configuration, mol, volume, density);// if volume == 0, will calculate from ConvexEnvelope1709 }1710 break;1711 case 'd':1712 ExitFlag = 1;1713 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {1714 ExitFlag = 255;1715 cerr << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl;1716 } else {1717 SaveFlag = true;1718 for (int axis = 1; axis <= NDIM; axis++) {1719 int faktor = atoi(argv[argptr++]);1720 int count;1721 element ** Elements;1722 Vector ** vectors;1723 if (faktor < 1) {1724 cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;1725 faktor = 1;1726 }1727 mol->CountAtoms((ofstream *)&cout);// recount atoms1728 if (mol->AtomCount != 0) {// if there is more than none1729 count = mol->AtomCount;// is changed becausing of adding, thus has to be stored away beforehand1730 Elements = new element *[count];1731 vectors = new Vector *[count];1732 j = 0;1733 first = mol->start;1734 while (first->next != mol->end) {// make a list of all atoms with coordinates and element1735 first = first->next;1736 Elements[j] = first->type;1737 vectors[j] = &first->x;1738 j++;1739 }1740 if (count != j)1741 cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;1742 x.Zero();1743 y.Zero();1744 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 magnitude1745 for (int i=1;i<faktor;i++) {// then add this list with respective translation factor times1746 x.AddVector(&y); // per factor one cell width further1747 for (int k=count;k--;) { // go through every atom of the original cell1748 first = new atom(); // create a new body1749 first->x.CopyVector(vectors[k]);// use coordinate of original atom1750 first->x.AddVector(&x);// translate the coordinates1751 first->type = Elements[k];// insert original element1752 mol->AddAtom(first);// and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)1753 }1754 }1755 // free memory1756 delete[](Elements);1757 delete[](vectors);1758 // correct cell size1759 if (axis < 0) { // if sign was negative, we have to translate everything1760 x.Zero();1761 x.AddVector(&y);1762 x.Scale(-(faktor-1));1763 mol->Translate(&x);1764 }1765 mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;1766 }1767 }1768 }1769 break;1770 default:// no match? Step on1771 if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step!1772 argptr++;1773 break;1774 }1775 }1776 } else argptr++;1777 } while (argptr < argc);1778 if (SaveFlag)1779 SaveConfig(ConfigFileName, &configuration, periode, molecules);1780 } else {// no arguments, hence scan the elements db1781 if (periode->LoadPeriodentafel(configuration.databasepath))1782 cout << Verbose(0) << "Element list loaded successfully." << endl;1783 else1784 cout << Verbose(0) << "Element list loading failed." << endl;1785 configuration.RetrieveConfigPathAndName("main_pcp_linux");1786 }1787 return(ExitFlag);1625 break; 1626 case 'r': 1627 ExitFlag = 1; 1628 SaveFlag = true; 1629 cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl; 1630 break; 1631 case 'F': 1632 case 'f': 1633 ExitFlag = 1; 1634 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) { 1635 ExitFlag = 255; 1636 cerr << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl; 1637 } else { 1638 cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl; 1639 cout << Verbose(0) << "Creating connection matrix..." << endl; 1640 start = clock(); 1641 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr++]), configuration.GetIsAngstroem()); 1642 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl; 1643 if (mol->first->next != mol->last) { 1644 ExitFlag = mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration); 1645 } 1646 end = clock(); 1647 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl; 1648 argptr+=2; 1649 } 1650 break; 1651 case 'm': 1652 ExitFlag = 1; 1653 j = atoi(argv[argptr++]); 1654 if ((j<0) || (j>1)) { 1655 cerr << Verbose(1) << "ERROR: Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl; 1656 j = 0; 1657 } 1658 if (j) { 1659 SaveFlag = true; 1660 cout << Verbose(0) << "Converting to prinicipal axis system." << endl; 1661 } else 1662 cout << Verbose(0) << "Evaluating prinicipal axis." << endl; 1663 mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j); 1664 break; 1665 case 'o': 1666 ExitFlag = 1; 1667 if ((argptr >= argc) || (argv[argptr][0] == '-')){ 1668 ExitFlag = 255; 1669 cerr << "Not enough or invalid arguments given for convex envelope: -o <tecplot output file>" << endl; 1670 } else { 1671 cout << Verbose(0) << "Evaluating volume of the convex envelope."; 1672 cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl; 1673 VolumeOfConvexEnvelope((ofstream *)&cout, argv[argptr], &configuration, NULL, mol); 1674 argptr+=1; 1675 } 1676 break; 1677 case 'U': 1678 ExitFlag = 1; 1679 if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) { 1680 ExitFlag = 255; 1681 cerr << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl; 1682 volume = -1; // for case 'u': don't print error again 1683 } else { 1684 volume = atof(argv[argptr++]); 1685 cout << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl; 1686 } 1687 case 'u': 1688 ExitFlag = 1; 1689 if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) { 1690 if (volume != -1) 1691 ExitFlag = 255; 1692 cerr << "Not enough arguments given for suspension: -u <density>" << endl; 1693 } else { 1694 double density; 1695 SaveFlag = true; 1696 cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water."; 1697 density = atof(argv[argptr++]); 1698 if (density < 1.0) { 1699 cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl; 1700 density = 1.3; 1701 } 1702 // for(int i=0;i<NDIM;i++) { 1703 // repetition[i] = atoi(argv[argptr++]); 1704 // if (repetition[i] < 1) 1705 // cerr << Verbose(0) << "ERROR: repetition value must be greater 1!" << endl; 1706 // repetition[i] = 1; 1707 // } 1708 PrepareClustersinWater((ofstream *)&cout, &configuration, mol, volume, density); // if volume == 0, will calculate from ConvexEnvelope 1709 } 1710 break; 1711 case 'd': 1712 ExitFlag = 1; 1713 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1714 ExitFlag = 255; 1715 cerr << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl; 1716 } else { 1717 SaveFlag = true; 1718 for (int axis = 1; axis <= NDIM; axis++) { 1719 int faktor = atoi(argv[argptr++]); 1720 int count; 1721 element ** Elements; 1722 Vector ** vectors; 1723 if (faktor < 1) { 1724 cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl; 1725 faktor = 1; 1726 } 1727 mol->CountAtoms((ofstream *)&cout); // recount atoms 1728 if (mol->AtomCount != 0) { // if there is more than none 1729 count = mol->AtomCount; // is changed becausing of adding, thus has to be stored away beforehand 1730 Elements = new element *[count]; 1731 vectors = new Vector *[count]; 1732 j = 0; 1733 first = mol->start; 1734 while (first->next != mol->end) { // make a list of all atoms with coordinates and element 1735 first = first->next; 1736 Elements[j] = first->type; 1737 vectors[j] = &first->x; 1738 j++; 1739 } 1740 if (count != j) 1741 cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl; 1742 x.Zero(); 1743 y.Zero(); 1744 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 1745 for (int i=1;i<faktor;i++) { // then add this list with respective translation factor times 1746 x.AddVector(&y); // per factor one cell width further 1747 for (int k=count;k--;) { // go through every atom of the original cell 1748 first = new atom(); // create a new body 1749 first->x.CopyVector(vectors[k]); // use coordinate of original atom 1750 first->x.AddVector(&x); // translate the coordinates 1751 first->type = Elements[k]; // insert original element 1752 mol->AddAtom(first); // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...) 1753 } 1754 } 1755 // free memory 1756 delete[](Elements); 1757 delete[](vectors); 1758 // correct cell size 1759 if (axis < 0) { // if sign was negative, we have to translate everything 1760 x.Zero(); 1761 x.AddVector(&y); 1762 x.Scale(-(faktor-1)); 1763 mol->Translate(&x); 1764 } 1765 mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor; 1766 } 1767 } 1768 } 1769 break; 1770 default: // no match? Step on 1771 if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step! 1772 argptr++; 1773 break; 1774 } 1775 } 1776 } else argptr++; 1777 } while (argptr < argc); 1778 if (SaveFlag) 1779 SaveConfig(ConfigFileName, &configuration, periode, molecules); 1780 } else { // no arguments, hence scan the elements db 1781 if (periode->LoadPeriodentafel(configuration.databasepath)) 1782 cout << Verbose(0) << "Element list loaded successfully." << endl; 1783 else 1784 cout << Verbose(0) << "Element list loading failed." << endl; 1785 configuration.RetrieveConfigPathAndName("main_pcp_linux"); 1786 } 1787 return(ExitFlag); 1788 1788 }; 1789 1789 … … 1792 1792 int main(int argc, char **argv) 1793 1793 { 1794 periodentafel *periode = new periodentafel; // and a period table of all elements1795 MoleculeListClass *molecules = new MoleculeListClass; // list of all molecules1796 molecule *mol = NULL;1797 config configuration;1798 char choice;// menu choice char1799 Vector x,y,z,n;// coordinates for absolute point in cell volume1800 ifstream test;1801 ofstream output;1802 string line;1803 char *ConfigFileName = NULL;1804 int j, count;1805 1806 // =========================== PARSE COMMAND LINE OPTIONS ====================================1807 j = ParseCommandLineOptions(argc, argv, molecules, periode, configuration, ConfigFileName);1808 switch(j) {1809 case 0: // something went wrong1810 delete(molecules); // also free's all molecules contained1811 delete(periode);1812 return j;1813 break;1814 case 1: // just for -v and -h options1815 delete(molecules); // also free's all molecules contained1816 delete(periode);1817 return 0;1818 break;1819 default:1820 break;1821 }1822 1823 // General stuff1824 if (molecules->ListOfMolecules.size() == 0) {1794 periodentafel *periode = new periodentafel; // and a period table of all elements 1795 MoleculeListClass *molecules = new MoleculeListClass; // list of all molecules 1796 molecule *mol = NULL; 1797 config configuration; 1798 char choice; // menu choice char 1799 Vector x,y,z,n; // coordinates for absolute point in cell volume 1800 ifstream test; 1801 ofstream output; 1802 string line; 1803 char *ConfigFileName = NULL; 1804 int j, count; 1805 1806 // =========================== PARSE COMMAND LINE OPTIONS ==================================== 1807 j = ParseCommandLineOptions(argc, argv, molecules, periode, configuration, ConfigFileName); 1808 switch(j) { 1809 case 0: // something went wrong 1810 delete(molecules); // also free's all molecules contained 1811 delete(periode); 1812 return j; 1813 break; 1814 case 1: // just for -v and -h options 1815 delete(molecules); // also free's all molecules contained 1816 delete(periode); 1817 return 0; 1818 break; 1819 default: 1820 break; 1821 } 1822 1823 // General stuff 1824 if (molecules->ListOfMolecules.size() == 0) { 1825 1825 mol = new molecule(periode); 1826 1826 if (mol->cell_size[0] == 0.) { … … 1832 1832 } 1833 1833 molecules->insert(mol); 1834 }1835 1836 // =========================== START INTERACTIVE SESSION ====================================1837 1838 // now the main construction loop1839 cout << Verbose(0) << endl << "Now comes the real construction..." << endl;1840 do {1841 cout << Verbose(0) << endl << endl;1842 cout << Verbose(0) << "============Molecule list=======================" << endl;1843 molecules->Enumerate((ofstream *)&cout);1844 cout << Verbose(0) << "============Menu===============================" << endl;1834 } 1835 1836 // =========================== START INTERACTIVE SESSION ==================================== 1837 1838 // now the main construction loop 1839 cout << Verbose(0) << endl << "Now comes the real construction..." << endl; 1840 do { 1841 cout << Verbose(0) << endl << endl; 1842 cout << Verbose(0) << "============Molecule list=======================" << endl; 1843 molecules->Enumerate((ofstream *)&cout); 1844 cout << Verbose(0) << "============Menu===============================" << endl; 1845 1845 cout << Verbose(0) << "a - set molecule (in)active" << endl; 1846 1846 cout << Verbose(0) << "e - edit new molecules" << endl; … … 1858 1858 cin >> choice; 1859 1859 1860 switch (choice) {1861 case 'a': // (in)activate molecule1860 switch (choice) { 1861 case 'a': // (in)activate molecule 1862 1862 { 1863 1863 cout << "Enter index of molecule: "; … … 1869 1869 (*ListRunner)->ActiveFlag = !(*ListRunner)->ActiveFlag; 1870 1870 } 1871 break;1872 1873 case 'c': // edit each field of the configuration1874 configuration.Edit();1875 break;1871 break; 1872 1873 case 'c': // edit each field of the configuration 1874 configuration.Edit(); 1875 break; 1876 1876 1877 1877 case 'e': // create molecule … … 1891 1891 break; 1892 1892 1893 case 'q': // quit1894 break;1895 1896 case 's': // save to config file1897 SaveConfig(ConfigFileName, &configuration, periode, molecules);1898 break;1899 1900 case 'T':1901 testroutine(molecules);1902 break;1903 1904 default:1905 break;1906 };1907 } while (choice != 'q');1908 1909 // save element data base1910 if (periode->StorePeriodentafel(configuration.databasepath)) //ElementsFileName1911 cout << Verbose(0) << "Saving of elements.db successful." << endl;1912 else1913 cout << Verbose(0) << "Saving of elements.db failed." << endl;1914 1915 delete(molecules); // also free's all molecules contained1916 delete(periode);1917 return (0);1893 case 'q': // quit 1894 break; 1895 1896 case 's': // save to config file 1897 SaveConfig(ConfigFileName, &configuration, periode, molecules); 1898 break; 1899 1900 case 'T': 1901 testroutine(molecules); 1902 break; 1903 1904 default: 1905 break; 1906 }; 1907 } while (choice != 'q'); 1908 1909 // save element data base 1910 if (periode->StorePeriodentafel(configuration.databasepath)) //ElementsFileName 1911 cout << Verbose(0) << "Saving of elements.db successful." << endl; 1912 else 1913 cout << Verbose(0) << "Saving of elements.db failed." << endl; 1914 1915 delete(molecules); // also free's all molecules contained 1916 delete(periode); 1917 return (0); 1918 1918 } 1919 1919
Note:
See TracChangeset
for help on using the changeset viewer.
