--- trunk/src/io/DumpWriter.cpp 2008/01/23 21:23:02 1217 +++ branches/development/src/io/DumpWriter.cpp 2012/09/13 14:10:11 1798 @@ -1,24 +1,15 @@ /* - * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. + * Copyright (c) 2009 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a * non-exclusive, royalty free, license to use, modify and * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,27 +28,55 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [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). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include "io/DumpWriter.hpp" #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 -namespace oopse { +using namespace std; +namespace OpenMD { DumpWriter::DumpWriter(SimInfo* info) : info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){ Globals* simParams = info->getSimParams(); - needCompression_ = simParams->getCompressDumpFile(); - needForceVector_ = simParams->getOutputForceVector(); + needCompression_ = simParams->getCompressDumpFile(); + needForceVector_ = simParams->getOutputForceVector(); + needParticlePot_ = simParams->getOutputParticlePotential(); + needFlucQ_ = simParams->getOutputFluctuatingCharges(); + needElectricField_ = simParams->getOutputElectricField(); + + if (needParticlePot_ || needFlucQ_ || needElectricField_) { + doSiteData_ = true; + } else { + doSiteData_ = false; + } + createDumpFile_ = true; #ifdef HAVE_LIBZ if (needCompression_) { @@ -95,8 +114,18 @@ namespace oopse { Globals* simParams = info->getSimParams(); eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor"; - needCompression_ = simParams->getCompressDumpFile(); - needForceVector_ = simParams->getOutputForceVector(); + needCompression_ = simParams->getCompressDumpFile(); + needForceVector_ = simParams->getOutputForceVector(); + needParticlePot_ = simParams->getOutputParticlePotential(); + needFlucQ_ = simParams->getOutputFluctuatingCharges(); + needElectricField_ = simParams->getOutputElectricField(); + + if (needParticlePot_ || needFlucQ_ || needElectricField_) { + doSiteData_ = true; + } else { + doSiteData_ = false; + } + createDumpFile_ = true; #ifdef HAVE_LIBZ if (needCompression_) { @@ -134,9 +163,18 @@ namespace oopse { Globals* simParams = info->getSimParams(); eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor"; - needCompression_ = simParams->getCompressDumpFile(); - needForceVector_ = simParams->getOutputForceVector(); - + needCompression_ = simParams->getCompressDumpFile(); + needForceVector_ = simParams->getOutputForceVector(); + needParticlePot_ = simParams->getOutputParticlePotential(); + needFlucQ_ = simParams->getOutputFluctuatingCharges(); + needElectricField_ = simParams->getOutputElectricField(); + + if (needParticlePot_ || needFlucQ_ || needElectricField_) { + doSiteData_ = true; + } else { + doSiteData_ = false; + } + #ifdef HAVE_LIBZ if (needCompression_) { filename_ += ".gz"; @@ -194,24 +232,64 @@ namespace oopse { os << " \n"; RealType currentTime = s->getTime(); + + if (isinf(currentTime) || isnan(currentTime)) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the time"); + painCave.isFatal = 1; + simError(); + } + sprintf(buffer, " Time: %.10g\n", currentTime); os << buffer; Mat3x3d hmat; hmat = s->getHmat(); + + for (unsigned int i = 0; i < 3; i++) { + for (unsigned int j = 0; j < 3; j++) { + if (isinf(hmat(i,j)) || isnan(hmat(i,j))) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the box"); + painCave.isFatal = 1; + simError(); + } + } + } + sprintf(buffer, " Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n", hmat(0, 0), hmat(1, 0), hmat(2, 0), hmat(0, 1), hmat(1, 1), hmat(2, 1), hmat(0, 2), hmat(1, 2), hmat(2, 2)); os << buffer; - RealType chi = s->getChi(); - RealType integralOfChiDt = s->getIntegralOfChiDt(); - sprintf(buffer, " Thermostat: %.10g , %.10g\n", chi, 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", 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++) { + if (isinf(eta(i,j)) || isnan(eta(i,j))) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the barostat"); + painCave.isFatal = 1; + simError(); + } + } + } + sprintf(buffer, " Barostat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n", eta(0, 0), eta(1, 0), eta(2, 0), eta(0, 1), eta(1, 1), eta(2, 1), @@ -224,13 +302,15 @@ namespace oopse { 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; + Atom* atom; #ifndef IS_MPI os << " \n"; @@ -241,91 +321,265 @@ namespace oopse { 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); } } os << " \n"; - + + if (doSiteData_) { + os << " \n"; + 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 + os << 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)) { + os << prepareSiteLine(atom, ioIndex, siteIndex); + siteIndex++; + } + } + } + } + os << " \n"; + } os << " \n"; 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 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; + } - int nProc; - MPI_Comm_size(MPI_COMM_WORLD, &nProc); + //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 - + // 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; - delete recvBuffer; + // get rid of the receive buffer: + delete [] recvBuffer; } } - os << " \n"; - - os << " \n"; - os.flush(); } else { int sendBufferLength = buffer.size() + 1; - 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); + 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"; + } -#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++; + } + } + } + } + + 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* integrableObject) { + 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(); - vel = integrableObject->getVel(); + pos = sd->getPos(); + + if (isinf(pos[0]) || isnan(pos[0]) || + isinf(pos[1]) || isnan(pos[1]) || + isinf(pos[2]) || isnan(pos[2]) ) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the position" + " for object %d", index); + painCave.isFatal = 1; + simError(); + } + + vel = sd->getVel(); + + if (isinf(vel[0]) || isnan(vel[0]) || + isinf(vel[1]) || isnan(vel[1]) || + isinf(vel[2]) || isnan(vel[2]) ) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the velocity" + " for object %d", index); + painCave.isFatal = 1; + simError(); + } + sprintf(tempBuffer, "%18.10g %18.10g %18.10g %13e %13e %13e", pos[0], pos[1], pos[2], vel[0], vel[1], vel[2]); line += tempBuffer; - if (integrableObject->isDirectional()) { + if (sd->isDirectional()) { type += "qj"; Quat4d q; Vector3d ji; - q = integrableObject->getQ(); - ji = integrableObject->getJ(); + q = sd->getQ(); + + if (isinf(q[0]) || isnan(q[0]) || + isinf(q[1]) || isnan(q[1]) || + isinf(q[2]) || isnan(q[2]) || + isinf(q[3]) || isnan(q[3]) ) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the quaternion" + " for object %d", index); + painCave.isFatal = 1; + simError(); + } + + ji = sd->getJ(); + + if (isinf(ji[0]) || isnan(ji[0]) || + isinf(ji[1]) || isnan(ji[1]) || + isinf(ji[2]) || isnan(ji[2]) ) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the angular" + " momentum for object %d", index); + painCave.isFatal = 1; + simError(); + } + sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e %13e", q[0], q[1], q[2], q[3], ji[0], ji[1], ji[2]); @@ -333,22 +587,135 @@ namespace oopse { } if (needForceVector_) { - type += "ft"; - Vector3d frc; - Vector3d trq; - frc = integrableObject->getFrc(); - trq = integrableObject->getTrq(); - - sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e", - frc[0], frc[1], frc[2], - trq[0], trq[1], trq[2]); + type += "f"; + Vector3d frc = sd->getFrc(); + if (isinf(frc[0]) || isnan(frc[0]) || + isinf(frc[1]) || isnan(frc[1]) || + isinf(frc[2]) || isnan(frc[2]) ) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the force" + " for object %d", index); + painCave.isFatal = 1; + simError(); + } + sprintf(tempBuffer, " %13e %13e %13e", + frc[0], frc[1], frc[2]); line += tempBuffer; + + if (sd->isDirectional()) { + type += "t"; + Vector3d trq = sd->getTrq(); + if (isinf(trq[0]) || isnan(trq[0]) || + isinf(trq[1]) || isnan(trq[1]) || + isinf(trq[2]) || isnan(trq[2]) ) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the torque" + " for object %d", index); + painCave.isFatal = 1; + simError(); + } + sprintf(tempBuffer, " %13e %13e %13e", + trq[0], trq[1], trq[2]); + line += tempBuffer; + } } - + sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str()); return std::string(tempBuffer); } + std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) { + + + std::string id; + std::string type; + std::string line; + char tempBuffer[4096]; + + if (sd->isRigidBody()) { + sprintf(tempBuffer, "%10d ", ioIndex); + id = std::string(tempBuffer); + } else { + sprintf(tempBuffer, "%10d %10d", ioIndex, siteIndex); + id = std::string(tempBuffer); + } + + if (needFlucQ_) { + type += "cw"; + RealType fqPos = sd->getFlucQPos(); + if (isinf(fqPos) || isnan(fqPos) ) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the" + " fluctuating charge for object %s", id.c_str()); + painCave.isFatal = 1; + simError(); + } + sprintf(tempBuffer, " %13e ", fqPos); + line += tempBuffer; + + RealType fqVel = sd->getFlucQVel(); + if (isinf(fqVel) || isnan(fqVel) ) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the" + " fluctuating charge velocity for object %s", id.c_str()); + painCave.isFatal = 1; + simError(); + } + sprintf(tempBuffer, " %13e ", fqVel); + line += tempBuffer; + + if (needForceVector_) { + type += "g"; + RealType fqFrc = sd->getFlucQFrc(); + if (isinf(fqFrc) || isnan(fqFrc) ) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the" + " fluctuating charge force for object %s", id.c_str()); + painCave.isFatal = 1; + simError(); + } + sprintf(tempBuffer, " %13e ", fqFrc); + line += tempBuffer; + } + } + + if (needElectricField_) { + type += "e"; + Vector3d eField= sd->getElectricField(); + if (isinf(eField[0]) || isnan(eField[0]) || + isinf(eField[1]) || isnan(eField[1]) || + isinf(eField[2]) || isnan(eField[2]) ) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the electric" + " field for object %s", id.c_str()); + painCave.isFatal = 1; + simError(); + } + sprintf(tempBuffer, " %13e %13e %13e", + eField[0], eField[1], eField[2]); + line += tempBuffer; + } + + + if (needParticlePot_) { + type += "u"; + RealType particlePot = sd->getParticlePot(); + if (isinf(particlePot) || isnan(particlePot)) { + sprintf( painCave.errMsg, + "DumpWriter detected a numerical error writing the particle " + " potential for object %s", id.c_str()); + painCave.isFatal = 1; + simError(); + } + sprintf(tempBuffer, " %13e", particlePot); + line += tempBuffer; + } + + + sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str()); + return std::string(tempBuffer); + } + void DumpWriter::writeDump() { writeFrame(*dumpFile_); } @@ -416,7 +783,7 @@ namespace oopse { 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 { @@ -426,7 +793,7 @@ namespace oopse { newOStream = new std::ofstream(filename.c_str()); #endif //write out MetaData first - (*newOStream) << "" << std::endl; + (*newOStream) << "" << std::endl; (*newOStream) << " " << std::endl; (*newOStream) << info_->getRawMetaData(); (*newOStream) << " " << std::endl; @@ -435,8 +802,8 @@ namespace oopse { void DumpWriter::writeClosing(std::ostream& os) { - os << "\n"; + os << "\n"; os.flush(); } -}//end namespace oopse +}//end namespace OpenMD