# | Line 35 | Line 35 | |
---|---|---|
35 | * | |
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). |
38 | > | * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
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 | |
44 | #define _FILE_OFFSET_BITS 64 | |
45 | ||
46 | + | #ifdef IS_MPI |
47 | + | #include <mpi.h> |
48 | + | #endif |
49 | + | |
50 | #include <sys/types.h> | |
51 | #include <sys/stat.h> | |
52 | ||
# | Line 57 | Line 62 | |
62 | #include "utils/simError.h" | |
63 | #include "utils/MemoryUtils.hpp" | |
64 | #include "utils/StringTokenizer.hpp" | |
65 | + | #include "brains/Thermo.hpp" |
66 | ||
61 | – | #ifdef IS_MPI |
67 | ||
63 | – | #include <mpi.h> |
64 | – | #define TAKE_THIS_TAG_CHAR 0 |
65 | – | #define TAKE_THIS_TAG_INT 1 |
66 | – | |
67 | – | #endif // is_mpi |
68 | – | |
69 | – | |
68 | namespace OpenMD { | |
69 | ||
70 | DumpReader::DumpReader(SimInfo* info, const std::string& filename) | |
# | Line 77 | Line 75 | namespace OpenMD { | |
75 | if (worldRank == 0) { | |
76 | #endif | |
77 | ||
78 | < | inFile_ = new std::ifstream(filename_.c_str()); |
78 | > | inFile_ = new std::ifstream(filename_.c_str(), |
79 | > | ifstream::in | ifstream::binary); |
80 | ||
81 | if (inFile_->fail()) { | |
82 | sprintf(painCave.errMsg, | |
# | Line 129 | Line 128 | namespace OpenMD { | |
128 | } | |
129 | ||
130 | void DumpReader::scanFile(void) { | |
131 | < | int lineNo = 0; |
131 | > | |
132 | std::streampos prevPos; | |
133 | std::streampos currPos; | |
134 | ||
# | Line 142 | Line 141 | namespace OpenMD { | |
141 | prevPos = currPos; | |
142 | bool foundOpenSnapshotTag = false; | |
143 | bool foundClosedSnapshotTag = false; | |
144 | + | |
145 | + | int lineNo = 0; |
146 | while(inFile_->getline(buffer, bufferSize)) { | |
147 | ++lineNo; | |
148 | ||
# | Line 227 | Line 228 | namespace OpenMD { | |
228 | needVel_ = false; | |
229 | } | |
230 | ||
231 | < | if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) { |
231 | > | if (storageLayout & DataStorage::dslAmat || |
232 | > | storageLayout & DataStorage::dslDipole || |
233 | > | storageLayout & DataStorage::dslQuadrupole) { |
234 | needQuaternion_ = true; | |
235 | } else { | |
236 | needQuaternion_ = false; | |
# | Line 243 | Line 246 | namespace OpenMD { | |
246 | ||
247 | if (needCOMprops_) { | |
248 | Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot(); | |
249 | + | Thermo thermo(info_); |
250 | Vector3d com; | |
251 | < | Vector3d comvel; |
252 | < | Vector3d comw; |
253 | < | if (needPos_ && needVel_){ |
254 | < | info_->getComAll(com, comvel); |
255 | < | comw = info_->getAngularMomentum(); |
256 | < | }else{ |
257 | < | com = info_->getCom(); |
258 | < | comvel = 0.0; |
259 | < | comw = 0.0; |
256 | < | } |
257 | < | s->setCOMprops(com, comvel, comw); |
251 | > | |
252 | > | if (needPos_ && needVel_) { |
253 | > | Vector3d comvel; |
254 | > | Vector3d comw; |
255 | > | thermo.getComAll(com, comvel); |
256 | > | comw = thermo.getAngularMomentum(); |
257 | > | } else { |
258 | > | com = thermo.getCom(); |
259 | > | } |
260 | } | |
259 | – | |
261 | } | |
262 | ||
263 | void DumpReader::readSet(int whichFrame) { | |
# | Line 289 | Line 290 | namespace OpenMD { | |
290 | ||
291 | int sendBufferSize = sendBuffer.size(); | |
292 | MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | |
293 | < | MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
293 | > | MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, |
294 | > | MPI_CHAR, masterNode, MPI_COMM_WORLD); |
295 | ||
296 | sstream.str(sendBuffer); | |
297 | } else { | |
298 | int sendBufferSize; | |
299 | < | MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
299 | > | MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
300 | char * recvBuffer = new char[sendBufferSize+1]; | |
301 | assert(recvBuffer); | |
302 | recvBuffer[sendBufferSize] = '\0'; | |
303 | < | MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
303 | > | MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
304 | sstream.str(recvBuffer); | |
305 | delete [] recvBuffer; | |
306 | } | |
# | Line 324 | Line 326 | namespace OpenMD { | |
326 | ||
327 | inputStream.getline(buffer, bufferSize); | |
328 | line = buffer; | |
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 | < | |
329 | > | |
330 | > | if (line.find("<SiteData>") != std::string::npos) { |
331 | > | //read SiteData |
332 | > | readSiteData(inputStream); |
333 | > | } else { |
334 | > | if (line.find("</Snapshot>") == std::string::npos) { |
335 | > | sprintf(painCave.errMsg, |
336 | > | "DumpReader Error: can not find </Snapshot>\n"); |
337 | > | painCave.isFatal = 1; |
338 | > | simError(); |
339 | > | } |
340 | > | } |
341 | } | |
342 | ||
343 | void DumpReader::parseDumpLine(const std::string& line) { | |
# | Line 350 | Line 357 | namespace OpenMD { | |
357 | ||
358 | int index = tokenizer.nextTokenAsInt(); | |
359 | ||
360 | < | StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index); |
360 | > | StuntDouble* sd = info_->getIOIndexToIntegrableObject(index); |
361 | ||
362 | < | if (integrableObject == NULL) { |
362 | > | if (sd == NULL) { |
363 | return; | |
364 | } | |
365 | std::string type = tokenizer.nextToken(); | |
# | Line 372 | Line 379 | namespace OpenMD { | |
379 | } | |
380 | } | |
381 | ||
382 | < | if (integrableObject->isDirectional()) { |
382 | > | if (sd->isDirectional()) { |
383 | if (needQuaternion_) { | |
384 | found = type.find("q"); | |
385 | if (found == std::string::npos) { | |
# | Line 395 | Line 402 | namespace OpenMD { | |
402 | pos[1] = tokenizer.nextTokenAsDouble(); | |
403 | pos[2] = tokenizer.nextTokenAsDouble(); | |
404 | if (needPos_) { | |
405 | < | integrableObject->setPos(pos); |
405 | > | sd->setPos(pos); |
406 | } | |
407 | break; | |
408 | } | |
# | Line 405 | Line 412 | namespace OpenMD { | |
412 | vel[1] = tokenizer.nextTokenAsDouble(); | |
413 | vel[2] = tokenizer.nextTokenAsDouble(); | |
414 | if (needVel_) { | |
415 | < | integrableObject->setVel(vel); |
415 | > | sd->setVel(vel); |
416 | } | |
417 | break; | |
418 | } | |
419 | ||
420 | case 'q' : { | |
421 | Quat4d q; | |
422 | < | if (integrableObject->isDirectional()) { |
422 | > | if (sd->isDirectional()) { |
423 | ||
424 | q[0] = tokenizer.nextTokenAsDouble(); | |
425 | q[1] = tokenizer.nextTokenAsDouble(); | |
# | Line 431 | Line 438 | namespace OpenMD { | |
438 | ||
439 | q.normalize(); | |
440 | if (needQuaternion_) { | |
441 | < | integrableObject->setQ(q); |
441 | > | sd->setQ(q); |
442 | } | |
443 | } | |
444 | break; | |
445 | } | |
446 | case 'j' : { | |
447 | Vector3d ji; | |
448 | < | if (integrableObject->isDirectional()) { |
448 | > | if (sd->isDirectional()) { |
449 | ji[0] = tokenizer.nextTokenAsDouble(); | |
450 | ji[1] = tokenizer.nextTokenAsDouble(); | |
451 | ji[2] = tokenizer.nextTokenAsDouble(); | |
452 | if (needAngMom_) { | |
453 | < | integrableObject->setJ(ji); |
453 | > | sd->setJ(ji); |
454 | } | |
455 | } | |
456 | break; | |
# | Line 454 | Line 461 | namespace OpenMD { | |
461 | force[0] = tokenizer.nextTokenAsDouble(); | |
462 | force[1] = tokenizer.nextTokenAsDouble(); | |
463 | force[2] = tokenizer.nextTokenAsDouble(); | |
464 | < | integrableObject->setFrc(force); |
464 | > | sd->setFrc(force); |
465 | break; | |
466 | } | |
467 | case 't' : { | |
# | Line 463 | Line 470 | namespace OpenMD { | |
470 | torque[0] = tokenizer.nextTokenAsDouble(); | |
471 | torque[1] = tokenizer.nextTokenAsDouble(); | |
472 | torque[2] = tokenizer.nextTokenAsDouble(); | |
473 | < | integrableObject->setTrq(torque); |
473 | > | sd->setTrq(torque); |
474 | > | break; |
475 | > | } |
476 | > | case 'u' : { |
477 | > | |
478 | > | RealType particlePot; |
479 | > | particlePot = tokenizer.nextTokenAsDouble(); |
480 | > | sd->setParticlePot(particlePot); |
481 | > | break; |
482 | > | } |
483 | > | case 'c' : { |
484 | > | |
485 | > | RealType flucQPos; |
486 | > | flucQPos = tokenizer.nextTokenAsDouble(); |
487 | > | sd->setFlucQPos(flucQPos); |
488 | > | break; |
489 | > | } |
490 | > | case 'w' : { |
491 | > | |
492 | > | RealType flucQVel; |
493 | > | flucQVel = tokenizer.nextTokenAsDouble(); |
494 | > | sd->setFlucQVel(flucQVel); |
495 | > | break; |
496 | > | } |
497 | > | case 'g' : { |
498 | > | |
499 | > | RealType flucQFrc; |
500 | > | flucQFrc = tokenizer.nextTokenAsDouble(); |
501 | > | sd->setFlucQFrc(flucQFrc); |
502 | > | break; |
503 | > | } |
504 | > | case 'e' : { |
505 | > | |
506 | > | Vector3d eField; |
507 | > | eField[0] = tokenizer.nextTokenAsDouble(); |
508 | > | eField[1] = tokenizer.nextTokenAsDouble(); |
509 | > | eField[2] = tokenizer.nextTokenAsDouble(); |
510 | > | sd->setElectricField(eField); |
511 | > | break; |
512 | > | } |
513 | > | case 's' : { |
514 | > | |
515 | > | RealType sPot; |
516 | > | sPot = tokenizer.nextTokenAsDouble(); |
517 | > | sd->setSitePotential(sPot); |
518 | break; | |
519 | } | |
520 | default: { | |
# | Line 476 | Line 527 | namespace OpenMD { | |
527 | ||
528 | } | |
529 | } | |
479 | – | |
530 | } | |
531 | ||
532 | ||
533 | < | void DumpReader::readStuntDoubles(std::istream& inputStream) { |
533 | > | void DumpReader::parseSiteLine(const std::string& line) { |
534 | ||
535 | + | StringTokenizer tokenizer(line); |
536 | + | int nTokens; |
537 | + | |
538 | + | nTokens = tokenizer.countTokens(); |
539 | + | |
540 | + | if (nTokens < 2) { |
541 | + | sprintf(painCave.errMsg, |
542 | + | "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str()); |
543 | + | painCave.isFatal = 1; |
544 | + | simError(); |
545 | + | } |
546 | + | |
547 | + | /** |
548 | + | * The first token is the global integrable object index. |
549 | + | */ |
550 | + | |
551 | + | int index = tokenizer.nextTokenAsInt(); |
552 | + | StuntDouble* sd = info_->getIOIndexToIntegrableObject(index); |
553 | + | if (sd == NULL) { |
554 | + | return; |
555 | + | } |
556 | + | |
557 | + | /** |
558 | + | * Test to see if the next token is an integer or not. If not, |
559 | + | * we've got data on the integrable object itself. If there is an |
560 | + | * integer, we're parsing data for a site on a rigid body. |
561 | + | */ |
562 | + | |
563 | + | std::string indexTest = tokenizer.peekNextToken(); |
564 | + | std::istringstream i(indexTest); |
565 | + | int siteIndex; |
566 | + | if (i >> siteIndex) { |
567 | + | // chew up this token and parse as an int: |
568 | + | siteIndex = tokenizer.nextTokenAsInt(); |
569 | + | |
570 | + | if (sd->isRigidBody()) { |
571 | + | RigidBody* rb = static_cast<RigidBody*>(sd); |
572 | + | sd = rb->getAtoms()[siteIndex]; |
573 | + | } |
574 | + | } |
575 | + | |
576 | + | /** |
577 | + | * The next token contains information on what follows. |
578 | + | */ |
579 | + | std::string type = tokenizer.nextToken(); |
580 | + | int size = type.size(); |
581 | + | |
582 | + | for(int i = 0; i < size; ++i) { |
583 | + | switch(type[i]) { |
584 | + | |
585 | + | case 'u' : { |
586 | + | |
587 | + | RealType particlePot; |
588 | + | particlePot = tokenizer.nextTokenAsDouble(); |
589 | + | sd->setParticlePot(particlePot); |
590 | + | break; |
591 | + | } |
592 | + | case 'c' : { |
593 | + | |
594 | + | RealType flucQPos; |
595 | + | flucQPos = tokenizer.nextTokenAsDouble(); |
596 | + | sd->setFlucQPos(flucQPos); |
597 | + | break; |
598 | + | } |
599 | + | case 'w' : { |
600 | + | |
601 | + | RealType flucQVel; |
602 | + | flucQVel = tokenizer.nextTokenAsDouble(); |
603 | + | sd->setFlucQVel(flucQVel); |
604 | + | break; |
605 | + | } |
606 | + | case 'g' : { |
607 | + | |
608 | + | RealType flucQFrc; |
609 | + | flucQFrc = tokenizer.nextTokenAsDouble(); |
610 | + | sd->setFlucQFrc(flucQFrc); |
611 | + | break; |
612 | + | } |
613 | + | case 'e' : { |
614 | + | |
615 | + | Vector3d eField; |
616 | + | eField[0] = tokenizer.nextTokenAsDouble(); |
617 | + | eField[1] = tokenizer.nextTokenAsDouble(); |
618 | + | eField[2] = tokenizer.nextTokenAsDouble(); |
619 | + | sd->setElectricField(eField); |
620 | + | break; |
621 | + | } |
622 | + | case 's' : { |
623 | + | |
624 | + | RealType sPot; |
625 | + | sPot = tokenizer.nextTokenAsDouble(); |
626 | + | sd->setSitePotential(sPot); |
627 | + | break; |
628 | + | } |
629 | + | default: { |
630 | + | sprintf(painCave.errMsg, |
631 | + | "DumpReader Error: %s is an unrecognized type\n", type.c_str()); |
632 | + | painCave.isFatal = 1; |
633 | + | simError(); |
634 | + | break; |
635 | + | } |
636 | + | } |
637 | + | } |
638 | + | } |
639 | + | |
640 | + | |
641 | + | void DumpReader::readStuntDoubles(std::istream& inputStream) { |
642 | + | |
643 | inputStream.getline(buffer, bufferSize); | |
644 | std::string line(buffer); | |
645 | ||
# | Line 504 | Line 662 | namespace OpenMD { | |
662 | ||
663 | } | |
664 | ||
665 | + | void DumpReader::readSiteData(std::istream& inputStream) { |
666 | + | |
667 | + | std::string line(buffer); |
668 | + | |
669 | + | // We already found the starting <SiteData> tag or we wouldn't be |
670 | + | // here, so just start parsing until we get to the ending |
671 | + | // </SiteData> tag: |
672 | + | |
673 | + | while(inputStream.getline(buffer, bufferSize)) { |
674 | + | line = buffer; |
675 | + | |
676 | + | if(line.find("</SiteData>") != std::string::npos) { |
677 | + | break; |
678 | + | } |
679 | + | |
680 | + | parseSiteLine(line); |
681 | + | } |
682 | + | |
683 | + | } |
684 | + | |
685 | void DumpReader::readFrameProperties(std::istream& inputStream) { | |
686 | ||
687 | Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot(); | |
# | Line 549 | Line 727 | namespace OpenMD { | |
727 | hmat(2, 2) = tokenizer.nextTokenAsDouble(); | |
728 | s->setHmat(hmat); | |
729 | } else if (propertyName == "Thermostat") { | |
730 | < | RealType chi = tokenizer.nextTokenAsDouble(); |
731 | < | RealType integralOfChiDt = tokenizer.nextTokenAsDouble(); |
732 | < | s->setChi(chi); |
733 | < | s->setIntegralOfChiDt(integralOfChiDt); |
730 | > | pair<RealType, RealType> thermostat; |
731 | > | thermostat.first = tokenizer.nextTokenAsDouble(); |
732 | > | thermostat.second = tokenizer.nextTokenAsDouble(); |
733 | > | s->setThermostat(thermostat); |
734 | } else if (propertyName == "Barostat") { | |
735 | Mat3x3d eta; | |
736 | eta(0, 0) = tokenizer.nextTokenAsDouble(); | |
# | Line 564 | Line 742 | namespace OpenMD { | |
742 | eta(2, 0) = tokenizer.nextTokenAsDouble(); | |
743 | eta(2, 1) = tokenizer.nextTokenAsDouble(); | |
744 | eta(2, 2) = tokenizer.nextTokenAsDouble(); | |
745 | < | s->setEta(eta); |
745 | > | s->setBarostat(eta); |
746 | } else { | |
747 | sprintf(painCave.errMsg, | |
748 | "DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str()); |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |