| 36 |
|
* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). |
| 37 |
|
* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). |
| 38 |
|
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). |
| 39 |
< |
* [4] Vardeman & Gezelter, in progress (2009). |
| 39 |
> |
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
| 40 |
> |
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
| 41 |
|
*/ |
| 42 |
|
|
| 43 |
|
#define _LARGEFILE_SOURCE64 |
| 58 |
|
#include "utils/simError.h" |
| 59 |
|
#include "utils/MemoryUtils.hpp" |
| 60 |
|
#include "utils/StringTokenizer.hpp" |
| 61 |
+ |
#include "brains/Thermo.hpp" |
| 62 |
|
|
| 63 |
|
#ifdef IS_MPI |
| 62 |
– |
|
| 64 |
|
#include <mpi.h> |
| 65 |
< |
#define TAKE_THIS_TAG_CHAR 0 |
| 65 |
< |
#define TAKE_THIS_TAG_INT 1 |
| 65 |
> |
#endif |
| 66 |
|
|
| 67 |
– |
#endif // is_mpi |
| 67 |
|
|
| 69 |
– |
|
| 68 |
|
namespace OpenMD { |
| 69 |
|
|
| 70 |
|
DumpReader::DumpReader(SimInfo* info, const std::string& filename) |
| 140 |
|
prevPos = currPos; |
| 141 |
|
bool foundOpenSnapshotTag = false; |
| 142 |
|
bool foundClosedSnapshotTag = false; |
| 143 |
+ |
bool foundOpenSiteDataTag = false; |
| 144 |
|
while(inFile_->getline(buffer, bufferSize)) { |
| 145 |
|
++lineNo; |
| 146 |
|
|
| 242 |
|
|
| 243 |
|
if (needCOMprops_) { |
| 244 |
|
Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot(); |
| 245 |
+ |
Thermo thermo(info_); |
| 246 |
|
Vector3d com; |
| 247 |
– |
Vector3d comvel; |
| 248 |
– |
Vector3d comw; |
| 249 |
– |
if (needPos_ && needVel_){ |
| 250 |
– |
info_->getComAll(com, comvel); |
| 251 |
– |
comw = info_->getAngularMomentum(); |
| 252 |
– |
}else{ |
| 253 |
– |
com = info_->getCom(); |
| 254 |
– |
comvel = 0.0; |
| 255 |
– |
comw = 0.0; |
| 256 |
– |
} |
| 257 |
– |
s->setCOMprops(com, comvel, comw); |
| 258 |
– |
} |
| 247 |
|
|
| 248 |
+ |
if (needPos_ && needVel_) { |
| 249 |
+ |
Vector3d comvel; |
| 250 |
+ |
Vector3d comw; |
| 251 |
+ |
thermo.getComAll(com, comvel); |
| 252 |
+ |
comw = thermo.getAngularMomentum(); |
| 253 |
+ |
} else { |
| 254 |
+ |
com = thermo.getCom(); |
| 255 |
+ |
} |
| 256 |
+ |
} |
| 257 |
|
} |
| 258 |
|
|
| 259 |
|
void DumpReader::readSet(int whichFrame) { |
| 321 |
|
|
| 322 |
|
inputStream.getline(buffer, bufferSize); |
| 323 |
|
line = buffer; |
| 324 |
< |
if (line.find("</Snapshot>") == std::string::npos) { |
| 325 |
< |
sprintf(painCave.errMsg, |
| 326 |
< |
"DumpReader Error: can not find </Snapshot>\n"); |
| 327 |
< |
painCave.isFatal = 1; |
| 328 |
< |
simError(); |
| 329 |
< |
} |
| 330 |
< |
|
| 324 |
> |
|
| 325 |
> |
if (line.find("<SiteData>") != std::string::npos) { |
| 326 |
> |
//read SiteData |
| 327 |
> |
readSiteData(inputStream); |
| 328 |
> |
} else { |
| 329 |
> |
if (line.find("</Snapshot>") == std::string::npos) { |
| 330 |
> |
sprintf(painCave.errMsg, |
| 331 |
> |
"DumpReader Error: can not find </Snapshot>\n"); |
| 332 |
> |
painCave.isFatal = 1; |
| 333 |
> |
simError(); |
| 334 |
> |
} |
| 335 |
> |
} |
| 336 |
|
} |
| 337 |
|
|
| 338 |
|
void DumpReader::parseDumpLine(const std::string& line) { |
| 468 |
|
integrableObject->setTrq(torque); |
| 469 |
|
break; |
| 470 |
|
} |
| 471 |
+ |
case 'u' : { |
| 472 |
+ |
|
| 473 |
+ |
RealType particlePot; |
| 474 |
+ |
particlePot = tokenizer.nextTokenAsDouble(); |
| 475 |
+ |
integrableObject->setParticlePot(particlePot); |
| 476 |
+ |
break; |
| 477 |
+ |
} |
| 478 |
+ |
case 'c' : { |
| 479 |
+ |
|
| 480 |
+ |
RealType flucQPos; |
| 481 |
+ |
flucQPos = tokenizer.nextTokenAsDouble(); |
| 482 |
+ |
integrableObject->setFlucQPos(flucQPos); |
| 483 |
+ |
break; |
| 484 |
+ |
} |
| 485 |
+ |
case 'w' : { |
| 486 |
+ |
|
| 487 |
+ |
RealType flucQVel; |
| 488 |
+ |
flucQVel = tokenizer.nextTokenAsDouble(); |
| 489 |
+ |
integrableObject->setFlucQVel(flucQVel); |
| 490 |
+ |
break; |
| 491 |
+ |
} |
| 492 |
+ |
case 'g' : { |
| 493 |
+ |
|
| 494 |
+ |
RealType flucQFrc; |
| 495 |
+ |
flucQFrc = tokenizer.nextTokenAsDouble(); |
| 496 |
+ |
integrableObject->setFlucQFrc(flucQFrc); |
| 497 |
+ |
break; |
| 498 |
+ |
} |
| 499 |
+ |
case 'e' : { |
| 500 |
+ |
|
| 501 |
+ |
Vector3d eField; |
| 502 |
+ |
eField[0] = tokenizer.nextTokenAsDouble(); |
| 503 |
+ |
eField[1] = tokenizer.nextTokenAsDouble(); |
| 504 |
+ |
eField[2] = tokenizer.nextTokenAsDouble(); |
| 505 |
+ |
integrableObject->setElectricField(eField); |
| 506 |
+ |
break; |
| 507 |
+ |
} |
| 508 |
|
default: { |
| 509 |
|
sprintf(painCave.errMsg, |
| 510 |
|
"DumpReader Error: %s is an unrecognized type\n", type.c_str()); |
| 519 |
|
} |
| 520 |
|
|
| 521 |
|
|
| 522 |
< |
void DumpReader::readStuntDoubles(std::istream& inputStream) { |
| 522 |
> |
void DumpReader::parseSiteLine(const std::string& line) { |
| 523 |
> |
|
| 524 |
> |
StringTokenizer tokenizer(line); |
| 525 |
> |
int nTokens; |
| 526 |
> |
|
| 527 |
> |
nTokens = tokenizer.countTokens(); |
| 528 |
> |
|
| 529 |
> |
if (nTokens < 2) { |
| 530 |
> |
sprintf(painCave.errMsg, |
| 531 |
> |
"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str()); |
| 532 |
> |
painCave.isFatal = 1; |
| 533 |
> |
simError(); |
| 534 |
> |
} |
| 535 |
> |
|
| 536 |
> |
/** |
| 537 |
> |
* The first token is the global integrable object index. |
| 538 |
> |
*/ |
| 539 |
> |
|
| 540 |
> |
int index = tokenizer.nextTokenAsInt(); |
| 541 |
> |
StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index); |
| 542 |
> |
if (integrableObject == NULL) { |
| 543 |
> |
return; |
| 544 |
> |
} |
| 545 |
> |
StuntDouble* sd = integrableObject; |
| 546 |
|
|
| 547 |
+ |
/** |
| 548 |
+ |
* Test to see if the next token is an integer or not. If not, |
| 549 |
+ |
* we've got data on the integrable object itself. If there is an |
| 550 |
+ |
* integer, we're parsing data for a site on a rigid body. |
| 551 |
+ |
*/ |
| 552 |
+ |
|
| 553 |
+ |
std::string indexTest = tokenizer.peekNextToken(); |
| 554 |
+ |
std::istringstream i(indexTest); |
| 555 |
+ |
int siteIndex; |
| 556 |
+ |
if (i >> siteIndex) { |
| 557 |
+ |
// chew up this token and parse as an int: |
| 558 |
+ |
siteIndex = tokenizer.nextTokenAsInt(); |
| 559 |
+ |
RigidBody* rb = static_cast<RigidBody*>(integrableObject); |
| 560 |
+ |
sd = rb->getAtoms()[siteIndex]; |
| 561 |
+ |
} |
| 562 |
+ |
|
| 563 |
+ |
/** |
| 564 |
+ |
* The next token contains information on what follows. |
| 565 |
+ |
*/ |
| 566 |
+ |
std::string type = tokenizer.nextToken(); |
| 567 |
+ |
int size = type.size(); |
| 568 |
+ |
|
| 569 |
+ |
for(int i = 0; i < size; ++i) { |
| 570 |
+ |
switch(type[i]) { |
| 571 |
+ |
|
| 572 |
+ |
case 'u' : { |
| 573 |
+ |
|
| 574 |
+ |
RealType particlePot; |
| 575 |
+ |
particlePot = tokenizer.nextTokenAsDouble(); |
| 576 |
+ |
sd->setParticlePot(particlePot); |
| 577 |
+ |
break; |
| 578 |
+ |
} |
| 579 |
+ |
case 'c' : { |
| 580 |
+ |
|
| 581 |
+ |
RealType flucQPos; |
| 582 |
+ |
flucQPos = tokenizer.nextTokenAsDouble(); |
| 583 |
+ |
sd->setFlucQPos(flucQPos); |
| 584 |
+ |
break; |
| 585 |
+ |
} |
| 586 |
+ |
case 'w' : { |
| 587 |
+ |
|
| 588 |
+ |
RealType flucQVel; |
| 589 |
+ |
flucQVel = tokenizer.nextTokenAsDouble(); |
| 590 |
+ |
sd->setFlucQVel(flucQVel); |
| 591 |
+ |
break; |
| 592 |
+ |
} |
| 593 |
+ |
case 'g' : { |
| 594 |
+ |
|
| 595 |
+ |
RealType flucQFrc; |
| 596 |
+ |
flucQFrc = tokenizer.nextTokenAsDouble(); |
| 597 |
+ |
sd->setFlucQFrc(flucQFrc); |
| 598 |
+ |
break; |
| 599 |
+ |
} |
| 600 |
+ |
case 'e' : { |
| 601 |
+ |
|
| 602 |
+ |
Vector3d eField; |
| 603 |
+ |
eField[0] = tokenizer.nextTokenAsDouble(); |
| 604 |
+ |
eField[1] = tokenizer.nextTokenAsDouble(); |
| 605 |
+ |
eField[2] = tokenizer.nextTokenAsDouble(); |
| 606 |
+ |
sd->setElectricField(eField); |
| 607 |
+ |
break; |
| 608 |
+ |
} |
| 609 |
+ |
default: { |
| 610 |
+ |
sprintf(painCave.errMsg, |
| 611 |
+ |
"DumpReader Error: %s is an unrecognized type\n", type.c_str()); |
| 612 |
+ |
painCave.isFatal = 1; |
| 613 |
+ |
simError(); |
| 614 |
+ |
break; |
| 615 |
+ |
} |
| 616 |
+ |
} |
| 617 |
+ |
} |
| 618 |
+ |
} |
| 619 |
+ |
|
| 620 |
+ |
|
| 621 |
+ |
void DumpReader::readStuntDoubles(std::istream& inputStream) { |
| 622 |
+ |
|
| 623 |
|
inputStream.getline(buffer, bufferSize); |
| 624 |
|
std::string line(buffer); |
| 625 |
|
|
| 642 |
|
|
| 643 |
|
} |
| 644 |
|
|
| 645 |
+ |
void DumpReader::readSiteData(std::istream& inputStream) { |
| 646 |
+ |
|
| 647 |
+ |
inputStream.getline(buffer, bufferSize); |
| 648 |
+ |
std::string line(buffer); |
| 649 |
+ |
|
| 650 |
+ |
if (line.find("<SiteData>") == std::string::npos) { |
| 651 |
+ |
// site data isn't required for a simulation, so skip |
| 652 |
+ |
return; |
| 653 |
+ |
} |
| 654 |
+ |
|
| 655 |
+ |
while(inputStream.getline(buffer, bufferSize)) { |
| 656 |
+ |
line = buffer; |
| 657 |
+ |
|
| 658 |
+ |
if(line.find("</SiteData>") != std::string::npos) { |
| 659 |
+ |
break; |
| 660 |
+ |
} |
| 661 |
+ |
|
| 662 |
+ |
parseSiteLine(line); |
| 663 |
+ |
} |
| 664 |
+ |
|
| 665 |
+ |
} |
| 666 |
+ |
|
| 667 |
|
void DumpReader::readFrameProperties(std::istream& inputStream) { |
| 668 |
|
|
| 669 |
|
Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot(); |
| 709 |
|
hmat(2, 2) = tokenizer.nextTokenAsDouble(); |
| 710 |
|
s->setHmat(hmat); |
| 711 |
|
} else if (propertyName == "Thermostat") { |
| 712 |
< |
RealType chi = tokenizer.nextTokenAsDouble(); |
| 713 |
< |
RealType integralOfChiDt = tokenizer.nextTokenAsDouble(); |
| 714 |
< |
s->setChi(chi); |
| 715 |
< |
s->setIntegralOfChiDt(integralOfChiDt); |
| 712 |
> |
pair<RealType, RealType> thermostat; |
| 713 |
> |
thermostat.first = tokenizer.nextTokenAsDouble(); |
| 714 |
> |
thermostat.second = tokenizer.nextTokenAsDouble(); |
| 715 |
> |
s->setThermostat(thermostat); |
| 716 |
|
} else if (propertyName == "Barostat") { |
| 717 |
|
Mat3x3d eta; |
| 718 |
|
eta(0, 0) = tokenizer.nextTokenAsDouble(); |
| 724 |
|
eta(2, 0) = tokenizer.nextTokenAsDouble(); |
| 725 |
|
eta(2, 1) = tokenizer.nextTokenAsDouble(); |
| 726 |
|
eta(2, 2) = tokenizer.nextTokenAsDouble(); |
| 727 |
< |
s->setEta(eta); |
| 727 |
> |
s->setBarostat(eta); |
| 728 |
|
} else { |
| 729 |
|
sprintf(painCave.errMsg, |
| 730 |
|
"DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str()); |