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