| 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). | 
| 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 | + | #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 |  |  | 
| 304 |  | void DumpWriter::writeFrame(std::ostream& os) { | 
| 305 |  |  | 
| 306 |  | #ifdef IS_MPI | 
| 307 | < | MPI::Status istatus; | 
| 307 | > | MPI_Status istatus; | 
| 308 |  | #endif | 
| 309 |  |  | 
| 310 |  | Molecule* mol; | 
| 312 |  | SimInfo::MoleculeIterator mi; | 
| 313 |  | Molecule::IntegrableObjectIterator ii; | 
| 314 |  | RigidBody::AtomIterator ai; | 
| 313 | – | Atom* atom; | 
| 315 |  |  | 
| 316 |  | #ifndef IS_MPI | 
| 317 |  | os << "  <Snapshot>\n"; | 
| 319 |  | writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot()); | 
| 320 |  |  | 
| 321 |  | os << "    <StuntDoubles>\n"; | 
| 322 | < | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { | 
| 323 | < |  | 
| 322 | > | for (mol = info_->beginMolecule(mi); mol != NULL; | 
| 323 | > | mol = info_->nextMolecule(mi)) { | 
| 324 |  |  | 
| 325 |  | for (sd = mol->beginIntegrableObject(ii); sd != NULL; | 
| 326 |  | sd = mol->nextIntegrableObject(ii)) { | 
| 336 |  | mol = info_->nextMolecule(mi)) { | 
| 337 |  |  | 
| 338 |  | for (sd = mol->beginIntegrableObject(ii); sd != NULL; | 
| 339 | < | sd = mol->nextIntegrableObject(ii)) { | 
| 340 | < |  | 
| 339 | > | sd = mol->nextIntegrableObject(ii)) { | 
| 340 | > |  | 
| 341 |  | int ioIndex = sd->getGlobalIntegrableObjectIndex(); | 
| 342 |  | // do one for the IO itself | 
| 343 |  | os << prepareSiteLine(sd, ioIndex, 0); | 
| 346 |  |  | 
| 347 |  | RigidBody* rb = static_cast<RigidBody*>(sd); | 
| 348 |  | int siteIndex = 0; | 
| 349 | < | for (atom = rb->beginAtom(ai); atom != NULL; | 
| 349 | > | for (Atom* atom = rb->beginAtom(ai); atom != NULL; | 
| 350 |  | atom = rb->nextAtom(ai)) { | 
| 351 |  | os << prepareSiteLine(atom, ioIndex, siteIndex); | 
| 352 |  | siteIndex++; | 
| 362 |  | #else | 
| 363 |  |  | 
| 364 |  | const int masterNode = 0; | 
| 365 | < | int worldRank = MPI::COMM_WORLD.Get_rank(); | 
| 366 | < | int nProc = MPI::COMM_WORLD.Get_size(); | 
| 365 | > | int worldRank; | 
| 366 | > | int nProc; | 
| 367 | > |  | 
| 368 | > | MPI_Comm_size( MPI_COMM_WORLD, &nProc); | 
| 369 | > | MPI_Comm_rank( MPI_COMM_WORLD, &worldRank); | 
| 370 |  |  | 
| 371 | + |  | 
| 372 |  | if (worldRank == masterNode) { | 
| 373 |  | os << "  <Snapshot>\n"; | 
| 374 |  | writeFrameProperties(os, | 
| 391 |  |  | 
| 392 |  | for (int i = 1; i < nProc; ++i) { | 
| 393 |  | // tell processor i to start sending us data: | 
| 394 | < | MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode); | 
| 394 | > | MPI_Bcast(&i, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 395 |  |  | 
| 396 |  | // receive the length of the string buffer that was | 
| 397 |  | // prepared by processor i: | 
| 398 |  | int recvLength; | 
| 399 | < | MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG, | 
| 400 | < | istatus); | 
| 399 | > | MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD, | 
| 400 | > | &istatus); | 
| 401 |  |  | 
| 402 |  | // create a buffer to receive the data | 
| 403 |  | char* recvBuffer = new char[recvLength]; | 
| 404 |  | if (recvBuffer == NULL) { | 
| 405 |  | } else { | 
| 406 |  | // receive the data: | 
| 407 | < | MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i, | 
| 408 | < | MPI::ANY_TAG, istatus); | 
| 407 | > | MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, | 
| 408 | > | MPI_ANY_TAG, MPI_COMM_WORLD, &istatus); | 
| 409 |  | // send it to the file: | 
| 410 |  | os << recvBuffer; | 
| 411 |  | // get rid of the receive buffer: | 
| 417 |  | int myturn = 0; | 
| 418 |  | for (int i = 1; i < nProc; ++i){ | 
| 419 |  | // wait for the master node to call our number: | 
| 420 | < | MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode); | 
| 420 | > | MPI_Bcast(&myturn, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 421 |  | if (myturn == worldRank){ | 
| 422 |  | // send the length of our buffer: | 
| 423 | < | MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0); | 
| 423 | > | MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD); | 
| 424 |  |  | 
| 425 |  | // send our buffer: | 
| 426 | < | MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength, | 
| 427 | < | MPI::CHAR, masterNode, 0); | 
| 426 | > | MPI_Send((void *)buffer.c_str(), sendBufferLength, | 
| 427 | > | MPI_CHAR, masterNode, 0, MPI_COMM_WORLD); | 
| 428 | > |  | 
| 429 |  | } | 
| 430 |  | } | 
| 431 |  | } | 
| 453 |  |  | 
| 454 |  | RigidBody* rb = static_cast<RigidBody*>(sd); | 
| 455 |  | int siteIndex = 0; | 
| 456 | < | for (atom = rb->beginAtom(ai); atom != NULL; | 
| 456 | > | for (Atom* atom = rb->beginAtom(ai); atom != NULL; | 
| 457 |  | atom = rb->nextAtom(ai)) { | 
| 458 |  | buffer += prepareSiteLine(atom, ioIndex, siteIndex); | 
| 459 |  | siteIndex++; | 
| 468 |  | for (int i = 1; i < nProc; ++i) { | 
| 469 |  |  | 
| 470 |  | // tell processor i to start sending us data: | 
| 471 | < | MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode); | 
| 471 | > | MPI_Bcast(&i, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 472 |  |  | 
| 473 |  | // receive the length of the string buffer that was | 
| 474 |  | // prepared by processor i: | 
| 475 |  | int recvLength; | 
| 476 | < | MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG, | 
| 477 | < | istatus); | 
| 476 | > | MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD, | 
| 477 | > | &istatus); | 
| 478 |  |  | 
| 479 |  | // create a buffer to receive the data | 
| 480 |  | char* recvBuffer = new char[recvLength]; | 
| 481 |  | if (recvBuffer == NULL) { | 
| 482 |  | } else { | 
| 483 |  | // receive the data: | 
| 484 | < | MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i, | 
| 485 | < | MPI::ANY_TAG, istatus); | 
| 484 | > | MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, | 
| 485 | > | MPI_ANY_TAG, MPI_COMM_WORLD, &istatus); | 
| 486 |  | // send it to the file: | 
| 487 |  | os << recvBuffer; | 
| 488 |  | // get rid of the receive buffer: | 
| 494 |  | int myturn = 0; | 
| 495 |  | for (int i = 1; i < nProc; ++i){ | 
| 496 |  | // wait for the master node to call our number: | 
| 497 | < | MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode); | 
| 497 | > | MPI_Bcast(&myturn, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 498 |  | if (myturn == worldRank){ | 
| 499 |  | // send the length of our buffer: | 
| 500 | < | MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0); | 
| 500 | > | MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD); | 
| 501 |  | // send our buffer: | 
| 502 | < | MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength, | 
| 503 | < | MPI::CHAR, masterNode, 0); | 
| 502 | > | MPI_Send((void *)buffer.c_str(), sendBufferLength, | 
| 503 | > | MPI_CHAR, masterNode, 0, MPI_COMM_WORLD); | 
| 504 |  | } | 
| 505 |  | } | 
| 506 |  | } | 
| 631 |  | } | 
| 632 |  |  | 
| 633 |  | std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) { | 
| 634 | < |  | 
| 634 | > | int storageLayout = info_->getSnapshotManager()->getStorageLayout(); | 
| 635 |  |  | 
| 636 |  | std::string id; | 
| 637 |  | std::string type; | 
| 647 |  | } | 
| 648 |  |  | 
| 649 |  | if (needFlucQ_) { | 
| 650 | < | type += "cw"; | 
| 651 | < | RealType fqPos = sd->getFlucQPos(); | 
| 652 | < | if (isinf(fqPos) || isnan(fqPos) ) { | 
| 653 | < | sprintf( painCave.errMsg, | 
| 654 | < | "DumpWriter detected a numerical error writing the" | 
| 655 | < | " fluctuating charge for object %s", id.c_str()); | 
| 656 | < | painCave.isFatal = 1; | 
| 657 | < | simError(); | 
| 658 | < | } | 
| 659 | < | sprintf(tempBuffer, " %13e ", fqPos); | 
| 660 | < | line += tempBuffer; | 
| 661 | < |  | 
| 662 | < | 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; | 
| 650 | > | if (storageLayout & DataStorage::dslFlucQPosition) { | 
| 651 | > | type += "c"; | 
| 652 | > | RealType fqPos = sd->getFlucQPos(); | 
| 653 | > | if (isinf(fqPos) || isnan(fqPos) ) { | 
| 654 | > | sprintf( painCave.errMsg, | 
| 655 | > | "DumpWriter detected a numerical error writing the" | 
| 656 | > | " fluctuating charge for object %s", id.c_str()); | 
| 657 | > | painCave.isFatal = 1; | 
| 658 | > | simError(); | 
| 659 | > | } | 
| 660 | > | sprintf(tempBuffer, " %13e ", fqPos); | 
| 661 | > | line += tempBuffer; | 
| 662 | > | } | 
| 663 |  |  | 
| 664 | < | if (needForceVector_) { | 
| 665 | < | type += "g"; | 
| 666 | < | RealType fqFrc = sd->getFlucQFrc(); | 
| 667 | < | if (isinf(fqFrc) || isnan(fqFrc) ) { | 
| 664 | > | if (storageLayout & DataStorage::dslFlucQVelocity) { | 
| 665 | > | type += "w"; | 
| 666 | > | RealType fqVel = sd->getFlucQVel(); | 
| 667 | > | if (isinf(fqVel) || isnan(fqVel) ) { | 
| 668 |  | sprintf( painCave.errMsg, | 
| 669 |  | "DumpWriter detected a numerical error writing the" | 
| 670 | < | " fluctuating charge force for object %s", id.c_str()); | 
| 670 | > | " fluctuating charge velocity for object %s", id.c_str()); | 
| 671 |  | painCave.isFatal = 1; | 
| 672 |  | simError(); | 
| 673 |  | } | 
| 674 | < | sprintf(tempBuffer, " %13e ", fqFrc); | 
| 674 | > | sprintf(tempBuffer, " %13e ", fqVel); | 
| 675 |  | line += tempBuffer; | 
| 676 |  | } | 
| 680 | – | } | 
| 677 |  |  | 
| 678 | + | if (needForceVector_) { | 
| 679 | + | if (storageLayout & DataStorage::dslFlucQForce) { | 
| 680 | + | type += "g"; | 
| 681 | + | RealType fqFrc = sd->getFlucQFrc(); | 
| 682 | + | if (isinf(fqFrc) || isnan(fqFrc) ) { | 
| 683 | + | sprintf( painCave.errMsg, | 
| 684 | + | "DumpWriter detected a numerical error writing the" | 
| 685 | + | " fluctuating charge force for object %s", id.c_str()); | 
| 686 | + | painCave.isFatal = 1; | 
| 687 | + | simError(); | 
| 688 | + | } | 
| 689 | + | sprintf(tempBuffer, " %13e ", fqFrc); | 
| 690 | + | line += tempBuffer; | 
| 691 | + | } | 
| 692 | + | } | 
| 693 | + | } | 
| 694 | + |  | 
| 695 |  | if (needElectricField_) { | 
| 696 | < | type += "e"; | 
| 697 | < | Vector3d eField= sd->getElectricField(); | 
| 698 | < | if (isinf(eField[0]) || isnan(eField[0]) || | 
| 699 | < | isinf(eField[1]) || isnan(eField[1]) || | 
| 700 | < | isinf(eField[2]) || isnan(eField[2]) ) { | 
| 701 | < | sprintf( painCave.errMsg, | 
| 702 | < | "DumpWriter detected a numerical error writing the electric" | 
| 703 | < | " field for object %s", id.c_str()); | 
| 704 | < | painCave.isFatal = 1; | 
| 705 | < | simError(); | 
| 696 | > | if (storageLayout & DataStorage::dslElectricField) { | 
| 697 | > | type += "e"; | 
| 698 | > | Vector3d eField= sd->getElectricField(); | 
| 699 | > | if (isinf(eField[0]) || isnan(eField[0]) || | 
| 700 | > | isinf(eField[1]) || isnan(eField[1]) || | 
| 701 | > | isinf(eField[2]) || isnan(eField[2]) ) { | 
| 702 | > | sprintf( painCave.errMsg, | 
| 703 | > | "DumpWriter detected a numerical error writing the electric" | 
| 704 | > | " field for object %s", id.c_str()); | 
| 705 | > | painCave.isFatal = 1; | 
| 706 | > | simError(); | 
| 707 | > | } | 
| 708 | > | sprintf(tempBuffer, " %13e %13e %13e", | 
| 709 | > | eField[0], eField[1], eField[2]); | 
| 710 | > | line += tempBuffer; | 
| 711 |  | } | 
| 694 | – | sprintf(tempBuffer, " %13e %13e %13e", | 
| 695 | – | eField[0], eField[1], eField[2]); | 
| 696 | – | line += tempBuffer; | 
| 712 |  | } | 
| 713 |  |  | 
| 714 |  |  | 
| 715 |  | if (needParticlePot_) { | 
| 716 | < | type += "u"; | 
| 717 | < | RealType particlePot = sd->getParticlePot(); | 
| 718 | < | if (isinf(particlePot) || isnan(particlePot)) { | 
| 719 | < | sprintf( painCave.errMsg, | 
| 720 | < | "DumpWriter detected a numerical error writing the particle " | 
| 721 | < | " potential for object %s", id.c_str()); | 
| 722 | < | painCave.isFatal = 1; | 
| 723 | < | simError(); | 
| 716 | > | if (storageLayout & DataStorage::dslParticlePot) { | 
| 717 | > | type += "u"; | 
| 718 | > | RealType particlePot = sd->getParticlePot(); | 
| 719 | > | if (isinf(particlePot) || isnan(particlePot)) { | 
| 720 | > | sprintf( painCave.errMsg, | 
| 721 | > | "DumpWriter detected a numerical error writing the particle " | 
| 722 | > | " potential for object %s", id.c_str()); | 
| 723 | > | painCave.isFatal = 1; | 
| 724 | > | simError(); | 
| 725 | > | } | 
| 726 | > | sprintf(tempBuffer, " %13e", particlePot); | 
| 727 | > | line += tempBuffer; | 
| 728 |  | } | 
| 710 | – | sprintf(tempBuffer, " %13e", particlePot); | 
| 711 | – | line += tempBuffer; | 
| 729 |  | } | 
| 730 | < |  | 
| 714 | < |  | 
| 730 | > |  | 
| 731 |  | sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str()); | 
| 732 |  | return std::string(tempBuffer); | 
| 733 |  | } | 
| 737 |  | } | 
| 738 |  |  | 
| 739 |  | void DumpWriter::writeEor() { | 
| 740 | < | std::ostream* eorStream; | 
| 741 | < |  | 
| 740 | > |  | 
| 741 | > | std::ostream* eorStream = NULL; | 
| 742 | > |  | 
| 743 |  | #ifdef IS_MPI | 
| 744 |  | if (worldRank == 0) { | 
| 745 |  | #endif // is_mpi | 
| 746 | < |  | 
| 746 | > |  | 
| 747 |  | eorStream = createOStream(eorFilename_); | 
| 748 |  |  | 
| 749 |  | #ifdef IS_MPI | 
| 750 |  | } | 
| 751 | < | #endif // is_mpi | 
| 752 | < |  | 
| 751 | > | #endif | 
| 752 | > |  | 
| 753 |  | writeFrame(*eorStream); | 
| 754 | < |  | 
| 754 | > |  | 
| 755 |  | #ifdef IS_MPI | 
| 756 |  | if (worldRank == 0) { | 
| 757 | < | #endif // is_mpi | 
| 757 | > | #endif | 
| 758 | > |  | 
| 759 |  | writeClosing(*eorStream); | 
| 760 |  | delete eorStream; | 
| 761 | + |  | 
| 762 |  | #ifdef IS_MPI | 
| 763 |  | } | 
| 764 |  | #endif // is_mpi | 
| 772 |  | #ifdef IS_MPI | 
| 773 |  | if (worldRank == 0) { | 
| 774 |  | #endif // is_mpi | 
| 756 | – |  | 
| 775 |  | buffers.push_back(dumpFile_->rdbuf()); | 
| 758 | – |  | 
| 776 |  | eorStream = createOStream(eorFilename_); | 
| 760 | – |  | 
| 777 |  | buffers.push_back(eorStream->rdbuf()); | 
| 762 | – |  | 
| 778 |  | #ifdef IS_MPI | 
| 779 |  | } | 
| 780 |  | #endif // is_mpi | 
| 781 |  |  | 
| 782 |  | TeeBuf tbuf(buffers.begin(), buffers.end()); | 
| 783 |  | std::ostream os(&tbuf); | 
| 769 | – |  | 
| 784 |  | writeFrame(os); | 
| 785 |  |  | 
| 786 |  | #ifdef IS_MPI | 
| 790 |  | delete eorStream; | 
| 791 |  | #ifdef IS_MPI | 
| 792 |  | } | 
| 793 | < | #endif // is_mpi | 
| 780 | < |  | 
| 793 | > | #endif // is_mpi | 
| 794 |  | } | 
| 795 |  |  | 
| 796 |  | std::ostream* DumpWriter::createOStream(const std::string& filename) { |