--- trunk/src/io/DumpWriter.cpp 2005/02/13 19:10:25 324 +++ branches/development/src/io/DumpWriter.cpp 2013/06/13 14:26:09 1878 @@ -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,516 +28,784 @@ * 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, 234107 (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, const std::string& filename) - : info_(info), filename_(filename){ + DumpWriter::DumpWriter(SimInfo* info) + : info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){ + + Globals* simParams = info->getSimParams(); + 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_) { + filename_ += ".gz"; + eorFilename_ += ".gz"; + } +#endif + #ifdef IS_MPI if (worldRank == 0) { #endif // is_mpi + + dumpFile_ = createOStream(filename_); - dumpFile_.open(filename_.c_str(), std::ios::out | std::ios::trunc); + if (!dumpFile_) { + sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n", + filename_.c_str()); + painCave.isFatal = 1; + simError(); + } - if (!dumpFile_) { - sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n", - filename_.c_str()); - painCave.isFatal = 1; - simError(); - } - #ifdef IS_MPI } - sprintf(checkPointMsg, "Sucessfully opened output file for dumping.\n"); - MPIcheckPoint(); - #endif // is_mpi -} + } -DumpWriter::~DumpWriter() { + DumpWriter::DumpWriter(SimInfo* info, const std::string& filename) + : info_(info), filename_(filename){ + + Globals* simParams = info->getSimParams(); + eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor"; + + 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_) { + filename_ += ".gz"; + eorFilename_ += ".gz"; + } +#endif + #ifdef IS_MPI if (worldRank == 0) { #endif // is_mpi - dumpFile_.close(); + + dumpFile_ = createOStream(filename_); + if (!dumpFile_) { + sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n", + filename_.c_str()); + painCave.isFatal = 1; + simError(); + } + #ifdef IS_MPI } #endif // is_mpi -} + } + + DumpWriter::DumpWriter(SimInfo* info, const std::string& filename, bool writeDumpFile) + : info_(info), filename_(filename){ + + Globals* simParams = info->getSimParams(); + eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor"; + + needCompression_ = simParams->getCompressDumpFile(); + needForceVector_ = simParams->getOutputForceVector(); + needParticlePot_ = simParams->getOutputParticlePotential(); + needFlucQ_ = simParams->getOutputFluctuatingCharges(); + needElectricField_ = simParams->getOutputElectricField(); -void DumpWriter::writeCommentLine(std::ostream& os, Snapshot* s) { + if (needParticlePot_ || needFlucQ_ || needElectricField_) { + doSiteData_ = true; + } else { + doSiteData_ = false; + } - double currentTime; - Mat3x3d hmat; - double chi; - double integralOfChiDt; - Mat3x3d eta; +#ifdef HAVE_LIBZ + if (needCompression_) { + filename_ += ".gz"; + eorFilename_ += ".gz"; + } +#endif - currentTime = s->getTime(); - hmat = s->getHmat(); - chi = s->getChi(); - integralOfChiDt = s->getIntegralOfChiDt(); - eta = s->getEta(); +#ifdef IS_MPI - os << currentTime << ";\t" - << hmat(0, 0) << "\t" << hmat(1, 0) << "\t" << hmat(2, 0) << ";\t" - << hmat(0, 1) << "\t" << hmat(1, 1) << "\t" << hmat(2, 1) << ";\t" - << hmat(0, 2) << "\t" << hmat(1, 2) << "\t" << hmat(2, 2) << ";\t"; + if (worldRank == 0) { +#endif // is_mpi + + createDumpFile_ = writeDumpFile; + if (createDumpFile_) { + dumpFile_ = createOStream(filename_); + + if (!dumpFile_) { + sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n", + filename_.c_str()); + painCave.isFatal = 1; + simError(); + } + } +#ifdef IS_MPI + + } - //write out additional parameters, such as chi and eta + +#endif // is_mpi + + } - os << chi << "\t" << integralOfChiDt << "\t;"; + DumpWriter::~DumpWriter() { - os << eta(0, 0) << "\t" << eta(1, 0) << "\t" << eta(2, 0) << ";\t" - << eta(0, 1) << "\t" << eta(1, 1) << "\t" << eta(2, 1) << ";\t" - << eta(0, 2) << "\t" << eta(1, 2) << "\t" << eta(2, 2) << ";"; - - os << std::endl; -} +#ifdef IS_MPI -void DumpWriter::writeFrame(std::ostream& os) { - const int BUFFERSIZE = 2000; - const int MINIBUFFERSIZE = 100; + if (worldRank == 0) { +#endif // is_mpi + if (createDumpFile_){ + writeClosing(*dumpFile_); + delete dumpFile_; + } +#ifdef IS_MPI - char tempBuffer[BUFFERSIZE]; - char writeLine[BUFFERSIZE]; + } - Quat4d q; - Vector3d ji; - Vector3d pos; - Vector3d vel; +#endif // is_mpi - Molecule* mol; - StuntDouble* integrableObject; - SimInfo::MoleculeIterator mi; - Molecule::IntegrableObjectIterator ii; - - int nTotObjects; - nTotObjects = info_->getNGlobalIntegrableObjects(); + } -#ifndef IS_MPI + void DumpWriter::writeFrameProperties(std::ostream& os, Snapshot* s) { + char buffer[1024]; - os << nTotObjects << "\n"; - - writeCommentLine(os, info_->getSnapshotManager()->getCurrentSnapshot()); + os << " \n"; - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { + RealType currentTime = s->getTime(); - for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { - + 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; - pos = integrableObject->getPos(); - vel = integrableObject->getVel(); + Mat3x3d hmat; + hmat = s->getHmat(); - sprintf(tempBuffer, "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t", - integrableObject->getType().c_str(), - pos[0], pos[1], pos[2], - vel[0], vel[1], vel[2]); + 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; - strcpy(writeLine, tempBuffer); + pair thermostat = s->getThermostat(); - if (integrableObject->isDirectional()) { - q = integrableObject->getQ(); - ji = integrableObject->getJ(); + 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; - sprintf(tempBuffer, "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", - q[0], q[1], q[2], q[3], - ji[0], ji[1], ji[2]); - strcat(writeLine, tempBuffer); - } else { - strcat(writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n"); - } + Mat3x3d eta; + eta = s->getBarostat(); - os << writeLine; - - } + 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(); + } + } } - os.flush(); -#else // is_mpi - /********************************************************************* - * Documentation? You want DOCUMENTATION? - * - * Why all the potatoes below? - * - * To make a long story short, the original version of DumpWriter - * worked in the most inefficient way possible. Node 0 would - * poke each of the node for an individual atom's formatted data - * as node 0 worked its way down the global index. This was particularly - * inefficient since the method blocked all processors at every atom - * (and did it twice!). - * - * An intermediate version of DumpWriter could be described from Node - * zero's perspective as follows: - * - * 1) Have 100 of your friends stand in a circle. - * 2) When you say go, have all of them start tossing potatoes at - * you (one at a time). - * 3) Catch the potatoes. - * - * It was an improvement, but MPI has buffers and caches that could - * best be described in this analogy as "potato nets", so there's no - * need to block the processors atom-by-atom. - * - * This new and improved DumpWriter works in an even more efficient - * way: - * - * 1) Have 100 of your friend stand in a circle. - * 2) When you say go, have them start tossing 5-pound bags of - * potatoes at you. - * 3) Once you've caught a friend's bag of potatoes, - * toss them a spud to let them know they can toss another bag. - * - * How's THAT for documentation? - * - *********************************************************************/ - const int masterNode = 0; + 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), + eta(0, 2), eta(1, 2), eta(2, 2)); + os << buffer; - int * potatoes; - int myPotato; - int nProc; - int which_node; - double atomData[13]; - int isDirectional; - const char * atomTypeString; - char MPIatomTypeString[MINIBUFFERSIZE]; - int msgLen; // the length of message actually recieved at master nodes - int haveError; - MPI_Status istatus; - int nCurObj; - - // code to find maximum tag value - int * tagub; - int flag; - int MAXTAG; - MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, &tagub, &flag); + os << " \n"; + } - if (flag) { - MAXTAG = *tagub; - } else { - MAXTAG = 32767; - } + void DumpWriter::writeFrame(std::ostream& os) { - if (worldRank == masterNode) { //master node (node 0) is responsible for writing the dump file +#ifdef IS_MPI + MPI::Status istatus; +#endif - // Node 0 needs a list of the magic potatoes for each processor; + Molecule* mol; + StuntDouble* sd; + SimInfo::MoleculeIterator mi; + Molecule::IntegrableObjectIterator ii; + RigidBody::AtomIterator ai; - MPI_Comm_size(MPI_COMM_WORLD, &nProc); - potatoes = new int[nProc]; +#ifndef IS_MPI + os << " \n"; + + writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot()); - //write out the comment lines - for(int i = 0; i < nProc; i++) { - potatoes[i] = 0; - } + os << " \n"; + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + + 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)) { - os << nTotObjects << "\n"; - writeCommentLine(os, info_->getSnapshotManager()->getCurrentSnapshot()); + int ioIndex = sd->getGlobalIntegrableObjectIndex(); + // do one for the IO itself + os << prepareSiteLine(sd, ioIndex, 0); - for(int i = 0; i < info_->getNGlobalMolecules(); i++) { + if (sd->isRigidBody()) { + + RigidBody* rb = static_cast(sd); + int siteIndex = 0; + for (Atom* atom = rb->beginAtom(ai); atom != NULL; + atom = rb->nextAtom(ai)) { + os << prepareSiteLine(atom, ioIndex, siteIndex); + siteIndex++; + } + } + } + } + os << " \n"; + } + os << " \n"; - // Get the Node number which has this atom; + os.flush(); +#else - which_node = info_->getMolToProc(i); + const int masterNode = 0; + int worldRank = MPI::COMM_WORLD.Get_rank(); + int nProc = MPI::COMM_WORLD.Get_size(); - if (which_node != masterNode) { //current molecule is in slave node - if (potatoes[which_node] + 1 >= MAXTAG) { - // The potato was going to exceed the maximum value, - // so wrap this processor potato back to 0: + if (worldRank == masterNode) { + os << " \n"; + writeFrameProperties(os, + info_->getSnapshotManager()->getCurrentSnapshot()); + os << " \n"; + } - potatoes[which_node] = 0; - MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, 0, - MPI_COMM_WORLD); - } + //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); - myPotato = potatoes[which_node]; + // 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); - //recieve the number of integrableObject in current molecule - MPI_Recv(&nCurObj, 1, MPI_INT, which_node, myPotato, - MPI_COMM_WORLD, &istatus); - myPotato++; + // 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); - for(int l = 0; l < nCurObj; l++) { - if (potatoes[which_node] + 2 >= MAXTAG) { - // The potato was going to exceed the maximum value, - // so wrap this processor potato back to 0: + // send our buffer: + MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength, + MPI::CHAR, masterNode, 0); - potatoes[which_node] = 0; - MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, - 0, MPI_COMM_WORLD); - } + } + } + } + + if (worldRank == masterNode) { + os << " \n"; + } - MPI_Recv(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, - which_node, myPotato, MPI_COMM_WORLD, - &istatus); + 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); - atomTypeString = MPIatomTypeString; + if (sd->isRigidBody()) { + + RigidBody* rb = static_cast(sd); + int siteIndex = 0; + for (Atom* atom = rb->beginAtom(ai); atom != NULL; + atom = rb->nextAtom(ai)) { + buffer += prepareSiteLine(atom, ioIndex, siteIndex); + siteIndex++; + } + } + } + } - myPotato++; - - MPI_Recv(atomData, 13, MPI_DOUBLE, which_node, myPotato, - MPI_COMM_WORLD, &istatus); - myPotato++; - - MPI_Get_count(&istatus, MPI_DOUBLE, &msgLen); - - if (msgLen == 13) - isDirectional = 1; - else - isDirectional = 0; - - // If we've survived to here, format the line: - - if (!isDirectional) { - sprintf(writeLine, "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t", - atomTypeString, atomData[0], - atomData[1], atomData[2], - atomData[3], atomData[4], - atomData[5]); + 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 + + } - strcat(writeLine, - "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n"); - } else { - sprintf(writeLine, - "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", - atomTypeString, - atomData[0], - atomData[1], - atomData[2], - atomData[3], - atomData[4], - atomData[5], - atomData[6], - atomData[7], - atomData[8], - atomData[9], - atomData[10], - atomData[11], - atomData[12]); - } + std::string DumpWriter::prepareDumpLine(StuntDouble* sd) { + + int index = sd->getGlobalIntegrableObjectIndex(); + std::string type("pv"); + std::string line; + char tempBuffer[4096]; - os << writeLine; + Vector3d pos; + Vector3d vel; + pos = sd->getPos(); - } // end for(int l =0) + 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(); + } - potatoes[which_node] = myPotato; - } else { //master node has current molecule + vel = sd->getVel(); - mol = info_->getMoleculeByGlobalIndex(i); + 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(); + } - if (mol == NULL) { - sprintf(painCave.errMsg, "Molecule not found on node %d!", worldRank); - painCave.isFatal = 1; - simError(); - } - - for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { - - atomTypeString = integrableObject->getType().c_str(); + 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; - pos = integrableObject->getPos(); - vel = integrableObject->getVel(); + if (sd->isDirectional()) { + type += "qj"; + Quat4d q; + Vector3d ji; + q = sd->getQ(); - atomData[0] = pos[0]; - atomData[1] = pos[1]; - atomData[2] = pos[2]; + 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(); + } - atomData[3] = vel[0]; - atomData[4] = vel[1]; - atomData[5] = vel[2]; + ji = sd->getJ(); - isDirectional = 0; + 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(); + } - if (integrableObject->isDirectional()) { - isDirectional = 1; + sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e %13e", + q[0], q[1], q[2], q[3], + ji[0], ji[1], ji[2]); + line += tempBuffer; + } - q = integrableObject->getQ(); - ji = integrableObject->getJ(); + if (needForceVector_) { + 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; + } + } - for(int j = 0; j < 6; j++) { - atomData[j] = atomData[j]; - } + sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str()); + return std::string(tempBuffer); + } - atomData[6] = q[0]; - atomData[7] = q[1]; - atomData[8] = q[2]; - atomData[9] = q[3]; + std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) { + - atomData[10] = ji[0]; - atomData[11] = ji[1]; - atomData[12] = ji[2]; - } + std::string id; + std::string type; + std::string line; + char tempBuffer[4096]; - // If we've survived to here, format the line: - - if (!isDirectional) { - sprintf(writeLine, "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t", - atomTypeString, atomData[0], - atomData[1], atomData[2], - atomData[3], atomData[4], - atomData[5]); - - strcat(writeLine, - "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n"); - } else { - sprintf(writeLine, - "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", - atomTypeString, - atomData[0], - atomData[1], - atomData[2], - atomData[3], - atomData[4], - atomData[5], - atomData[6], - atomData[7], - atomData[8], - atomData[9], - atomData[10], - atomData[11], - atomData[12]); - } - - - os << writeLine; - - } //end for(iter = integrableObject.begin()) - } - } //end for(i = 0; i < mpiSim->getNmol()) - - os.flush(); - - sprintf(checkPointMsg, "Sucessfully took a dump.\n"); - MPIcheckPoint(); + 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; - delete [] potatoes; - } else { + 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; + } + } - // worldRank != 0, so I'm a remote node. + 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; + } - // Set my magic potato to 0: - myPotato = 0; + 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; + } + - for(int i = 0; i < info_->getNGlobalMolecules(); i++) { + sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str()); + return std::string(tempBuffer); + } - // Am I the node which has this integrableObject? - int whichNode = info_->getMolToProc(i); - if (whichNode == worldRank) { - if (myPotato + 1 >= MAXTAG) { + void DumpWriter::writeDump() { + writeFrame(*dumpFile_); + } - // The potato was going to exceed the maximum value, - // so wrap this processor potato back to 0 (and block until - // node 0 says we can go: + void DumpWriter::writeEor() { + std::ostream* eorStream; + +#ifdef IS_MPI + if (worldRank == 0) { +#endif // is_mpi + + eorStream = createOStream(eorFilename_); - MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, - &istatus); - } +#ifdef IS_MPI + } +#endif + + writeFrame(*eorStream); + +#ifdef IS_MPI + if (worldRank == 0) { +#endif + + writeClosing(*eorStream); + delete eorStream; + +#ifdef IS_MPI + } +#endif // is_mpi - mol = info_->getMoleculeByGlobalIndex(i); + } - - nCurObj = mol->getNIntegrableObjects(); - MPI_Send(&nCurObj, 1, MPI_INT, 0, myPotato, MPI_COMM_WORLD); - myPotato++; + void DumpWriter::writeDumpAndEor() { + std::vector buffers; + std::ostream* eorStream; +#ifdef IS_MPI + if (worldRank == 0) { +#endif // is_mpi - for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { + buffers.push_back(dumpFile_->rdbuf()); - if (myPotato + 2 >= MAXTAG) { + eorStream = createOStream(eorFilename_); - // The potato was going to exceed the maximum value, - // so wrap this processor potato back to 0 (and block until - // node 0 says we can go: + buffers.push_back(eorStream->rdbuf()); + +#ifdef IS_MPI + } +#endif // is_mpi - MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, - &istatus); - } + TeeBuf tbuf(buffers.begin(), buffers.end()); + std::ostream os(&tbuf); - atomTypeString = integrableObject->getType().c_str(); + writeFrame(os); - pos = integrableObject->getPos(); - vel = integrableObject->getVel(); +#ifdef IS_MPI + if (worldRank == 0) { +#endif // is_mpi + writeClosing(*eorStream); + delete eorStream; +#ifdef IS_MPI + } +#endif // is_mpi + + } - atomData[0] = pos[0]; - atomData[1] = pos[1]; - atomData[2] = pos[2]; + std::ostream* DumpWriter::createOStream(const std::string& filename) { - atomData[3] = vel[0]; - atomData[4] = vel[1]; - atomData[5] = vel[2]; - - isDirectional = 0; - - if (integrableObject->isDirectional()) { - isDirectional = 1; - - q = integrableObject->getQ(); - ji = integrableObject->getJ(); - - atomData[6] = q[0]; - atomData[7] = q[1]; - atomData[8] = q[2]; - atomData[9] = q[3]; - - atomData[10] = ji[0]; - atomData[11] = ji[1]; - atomData[12] = ji[2]; - } - - strncpy(MPIatomTypeString, atomTypeString, MINIBUFFERSIZE); - - // null terminate the std::string before sending (just in case): - MPIatomTypeString[MINIBUFFERSIZE - 1] = '\0'; - - MPI_Send(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0, - myPotato, MPI_COMM_WORLD); - - myPotato++; - - if (isDirectional) { - MPI_Send(atomData, 13, MPI_DOUBLE, 0, myPotato, - MPI_COMM_WORLD); - } else { - MPI_Send(atomData, 6, MPI_DOUBLE, 0, myPotato, - MPI_COMM_WORLD); - } - - myPotato++; - } - - } - - } - sprintf(checkPointMsg, "Sucessfully took a dump.\n"); - MPIcheckPoint(); + std::ostream* newOStream; +#ifdef HAVE_ZLIB + if (needCompression_) { + newOStream = new ogzstream(filename.c_str()); + } else { + newOStream = new std::ofstream(filename.c_str()); } +#else + newOStream = new std::ofstream(filename.c_str()); +#endif + //write out MetaData first + (*newOStream) << "" << std::endl; + (*newOStream) << " " << std::endl; + (*newOStream) << info_->getRawMetaData(); + (*newOStream) << " " << std::endl; + return newOStream; + } -#endif // is_mpi + void DumpWriter::writeClosing(std::ostream& os) { -} + os << "\n"; + os.flush(); + } -}//end namespace oopse +}//end namespace OpenMD