| 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 | + | #include "config.h" | 
| 44 | + |  | 
| 45 | + | #ifdef IS_MPI | 
| 46 | + | #include <mpi.h> | 
| 47 | + | #endif | 
| 48 |  |  | 
| 49 |  | #include "io/DumpWriter.hpp" | 
| 50 |  | #include "primitives/Molecule.hpp" | 
| 60 |  | #define isinf(x) (!_finite(x) && !_isnan(x)) | 
| 61 |  | #endif | 
| 62 |  |  | 
| 57 | – | #ifdef IS_MPI | 
| 58 | – | #include <mpi.h> | 
| 59 | – | #endif | 
| 60 | – |  | 
| 63 |  | using namespace std; | 
| 64 |  | namespace OpenMD { | 
| 65 |  |  | 
| 66 |  | DumpWriter::DumpWriter(SimInfo* info) | 
| 67 | < | : info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){ | 
| 67 | > | : info_(info), filename_(info->getDumpFileName()), | 
| 68 | > | eorFilename_(info->getFinalConfigFileName()){ | 
| 69 |  |  | 
| 70 |  | Globals* simParams = info->getSimParams(); | 
| 71 |  | needCompression_   = simParams->getCompressDumpFile(); | 
| 73 |  | needParticlePot_   = simParams->getOutputParticlePotential(); | 
| 74 |  | needFlucQ_         = simParams->getOutputFluctuatingCharges(); | 
| 75 |  | needElectricField_ = simParams->getOutputElectricField(); | 
| 76 | + | needSitePotential_ = simParams->getOutputSitePotential(); | 
| 77 |  |  | 
| 78 | < | if (needParticlePot_ || needFlucQ_ || needElectricField_) { | 
| 78 | > | if (needParticlePot_ || needFlucQ_ || needElectricField_ || | 
| 79 | > | needSitePotential_) { | 
| 80 |  | doSiteData_ = true; | 
| 81 |  | } else { | 
| 82 |  | doSiteData_ = false; | 
| 124 |  | needParticlePot_   = simParams->getOutputParticlePotential(); | 
| 125 |  | needFlucQ_         = simParams->getOutputFluctuatingCharges(); | 
| 126 |  | needElectricField_ = simParams->getOutputElectricField(); | 
| 127 | + | needSitePotential_ = simParams->getOutputSitePotential(); | 
| 128 |  |  | 
| 129 | < | if (needParticlePot_ || needFlucQ_ || needElectricField_) { | 
| 129 | > | if (needParticlePot_ || needFlucQ_ || needElectricField_ || | 
| 130 | > | needSitePotential_) { | 
| 131 |  | doSiteData_ = true; | 
| 132 |  | } else { | 
| 133 |  | doSiteData_ = false; | 
| 175 |  | needParticlePot_   = simParams->getOutputParticlePotential(); | 
| 176 |  | needFlucQ_         = simParams->getOutputFluctuatingCharges(); | 
| 177 |  | needElectricField_ = simParams->getOutputElectricField(); | 
| 178 | + | needSitePotential_ = simParams->getOutputSitePotential(); | 
| 179 |  |  | 
| 180 | < | if (needParticlePot_ || needFlucQ_ || needElectricField_) { | 
| 180 | > | if (needParticlePot_ || needFlucQ_ || needElectricField_ || | 
| 181 | > | needSitePotential_) { | 
| 182 |  | doSiteData_ = true; | 
| 183 |  | } else { | 
| 184 |  | doSiteData_ = false; | 
| 311 |  | void DumpWriter::writeFrame(std::ostream& os) { | 
| 312 |  |  | 
| 313 |  | #ifdef IS_MPI | 
| 314 | < | MPI::Status istatus; | 
| 314 | > | MPI_Status istatus; | 
| 315 |  | #endif | 
| 316 |  |  | 
| 317 |  | Molecule* mol; | 
| 343 |  | mol = info_->nextMolecule(mi)) { | 
| 344 |  |  | 
| 345 |  | for (sd = mol->beginIntegrableObject(ii); sd != NULL; | 
| 346 | < | sd = mol->nextIntegrableObject(ii)) { | 
| 347 | < |  | 
| 346 | > | sd = mol->nextIntegrableObject(ii)) { | 
| 347 | > |  | 
| 348 |  | int ioIndex = sd->getGlobalIntegrableObjectIndex(); | 
| 349 |  | // do one for the IO itself | 
| 350 |  | os << prepareSiteLine(sd, ioIndex, 0); | 
| 369 |  | #else | 
| 370 |  |  | 
| 371 |  | const int masterNode = 0; | 
| 372 | < | int worldRank = MPI::COMM_WORLD.Get_rank(); | 
| 373 | < | int nProc = MPI::COMM_WORLD.Get_size(); | 
| 372 | > | int worldRank; | 
| 373 | > | int nProc; | 
| 374 |  |  | 
| 375 | + | MPI_Comm_size( MPI_COMM_WORLD, &nProc); | 
| 376 | + | MPI_Comm_rank( MPI_COMM_WORLD, &worldRank); | 
| 377 | + |  | 
| 378 | + |  | 
| 379 |  | if (worldRank == masterNode) { | 
| 380 |  | os << "  <Snapshot>\n"; | 
| 381 |  | writeFrameProperties(os, | 
| 398 |  |  | 
| 399 |  | for (int i = 1; i < nProc; ++i) { | 
| 400 |  | // tell processor i to start sending us data: | 
| 388 | – | MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode); | 
| 401 |  |  | 
| 402 | + | MPI_Bcast(&i, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 403 | + |  | 
| 404 |  | // receive the length of the string buffer that was | 
| 405 |  | // prepared by processor i: | 
| 406 |  | int recvLength; | 
| 407 | < | MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG, | 
| 408 | < | istatus); | 
| 407 | > | MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD, | 
| 408 | > | &istatus); | 
| 409 |  |  | 
| 410 |  | // create a buffer to receive the data | 
| 411 |  | char* recvBuffer = new char[recvLength]; | 
| 412 |  | if (recvBuffer == NULL) { | 
| 413 |  | } else { | 
| 414 |  | // receive the data: | 
| 415 | < | MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i, | 
| 416 | < | MPI::ANY_TAG, istatus); | 
| 415 | > | MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, | 
| 416 | > | MPI_ANY_TAG, MPI_COMM_WORLD, &istatus); | 
| 417 |  | // send it to the file: | 
| 418 |  | os << recvBuffer; | 
| 419 |  | // get rid of the receive buffer: | 
| 425 |  | int myturn = 0; | 
| 426 |  | for (int i = 1; i < nProc; ++i){ | 
| 427 |  | // wait for the master node to call our number: | 
| 428 | < | MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode); | 
| 428 | > | MPI_Bcast(&myturn, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 429 |  | if (myturn == worldRank){ | 
| 430 |  | // send the length of our buffer: | 
| 417 | – | MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0); | 
| 431 |  |  | 
| 432 | + | MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD); | 
| 433 | + |  | 
| 434 |  | // send our buffer: | 
| 435 | < | MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength, | 
| 436 | < | MPI::CHAR, masterNode, 0); | 
| 435 | > | MPI_Send((void *)buffer.c_str(), sendBufferLength, | 
| 436 | > | MPI_CHAR, masterNode, 0, MPI_COMM_WORLD); | 
| 437 |  |  | 
| 438 |  | } | 
| 439 |  | } | 
| 477 |  | for (int i = 1; i < nProc; ++i) { | 
| 478 |  |  | 
| 479 |  | // tell processor i to start sending us data: | 
| 480 | < | MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode); | 
| 480 | > | MPI_Bcast(&i, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 481 |  |  | 
| 482 |  | // receive the length of the string buffer that was | 
| 483 |  | // prepared by processor i: | 
| 484 |  | int recvLength; | 
| 485 | < | MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG, | 
| 486 | < | istatus); | 
| 485 | > | MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD, | 
| 486 | > | &istatus); | 
| 487 |  |  | 
| 488 |  | // create a buffer to receive the data | 
| 489 |  | char* recvBuffer = new char[recvLength]; | 
| 490 |  | if (recvBuffer == NULL) { | 
| 491 |  | } else { | 
| 492 |  | // receive the data: | 
| 493 | < | MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i, | 
| 494 | < | MPI::ANY_TAG, istatus); | 
| 493 | > | MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, | 
| 494 | > | MPI_ANY_TAG, MPI_COMM_WORLD, &istatus); | 
| 495 |  | // send it to the file: | 
| 496 |  | os << recvBuffer; | 
| 497 |  | // get rid of the receive buffer: | 
| 503 |  | int myturn = 0; | 
| 504 |  | for (int i = 1; i < nProc; ++i){ | 
| 505 |  | // wait for the master node to call our number: | 
| 506 | < | MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode); | 
| 506 | > | MPI_Bcast(&myturn, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 507 |  | if (myturn == worldRank){ | 
| 508 |  | // send the length of our buffer: | 
| 509 | < | MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0); | 
| 509 | > | MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD); | 
| 510 |  | // send our buffer: | 
| 511 | < | MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength, | 
| 512 | < | MPI::CHAR, masterNode, 0); | 
| 511 | > | MPI_Send((void *)buffer.c_str(), sendBufferLength, | 
| 512 | > | MPI_CHAR, masterNode, 0, MPI_COMM_WORLD); | 
| 513 |  | } | 
| 514 |  | } | 
| 515 |  | } | 
| 640 |  | } | 
| 641 |  |  | 
| 642 |  | std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) { | 
| 643 | < |  | 
| 643 | > | int storageLayout = info_->getSnapshotManager()->getStorageLayout(); | 
| 644 |  |  | 
| 645 |  | std::string id; | 
| 646 |  | std::string type; | 
| 656 |  | } | 
| 657 |  |  | 
| 658 |  | if (needFlucQ_) { | 
| 659 | < | type += "cw"; | 
| 660 | < | RealType fqPos = sd->getFlucQPos(); | 
| 661 | < | if (isinf(fqPos) || isnan(fqPos) ) { | 
| 662 | < | sprintf( painCave.errMsg, | 
| 663 | < | "DumpWriter detected a numerical error writing the" | 
| 664 | < | " fluctuating charge for object %s", id.c_str()); | 
| 665 | < | painCave.isFatal = 1; | 
| 666 | < | simError(); | 
| 667 | < | } | 
| 668 | < | sprintf(tempBuffer, " %13e ", fqPos); | 
| 669 | < | line += tempBuffer; | 
| 670 | < |  | 
| 671 | < | RealType fqVel = sd->getFlucQVel(); | 
| 657 | < | if (isinf(fqVel) || isnan(fqVel) ) { | 
| 658 | < | sprintf( painCave.errMsg, | 
| 659 | < | "DumpWriter detected a numerical error writing the" | 
| 660 | < | " fluctuating charge velocity for object %s", id.c_str()); | 
| 661 | < | painCave.isFatal = 1; | 
| 662 | < | simError(); | 
| 663 | < | } | 
| 664 | < | sprintf(tempBuffer, " %13e ", fqVel); | 
| 665 | < | line += tempBuffer; | 
| 659 | > | if (storageLayout & DataStorage::dslFlucQPosition) { | 
| 660 | > | type += "c"; | 
| 661 | > | RealType fqPos = sd->getFlucQPos(); | 
| 662 | > | if (isinf(fqPos) || isnan(fqPos) ) { | 
| 663 | > | sprintf( painCave.errMsg, | 
| 664 | > | "DumpWriter detected a numerical error writing the" | 
| 665 | > | " fluctuating charge for object %s", id.c_str()); | 
| 666 | > | painCave.isFatal = 1; | 
| 667 | > | simError(); | 
| 668 | > | } | 
| 669 | > | sprintf(tempBuffer, " %13e ", fqPos); | 
| 670 | > | line += tempBuffer; | 
| 671 | > | } | 
| 672 |  |  | 
| 673 | < | if (needForceVector_) { | 
| 674 | < | type += "g"; | 
| 675 | < | RealType fqFrc = sd->getFlucQFrc(); | 
| 676 | < | if (isinf(fqFrc) || isnan(fqFrc) ) { | 
| 673 | > | if (storageLayout & DataStorage::dslFlucQVelocity) { | 
| 674 | > | type += "w"; | 
| 675 | > | RealType fqVel = sd->getFlucQVel(); | 
| 676 | > | if (isinf(fqVel) || isnan(fqVel) ) { | 
| 677 |  | sprintf( painCave.errMsg, | 
| 678 |  | "DumpWriter detected a numerical error writing the" | 
| 679 | < | " fluctuating charge force for object %s", id.c_str()); | 
| 679 | > | " fluctuating charge velocity for object %s", id.c_str()); | 
| 680 |  | painCave.isFatal = 1; | 
| 681 |  | simError(); | 
| 682 |  | } | 
| 683 | < | sprintf(tempBuffer, " %13e ", fqFrc); | 
| 683 | > | sprintf(tempBuffer, " %13e ", fqVel); | 
| 684 |  | line += tempBuffer; | 
| 685 |  | } | 
| 680 | – | } | 
| 686 |  |  | 
| 687 | + | if (needForceVector_) { | 
| 688 | + | if (storageLayout & DataStorage::dslFlucQForce) { | 
| 689 | + | type += "g"; | 
| 690 | + | RealType fqFrc = sd->getFlucQFrc(); | 
| 691 | + | if (isinf(fqFrc) || isnan(fqFrc) ) { | 
| 692 | + | sprintf( painCave.errMsg, | 
| 693 | + | "DumpWriter detected a numerical error writing the" | 
| 694 | + | " fluctuating charge force for object %s", id.c_str()); | 
| 695 | + | painCave.isFatal = 1; | 
| 696 | + | simError(); | 
| 697 | + | } | 
| 698 | + | sprintf(tempBuffer, " %13e ", fqFrc); | 
| 699 | + | line += tempBuffer; | 
| 700 | + | } | 
| 701 | + | } | 
| 702 | + | } | 
| 703 | + |  | 
| 704 |  | if (needElectricField_) { | 
| 705 | < | type += "e"; | 
| 706 | < | Vector3d eField= sd->getElectricField(); | 
| 707 | < | if (isinf(eField[0]) || isnan(eField[0]) || | 
| 708 | < | isinf(eField[1]) || isnan(eField[1]) || | 
| 709 | < | isinf(eField[2]) || isnan(eField[2]) ) { | 
| 710 | < | sprintf( painCave.errMsg, | 
| 711 | < | "DumpWriter detected a numerical error writing the electric" | 
| 712 | < | " field for object %s", id.c_str()); | 
| 713 | < | painCave.isFatal = 1; | 
| 714 | < | simError(); | 
| 705 | > | if (storageLayout & DataStorage::dslElectricField) { | 
| 706 | > | type += "e"; | 
| 707 | > | Vector3d eField= sd->getElectricField(); | 
| 708 | > | if (isinf(eField[0]) || isnan(eField[0]) || | 
| 709 | > | isinf(eField[1]) || isnan(eField[1]) || | 
| 710 | > | isinf(eField[2]) || isnan(eField[2]) ) { | 
| 711 | > | sprintf( painCave.errMsg, | 
| 712 | > | "DumpWriter detected a numerical error writing the electric" | 
| 713 | > | " field for object %s", id.c_str()); | 
| 714 | > | painCave.isFatal = 1; | 
| 715 | > | simError(); | 
| 716 | > | } | 
| 717 | > | sprintf(tempBuffer, " %13e %13e %13e", | 
| 718 | > | eField[0], eField[1], eField[2]); | 
| 719 | > | line += tempBuffer; | 
| 720 |  | } | 
| 694 | – | sprintf(tempBuffer, " %13e %13e %13e", | 
| 695 | – | eField[0], eField[1], eField[2]); | 
| 696 | – | line += tempBuffer; | 
| 721 |  | } | 
| 722 |  |  | 
| 723 | < |  | 
| 723 | > | if (needSitePotential_) { | 
| 724 | > | if (storageLayout & DataStorage::dslSitePotential) { | 
| 725 | > | type += "s"; | 
| 726 | > | RealType sPot = sd->getSitePotential(); | 
| 727 | > | if (isinf(sPot) || isnan(sPot) ) { | 
| 728 | > | sprintf( painCave.errMsg, | 
| 729 | > | "DumpWriter detected a numerical error writing the" | 
| 730 | > | " site potential for object %s", id.c_str()); | 
| 731 | > | painCave.isFatal = 1; | 
| 732 | > | simError(); | 
| 733 | > | } | 
| 734 | > | sprintf(tempBuffer, " %13e ", sPot); | 
| 735 | > | line += tempBuffer; | 
| 736 | > | } | 
| 737 | > | } | 
| 738 | > |  | 
| 739 |  | if (needParticlePot_) { | 
| 740 | < | type += "u"; | 
| 741 | < | RealType particlePot = sd->getParticlePot(); | 
| 742 | < | if (isinf(particlePot) || isnan(particlePot)) { | 
| 743 | < | sprintf( painCave.errMsg, | 
| 744 | < | "DumpWriter detected a numerical error writing the particle " | 
| 745 | < | " potential for object %s", id.c_str()); | 
| 746 | < | painCave.isFatal = 1; | 
| 747 | < | simError(); | 
| 740 | > | if (storageLayout & DataStorage::dslParticlePot) { | 
| 741 | > | type += "u"; | 
| 742 | > | RealType particlePot = sd->getParticlePot(); | 
| 743 | > | if (isinf(particlePot) || isnan(particlePot)) { | 
| 744 | > | sprintf( painCave.errMsg, | 
| 745 | > | "DumpWriter detected a numerical error writing the particle " | 
| 746 | > | " potential for object %s", id.c_str()); | 
| 747 | > | painCave.isFatal = 1; | 
| 748 | > | simError(); | 
| 749 | > | } | 
| 750 | > | sprintf(tempBuffer, " %13e", particlePot); | 
| 751 | > | line += tempBuffer; | 
| 752 |  | } | 
| 710 | – | sprintf(tempBuffer, " %13e", particlePot); | 
| 711 | – | line += tempBuffer; | 
| 753 |  | } | 
| 754 | < |  | 
| 714 | < |  | 
| 754 | > |  | 
| 755 |  | sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str()); | 
| 756 |  | return std::string(tempBuffer); | 
| 757 |  | } | 
| 761 |  | } | 
| 762 |  |  | 
| 763 |  | void DumpWriter::writeEor() { | 
| 764 | < | std::ostream* eorStream; | 
| 765 | < |  | 
| 764 | > |  | 
| 765 | > | std::ostream* eorStream = NULL; | 
| 766 | > |  | 
| 767 |  | #ifdef IS_MPI | 
| 768 |  | if (worldRank == 0) { | 
| 769 |  | #endif // is_mpi | 
| 792 |  |  | 
| 793 |  | void DumpWriter::writeDumpAndEor() { | 
| 794 |  | std::vector<std::streambuf*> buffers; | 
| 795 | < | std::ostream* eorStream; | 
| 795 | > | std::ostream* eorStream = NULL; | 
| 796 |  | #ifdef IS_MPI | 
| 797 |  | if (worldRank == 0) { | 
| 798 |  | #endif // is_mpi | 
| 758 | – |  | 
| 799 |  | buffers.push_back(dumpFile_->rdbuf()); | 
| 760 | – |  | 
| 800 |  | eorStream = createOStream(eorFilename_); | 
| 762 | – |  | 
| 801 |  | buffers.push_back(eorStream->rdbuf()); | 
| 764 | – |  | 
| 802 |  | #ifdef IS_MPI | 
| 803 |  | } | 
| 804 |  | #endif // is_mpi | 
| 805 |  |  | 
| 806 |  | TeeBuf tbuf(buffers.begin(), buffers.end()); | 
| 807 |  | std::ostream os(&tbuf); | 
| 771 | – |  | 
| 808 |  | writeFrame(os); | 
| 809 |  |  | 
| 810 |  | #ifdef IS_MPI | 
| 814 |  | delete eorStream; | 
| 815 |  | #ifdef IS_MPI | 
| 816 |  | } | 
| 817 | < | #endif // is_mpi | 
| 782 | < |  | 
| 817 | > | #endif // is_mpi | 
| 818 |  | } | 
| 819 |  |  | 
| 820 |  | std::ostream* DumpWriter::createOStream(const std::string& filename) { |