--- branches/development/src/io/DumpWriter.cpp 2012/06/10 14:05:02 1752 +++ branches/development/src/io/DumpWriter.cpp 2013/05/15 15:09:35 1874 @@ -35,7 +35,7 @@ * * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). - * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ @@ -44,13 +44,19 @@ #include "primitives/Molecule.hpp" #include "utils/simError.h" #include "io/basic_teebuf.hpp" +#ifdef HAVE_ZLIB #include "io/gzstream.hpp" +#endif #include "io/Globals.hpp" +#ifdef _MSC_VER +#define isnan(x) _isnan((x)) +#define isinf(x) (!_finite(x) && !_isnan(x)) +#endif #ifdef IS_MPI #include -#endif //is_mpi +#endif using namespace std; namespace OpenMD { @@ -257,20 +263,21 @@ namespace OpenMD { hmat(0, 2), hmat(1, 2), hmat(2, 2)); os << buffer; - RealType chi = s->getChi(); - RealType integralOfChiDt = s->getIntegralOfChiDt(); - if (isinf(chi) || isnan(chi) || - isinf(integralOfChiDt) || isnan(integralOfChiDt)) { + pair thermostat = s->getThermostat(); + + if (isinf(thermostat.first) || isnan(thermostat.first) || + isinf(thermostat.second) || isnan(thermostat.second)) { sprintf( painCave.errMsg, "DumpWriter detected a numerical error writing the thermostat"); painCave.isFatal = 1; simError(); } - sprintf(buffer, " Thermostat: %.10g , %.10g\n", chi, integralOfChiDt); + sprintf(buffer, " Thermostat: %.10g , %.10g\n", thermostat.first, + thermostat.second); os << buffer; Mat3x3d eta; - eta = s->getEta(); + eta = s->getBarostat(); for (unsigned int i = 0; i < 3; i++) { for (unsigned int j = 0; j < 3; j++) { @@ -295,11 +302,11 @@ namespace OpenMD { void DumpWriter::writeFrame(std::ostream& os) { #ifdef IS_MPI - MPI_Status istatus; + MPI::Status istatus; #endif Molecule* mol; - StuntDouble* integrableObject; + StuntDouble* sd; SimInfo::MoleculeIterator mi; Molecule::IntegrableObjectIterator ii; RigidBody::AtomIterator ai; @@ -314,9 +321,9 @@ namespace OpenMD { for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { - for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { - os << prepareDumpLine(integrableObject); + for (sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { + os << prepareDumpLine(sd); } } @@ -324,18 +331,19 @@ namespace OpenMD { if (doSiteData_) { os << " \n"; - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { - for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { + for (sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { - int ioIndex = integrableObject->getGlobalIntegrableObjectIndex(); + int ioIndex = sd->getGlobalIntegrableObjectIndex(); // do one for the IO itself - os << prepareSiteLine(integrableObject, ioIndex, 0); + os << prepareSiteLine(sd, ioIndex, 0); - if (integrableObject->isRigidBody()) { + if (sd->isRigidBody()) { - RigidBody* rb = static_cast(integrableObject); + RigidBody* rb = static_cast(sd); int siteIndex = 0; for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) { @@ -351,73 +359,170 @@ namespace OpenMD { os.flush(); #else - //every node prepares the dump lines for integrable objects belong to itself - std::string buffer; - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { - - for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { - buffer += prepareDumpLine(integrableObject); - } - } - const int masterNode = 0; - int nProc; - MPI_Comm_size(MPI_COMM_WORLD, &nProc); + int worldRank = MPI::COMM_WORLD.Get_rank(); + int nProc = MPI::COMM_WORLD.Get_size(); + if (worldRank == masterNode) { os << " \n"; - writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot()); + writeFrameProperties(os, + info_->getSnapshotManager()->getCurrentSnapshot()); os << " \n"; - - os << buffer; + } + //every node prepares the dump lines for integrable objects belong to itself + std::string buffer; + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for (sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { + buffer += prepareDumpLine(sd); + } + } + + if (worldRank == masterNode) { + os << buffer; + for (int i = 1; i < nProc; ++i) { + // tell processor i to start sending us data: + MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode); // receive the length of the string buffer that was - // prepared by processor i - - MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD); + // prepared by processor i: int recvLength; - MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus); + MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG, + istatus); + + // create a buffer to receive the data char* recvBuffer = new char[recvLength]; if (recvBuffer == NULL) { } else { - MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus); + // receive the data: + MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i, + MPI::ANY_TAG, istatus); + // send it to the file: os << recvBuffer; + // get rid of the receive buffer: delete [] recvBuffer; } } - os << " \n"; - - os << " \n"; - os.flush(); } else { int sendBufferLength = buffer.size() + 1; int myturn = 0; for (int i = 1; i < nProc; ++i){ - MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD); + // wait for the master node to call our number: + MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode); if (myturn == worldRank){ - MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD); - MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD); + // send the length of our buffer: + MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0); + + // send our buffer: + MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength, + MPI::CHAR, masterNode, 0); } } } + + if (worldRank == masterNode) { + os << " \n"; + } -#endif // is_mpi + if (doSiteData_) { + if (worldRank == masterNode) { + os << " \n"; + } + buffer.clear(); + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + + for (sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { + + int ioIndex = sd->getGlobalIntegrableObjectIndex(); + // do one for the IO itself + buffer += prepareSiteLine(sd, ioIndex, 0); - } + if (sd->isRigidBody()) { + + RigidBody* rb = static_cast(sd); + int siteIndex = 0; + for (atom = rb->beginAtom(ai); atom != NULL; + atom = rb->nextAtom(ai)) { + buffer += prepareSiteLine(atom, ioIndex, siteIndex); + siteIndex++; + } + } + } + } - std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) { + if (worldRank == masterNode) { + os << buffer; + + for (int i = 1; i < nProc; ++i) { + + // tell processor i to start sending us data: + MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode); + + // receive the length of the string buffer that was + // prepared by processor i: + int recvLength; + MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG, + istatus); + + // create a buffer to receive the data + char* recvBuffer = new char[recvLength]; + if (recvBuffer == NULL) { + } else { + // receive the data: + MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i, + MPI::ANY_TAG, istatus); + // send it to the file: + os << recvBuffer; + // get rid of the receive buffer: + delete [] recvBuffer; + } + } + } else { + int sendBufferLength = buffer.size() + 1; + int myturn = 0; + for (int i = 1; i < nProc; ++i){ + // wait for the master node to call our number: + MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode); + if (myturn == worldRank){ + // send the length of our buffer: + MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0); + // send our buffer: + MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength, + MPI::CHAR, masterNode, 0); + } + } + } + + if (worldRank == masterNode) { + os << " \n"; + } + } + + if (worldRank == masterNode) { + os << " \n"; + os.flush(); + } + +#endif // is_mpi + + } + + std::string DumpWriter::prepareDumpLine(StuntDouble* sd) { - int index = integrableObject->getGlobalIntegrableObjectIndex(); + int index = sd->getGlobalIntegrableObjectIndex(); std::string type("pv"); std::string line; char tempBuffer[4096]; Vector3d pos; Vector3d vel; - pos = integrableObject->getPos(); + pos = sd->getPos(); if (isinf(pos[0]) || isnan(pos[0]) || isinf(pos[1]) || isnan(pos[1]) || @@ -429,7 +534,7 @@ namespace OpenMD { simError(); } - vel = integrableObject->getVel(); + vel = sd->getVel(); if (isinf(vel[0]) || isnan(vel[0]) || isinf(vel[1]) || isnan(vel[1]) || @@ -446,11 +551,11 @@ namespace OpenMD { vel[0], vel[1], vel[2]); line += tempBuffer; - if (integrableObject->isDirectional()) { + if (sd->isDirectional()) { type += "qj"; Quat4d q; Vector3d ji; - q = integrableObject->getQ(); + q = sd->getQ(); if (isinf(q[0]) || isnan(q[0]) || isinf(q[1]) || isnan(q[1]) || @@ -463,7 +568,7 @@ namespace OpenMD { simError(); } - ji = integrableObject->getJ(); + ji = sd->getJ(); if (isinf(ji[0]) || isnan(ji[0]) || isinf(ji[1]) || isnan(ji[1]) || @@ -483,7 +588,7 @@ namespace OpenMD { if (needForceVector_) { type += "f"; - Vector3d frc = integrableObject->getFrc(); + Vector3d frc = sd->getFrc(); if (isinf(frc[0]) || isnan(frc[0]) || isinf(frc[1]) || isnan(frc[1]) || isinf(frc[2]) || isnan(frc[2]) ) { @@ -497,9 +602,9 @@ namespace OpenMD { frc[0], frc[1], frc[2]); line += tempBuffer; - if (integrableObject->isDirectional()) { + if (sd->isDirectional()) { type += "t"; - Vector3d trq = integrableObject->getTrq(); + Vector3d trq = sd->getTrq(); if (isinf(trq[0]) || isnan(trq[0]) || isinf(trq[1]) || isnan(trq[1]) || isinf(trq[2]) || isnan(trq[2]) ) { @@ -519,7 +624,7 @@ namespace OpenMD { return std::string(tempBuffer); } - std::string DumpWriter::prepareSiteLine(StuntDouble* integrableObject, int ioIndex, int siteIndex) { + std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) { std::string id; @@ -527,7 +632,7 @@ namespace OpenMD { std::string line; char tempBuffer[4096]; - if (integrableObject->isRigidBody()) { + if (sd->isRigidBody()) { sprintf(tempBuffer, "%10d ", ioIndex); id = std::string(tempBuffer); } else { @@ -537,7 +642,7 @@ namespace OpenMD { if (needFlucQ_) { type += "cw"; - RealType fqPos = integrableObject->getFlucQPos(); + RealType fqPos = sd->getFlucQPos(); if (isinf(fqPos) || isnan(fqPos) ) { sprintf( painCave.errMsg, "DumpWriter detected a numerical error writing the" @@ -548,7 +653,7 @@ namespace OpenMD { sprintf(tempBuffer, " %13e ", fqPos); line += tempBuffer; - RealType fqVel = integrableObject->getFlucQVel(); + RealType fqVel = sd->getFlucQVel(); if (isinf(fqVel) || isnan(fqVel) ) { sprintf( painCave.errMsg, "DumpWriter detected a numerical error writing the" @@ -561,7 +666,7 @@ namespace OpenMD { if (needForceVector_) { type += "g"; - RealType fqFrc = integrableObject->getFlucQFrc(); + RealType fqFrc = sd->getFlucQFrc(); if (isinf(fqFrc) || isnan(fqFrc) ) { sprintf( painCave.errMsg, "DumpWriter detected a numerical error writing the" @@ -576,7 +681,7 @@ namespace OpenMD { if (needElectricField_) { type += "e"; - Vector3d eField= integrableObject->getElectricField(); + Vector3d eField= sd->getElectricField(); if (isinf(eField[0]) || isnan(eField[0]) || isinf(eField[1]) || isnan(eField[1]) || isinf(eField[2]) || isnan(eField[2]) ) { @@ -594,7 +699,7 @@ namespace OpenMD { if (needParticlePot_) { type += "u"; - RealType particlePot = integrableObject->getParticlePot(); + RealType particlePot = sd->getParticlePot(); if (isinf(particlePot) || isnan(particlePot)) { sprintf( painCave.errMsg, "DumpWriter detected a numerical error writing the particle " @@ -621,20 +726,18 @@ namespace OpenMD { #ifdef IS_MPI if (worldRank == 0) { #endif // is_mpi - + eorStream = createOStream(eorFilename_); - + writeFrame(*eorStream); + #ifdef IS_MPI } -#endif // is_mpi - - writeFrame(*eorStream); - -#ifdef IS_MPI if (worldRank == 0) { #endif // is_mpi + writeClosing(*eorStream); delete eorStream; + #ifdef IS_MPI } #endif // is_mpi @@ -678,7 +781,7 @@ namespace OpenMD { std::ostream* DumpWriter::createOStream(const std::string& filename) { std::ostream* newOStream; -#ifdef HAVE_LIBZ +#ifdef HAVE_ZLIB if (needCompression_) { newOStream = new ogzstream(filename.c_str()); } else {