| 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 |
| 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) { |