src/io/DumpWriter.cpp

/*
 * 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. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 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.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * 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"
#include "io/gzstream.hpp"
#include "io/Globals.hpp"


#ifdef IS_MPI
#include <mpi.h>
#endif //is_mpi

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();
    needParticlePot_   = simParams->getOutputParticlePotential();
    needFlucQ_         = simParams->getOutputFluctuatingCharges();
    needElectricField_ = simParams->getOutputElectricField();

    createDumpFile_ = true;
#ifdef HAVE_LIBZ
    if (needCompression_) {
      filename_ += ".gz";
      eorFilename_ += ".gz";
    }
#endif
    
#ifdef IS_MPI

    if (worldRank == 0) {
#endif // is_mpi
        
      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) 
    : 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();

    createDumpFile_ = true;
#ifdef HAVE_LIBZ
    if (needCompression_) {
      filename_ += ".gz";
      eorFilename_ += ".gz";
    }
#endif
    
#ifdef IS_MPI

    if (worldRank == 0) {
#endif // is_mpi

      
      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();

#ifdef HAVE_LIBZ
    if (needCompression_) {
      filename_ += ".gz";
      eorFilename_ += ".gz";
    }
#endif
    
#ifdef IS_MPI
    
    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
      
    }

    
#endif // is_mpi
    
  }

  DumpWriter::~DumpWriter() {

#ifdef IS_MPI

    if (worldRank == 0) {
#endif // is_mpi
      if (createDumpFile_){
        writeClosing(*dumpFile_);
        delete dumpFile_;
      }
#ifdef IS_MPI

    }

#endif // is_mpi

  }

  void DumpWriter::writeFrameProperties(std::ostream& os, Snapshot* s) {

    char buffer[1024];

    os << "    <FrameData>\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();
    if (isinf(chi) || isnan(chi) || 
        isinf(integralOfChiDt) || isnan(integralOfChiDt)) {      
      sprintf( painCave.errMsg,
               "DumpWriter detected a numerical error writing the thermostat");
      painCave.isFatal = 1;
      simError();
    }
    sprintf(buffer, "  Thermostat: %.10g , %.10g\n", chi, integralOfChiDt);
    os << buffer;

    Mat3x3d eta;
    eta = s->getEta();

    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),
            eta(0, 2), eta(1, 2), eta(2, 2));
    os << buffer;

    os << "    </FrameData>\n";
  }

  void DumpWriter::writeFrame(std::ostream& os) {

#ifdef IS_MPI
    MPI_Status istatus;
#endif

    Molecule* mol;
    StuntDouble* integrableObject;
    SimInfo::MoleculeIterator mi;
    Molecule::IntegrableObjectIterator ii;

#ifndef IS_MPI
    os << "  <Snapshot>\n";
 
    writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());

    os << "    <StuntDoubles>\n";
    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);
          
      }
    }    
    os << "    </StuntDoubles>\n";
    
    os << "  </Snapshot>\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 nProc;
    MPI_Comm_size(MPI_COMM_WORLD, &nProc);
    if (worldRank == masterNode) {      
      os << "  <Snapshot>\n";   
      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
      os << "    <StuntDoubles>\n";
        
      os << buffer;

      for (int i = 1; i < nProc; ++i) {

        // receive the length of the string buffer that was
        // prepared by processor i

        MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
        int recvLength;
        MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
        char* recvBuffer = new char[recvLength];
        if (recvBuffer == NULL) {
        } else {
          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
          os << recvBuffer;
          delete [] recvBuffer;
        }
      } 
      os << "    </StuntDoubles>\n";
      
      os << "  </Snapshot>\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);
        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);
        }
      }
    }

#endif // is_mpi

  }

  std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
        
    int index = integrableObject->getGlobalIntegrableObjectIndex();
    std::string type("pv");
    std::string line;
    char tempBuffer[4096];

    Vector3d pos;
    Vector3d vel;
    pos = integrableObject->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 = integrableObject->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()) {
      type += "qj";
      Quat4d q;
      Vector3d ji;
      q = integrableObject->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 = integrableObject->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]);
      line += tempBuffer;
    }

    if (needForceVector_) {
      type += "f";
      Vector3d frc = integrableObject->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 (integrableObject->isDirectional()) {
        type += "t";
        Vector3d trq = integrableObject->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;
      }      
    }

    if (needParticlePot_) {
      type += "u";
      RealType particlePot = integrableObject->getParticlePot();
      if (isinf(particlePot) || isnan(particlePot)) {      
        sprintf( painCave.errMsg,
                 "DumpWriter detected a numerical error writing the particle "
                 " potential for object %d", index);      
        painCave.isFatal = 1;
        simError();
      }
      sprintf(tempBuffer, " %13e", particlePot);
      line += tempBuffer;
    }
    
    if (needFlucQ_) {
      type += "cw";
      RealType fqPos = integrableObject->getFlucQPos();
      if (isinf(fqPos) || isnan(fqPos) ) {      
        sprintf( painCave.errMsg,
                 "DumpWriter detected a numerical error writing the"
                 " fluctuating charge for object %d", index);      
        painCave.isFatal = 1;
        simError();
      }
      sprintf(tempBuffer, " %13e ", fqPos);
      line += tempBuffer;
    
      RealType fqVel = integrableObject->getFlucQVel();
      if (isinf(fqVel) || isnan(fqVel) ) {      
        sprintf( painCave.errMsg,
                 "DumpWriter detected a numerical error writing the"
                 " fluctuating charge velocity for object %d", index);      
        painCave.isFatal = 1;
        simError();
      }
      sprintf(tempBuffer, " %13e ", fqVel);
      line += tempBuffer;

      if (needForceVector_) {
        type += "g";
        RealType fqFrc = integrableObject->getFlucQFrc();        
        if (isinf(fqFrc) || isnan(fqFrc) ) {      
          sprintf( painCave.errMsg,
                   "DumpWriter detected a numerical error writing the"
                   " fluctuating charge force for object %d", index);      
          painCave.isFatal = 1;
          simError();
        }
        sprintf(tempBuffer, " %13e ", fqFrc);        
        line += tempBuffer;
      }
    }

    if (needElectricField_) {
      type += "e";
      Vector3d eField= integrableObject->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 %d", index);      
        painCave.isFatal = 1;
        simError();
      }
      sprintf(tempBuffer, " %13e %13e %13e",
              eField[0], eField[1], eField[2]);
      line += tempBuffer;
    }

    sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
    return std::string(tempBuffer);
  }

  void DumpWriter::writeDump() {
    writeFrame(*dumpFile_);
  }

  void DumpWriter::writeEor() {
    std::ostream* eorStream;
    
#ifdef IS_MPI
    if (worldRank == 0) {
#endif // is_mpi

      eorStream = createOStream(eorFilename_);

#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  

  }


  void DumpWriter::writeDumpAndEor() {
    std::vector<std::streambuf*> buffers;
    std::ostream* eorStream;
#ifdef IS_MPI
    if (worldRank == 0) {
#endif // is_mpi

      buffers.push_back(dumpFile_->rdbuf());

      eorStream = createOStream(eorFilename_);

      buffers.push_back(eorStream->rdbuf());
        
#ifdef IS_MPI
    }
#endif // is_mpi    

    TeeBuf tbuf(buffers.begin(), buffers.end());
    std::ostream os(&tbuf);

    writeFrame(os);

#ifdef IS_MPI
    if (worldRank == 0) {
#endif // is_mpi
      writeClosing(*eorStream);
      delete eorStream;
#ifdef IS_MPI
    }
#endif // is_mpi  
    
  }

  std::ostream* DumpWriter::createOStream(const std::string& filename) {

    std::ostream* newOStream;
#ifdef HAVE_LIBZ 
    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) << "<OpenMD version=2>" << std::endl;
    (*newOStream) << "  <MetaData>" << std::endl;
    (*newOStream) << info_->getRawMetaData();
    (*newOStream) << "  </MetaData>" << std::endl;
    return newOStream;
  }

  void DumpWriter::writeClosing(std::ostream& os) {

    os << "</OpenMD>\n";
    os.flush();
  }

}//end namespace OpenMD
Revision:	1746
Committed:	Wed Jun 6 02:18:54 2012 UTC (12 years, 10 months ago) by gezelter
File size:	18773 byte(s)
Log Message:	added a minimizer parsing block
#	User	Rev	Content
1	gezelter	507	/*
2	gezelter	1390	* Copyright (c) 2009 The University of Notre Dame. All Rights Reserved.
3	gezelter	246	*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	gezelter	246	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	gezelter	246	* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
28			* arising out of the use of or inability to use software, even if the
29			* University of Notre Dame has been advised of the possibility of
30			* such damages.
31	gezelter	1390	*
32			* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			* research, please cite the appropriate papers when you publish your
34			* work. Good starting points are:
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).
39	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	gezelter	246	*/
42
43			#include "io/DumpWriter.hpp"
44			#include "primitives/Molecule.hpp"
45			#include "utils/simError.h"
46	tim	376	#include "io/basic_teebuf.hpp"
47	tim	615	#include "io/gzstream.hpp"
48			#include "io/Globals.hpp"
49
50	gezelter	1437
51	gezelter	2	#ifdef IS_MPI
52			#include <mpi.h>
53	gezelter	246	#endif //is_mpi
54	gezelter	2
55	gezelter	1437	using namespace std;
56	gezelter	1390	namespace OpenMD {
57	gezelter	2
58	gezelter	507	DumpWriter::DumpWriter(SimInfo* info)
59			: info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
60	tim	615
61			Globals* simParams = info->getSimParams();
62	gezelter	1714	needCompression_ = simParams->getCompressDumpFile();
63			needForceVector_ = simParams->getOutputForceVector();
64			needParticlePot_ = simParams->getOutputParticlePotential();
65			needFlucQ_ = simParams->getOutputFluctuatingCharges();
66			needElectricField_ = simParams->getOutputElectricField();
67
68	chuckv	791	createDumpFile_ = true;
69	tim	619	#ifdef HAVE_LIBZ
70	tim	615	if (needCompression_) {
71	tim	1024	filename_ += ".gz";
72			eorFilename_ += ".gz";
73	tim	615	}
74	tim	619	#endif
75	tim	615
76	tim	376	#ifdef IS_MPI
77
78	tim	1024	if (worldRank == 0) {
79	tim	376	#endif // is_mpi
80	chuckv	791
81	tim	1024	dumpFile_ = createOStream(filename_);
82	tim	615
83	tim	1024	if (!dumpFile_) {
84			sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
85			filename_.c_str());
86			painCave.isFatal = 1;
87			simError();
88			}
89	tim	376
90			#ifdef IS_MPI
91
92	tim	1024	}
93	tim	376
94			#endif // is_mpi
95
96	tim	1024	}
97	tim	376
98
99	gezelter	507	DumpWriter::DumpWriter(SimInfo* info, const std::string& filename)
100			: info_(info), filename_(filename){
101	tim	615
102			Globals* simParams = info->getSimParams();
103			eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
104
105	gezelter	1714	needCompression_ = simParams->getCompressDumpFile();
106			needForceVector_ = simParams->getOutputForceVector();
107			needParticlePot_ = simParams->getOutputParticlePotential();
108			needFlucQ_ = simParams->getOutputFluctuatingCharges();
109			needElectricField_ = simParams->getOutputElectricField();
110
111	chuckv	791	createDumpFile_ = true;
112	tim	619	#ifdef HAVE_LIBZ
113	tim	615	if (needCompression_) {
114	tim	1024	filename_ += ".gz";
115			eorFilename_ += ".gz";
116	tim	615	}
117	tim	619	#endif
118	tim	615
119	gezelter	246	#ifdef IS_MPI
120	gezelter	2
121	tim	1024	if (worldRank == 0) {
122	gezelter	246	#endif // is_mpi
123	gezelter	2
124	chuckv	791
125	tim	1024	dumpFile_ = createOStream(filename_);
126	tim	615
127	tim	1024	if (!dumpFile_) {
128			sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
129			filename_.c_str());
130			painCave.isFatal = 1;
131			simError();
132			}
133	gezelter	2
134			#ifdef IS_MPI
135
136	tim	1024	}
137	gezelter	2
138			#endif // is_mpi
139	gezelter	246
140	tim	1024	}
141	chuckv	791
142			DumpWriter::DumpWriter(SimInfo* info, const std::string& filename, bool writeDumpFile)
143	tim	1024	: info_(info), filename_(filename){
144	chuckv	791
145			Globals* simParams = info->getSimParams();
146			eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
147
148	gezelter	1714	needCompression_ = simParams->getCompressDumpFile();
149			needForceVector_ = simParams->getOutputForceVector();
150			needParticlePot_ = simParams->getOutputParticlePotential();
151			needFlucQ_ = simParams->getOutputFluctuatingCharges();
152			needElectricField_ = simParams->getOutputElectricField();
153
154	chuckv	791	#ifdef HAVE_LIBZ
155			if (needCompression_) {
156			filename_ += ".gz";
157			eorFilename_ += ".gz";
158			}
159			#endif
160
161			#ifdef IS_MPI
162
163			if (worldRank == 0) {
164			#endif // is_mpi
165
166			createDumpFile_ = writeDumpFile;
167			if (createDumpFile_) {
168			dumpFile_ = createOStream(filename_);
169
170			if (!dumpFile_) {
171			sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
172			filename_.c_str());
173			painCave.isFatal = 1;
174			simError();
175			}
176			}
177			#ifdef IS_MPI
178
179			}
180	tim	1024
181	chuckv	791
182			#endif // is_mpi
183
184			}
185	gezelter	2
186	gezelter	507	DumpWriter::~DumpWriter() {
187	gezelter	2
188			#ifdef IS_MPI
189	gezelter	246
190			if (worldRank == 0) {
191	gezelter	2	#endif // is_mpi
192	chuckv	791	if (createDumpFile_){
193	tim	1024	writeClosing(*dumpFile_);
194	chuckv	791	delete dumpFile_;
195			}
196	gezelter	2	#ifdef IS_MPI
197	gezelter	246
198			}
199
200	gezelter	2	#endif // is_mpi
201	gezelter	246
202	gezelter	507	}
203	gezelter	2
204	tim	1024	void DumpWriter::writeFrameProperties(std::ostream& os, Snapshot* s) {
205	gezelter	2
206	tim	1024	char buffer[1024];
207
208			os << " <FrameData>\n";
209
210			RealType currentTime = s->getTime();
211	gezelter	1437
212			if (isinf(currentTime) \|\| isnan(currentTime)) {
213			sprintf( painCave.errMsg,
214			"DumpWriter detected a numerical error writing the time");
215			painCave.isFatal = 1;
216			simError();
217			}
218
219	gezelter	1025	sprintf(buffer, " Time: %.10g\n", currentTime);
220	tim	1024	os << buffer;
221
222	gezelter	246	Mat3x3d hmat;
223			hmat = s->getHmat();
224	gezelter	1437
225			for (unsigned int i = 0; i < 3; i++) {
226			for (unsigned int j = 0; j < 3; j++) {
227			if (isinf(hmat(i,j)) \|\| isnan(hmat(i,j))) {
228			sprintf( painCave.errMsg,
229			"DumpWriter detected a numerical error writing the box");
230			painCave.isFatal = 1;
231			simError();
232			}
233			}
234			}
235
236	tim	1024	sprintf(buffer, " Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
237			hmat(0, 0), hmat(1, 0), hmat(2, 0),
238			hmat(0, 1), hmat(1, 1), hmat(2, 1),
239			hmat(0, 2), hmat(1, 2), hmat(2, 2));
240			os << buffer;
241	gezelter	2
242	tim	1024	RealType chi = s->getChi();
243			RealType integralOfChiDt = s->getIntegralOfChiDt();
244	gezelter	1437	if (isinf(chi) \|\| isnan(chi) \|\|
245			isinf(integralOfChiDt) \|\| isnan(integralOfChiDt)) {
246			sprintf( painCave.errMsg,
247			"DumpWriter detected a numerical error writing the thermostat");
248			painCave.isFatal = 1;
249			simError();
250			}
251	tim	1024	sprintf(buffer, " Thermostat: %.10g , %.10g\n", chi, integralOfChiDt);
252			os << buffer;
253	gezelter	246
254	tim	1024	Mat3x3d eta;
255			eta = s->getEta();
256	gezelter	1437
257			for (unsigned int i = 0; i < 3; i++) {
258			for (unsigned int j = 0; j < 3; j++) {
259			if (isinf(eta(i,j)) \|\| isnan(eta(i,j))) {
260			sprintf( painCave.errMsg,
261			"DumpWriter detected a numerical error writing the barostat");
262			painCave.isFatal = 1;
263			simError();
264			}
265			}
266			}
267
268	tim	1024	sprintf(buffer, " Barostat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
269			eta(0, 0), eta(1, 0), eta(2, 0),
270			eta(0, 1), eta(1, 1), eta(2, 1),
271			eta(0, 2), eta(1, 2), eta(2, 2));
272			os << buffer;
273	gezelter	246
274	tim	1024	os << " </FrameData>\n";
275	gezelter	507	}
276	gezelter	2
277	gezelter	507	void DumpWriter::writeFrame(std::ostream& os) {
278	gezelter	246
279	tim	1024	#ifdef IS_MPI
280			MPI_Status istatus;
281			#endif
282	gezelter	246
283			Molecule* mol;
284			StuntDouble* integrableObject;
285			SimInfo::MoleculeIterator mi;
286			Molecule::IntegrableObjectIterator ii;
287	gezelter	2
288	gezelter	246	#ifndef IS_MPI
289	tim	1024	os << " <Snapshot>\n";
290	gezelter	1025
291	tim	1024	writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
292	gezelter	2
293	tim	1024	os << " <StuntDoubles>\n";
294	gezelter	246	for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
295	gezelter	2
296	xsun	1217
297			for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
298	tim	1024	integrableObject = mol->nextIntegrableObject(ii)) {
299	xsun	1217	os << prepareDumpLine(integrableObject);
300
301	tim	1024	}
302			}
303			os << " </StuntDoubles>\n";
304	xsun	1217
305	tim	1024	os << " </Snapshot>\n";
306	gezelter	2
307	tim	1024	os.flush();
308			#else
309			//every node prepares the dump lines for integrable objects belong to itself
310			std::string buffer;
311			for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
312	xsun	1217
313
314	tim	1024	for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
315			integrableObject = mol->nextIntegrableObject(ii)) {
316	xsun	1217	buffer += prepareDumpLine(integrableObject);
317	gezelter	507	}
318	gezelter	2	}
319	tim	1024
320	gezelter	246	const int masterNode = 0;
321	gezelter	1629	int nProc;
322			MPI_Comm_size(MPI_COMM_WORLD, &nProc);
323	tim	1024	if (worldRank == masterNode) {
324			os << " <Snapshot>\n";
325			writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
326			os << " <StuntDoubles>\n";
327
328	gezelter	1025	os << buffer;
329
330	tim	1024	for (int i = 1; i < nProc; ++i) {
331	gezelter	2
332	tim	1024	// receive the length of the string buffer that was
333			// prepared by processor i
334	gezelter	2
335	gezelter	1629	MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
336	tim	1024	int recvLength;
337			MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
338			char* recvBuffer = new char[recvLength];
339			if (recvBuffer == NULL) {
340			} else {
341	gezelter	1025	MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
342	tim	1024	os << recvBuffer;
343	gezelter	1313	delete [] recvBuffer;
344	tim	1024	}
345			}
346			os << " </StuntDoubles>\n";
347
348			os << " </Snapshot>\n";
349	gezelter	507	os.flush();
350	gezelter	246	} else {
351	gezelter	1025	int sendBufferLength = buffer.size() + 1;
352	gezelter	1629	int myturn = 0;
353			for (int i = 1; i < nProc; ++i){
354			MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
355			if (myturn == worldRank){
356			MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
357			MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
358			}
359			}
360	tim	1024	}
361	gezelter	2
362	tim	1024	#endif // is_mpi
363	gezelter	2
364	tim	1024	}
365	gezelter	2
366	tim	1024	std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
367
368			int index = integrableObject->getGlobalIntegrableObjectIndex();
369			std::string type("pv");
370			std::string line;
371			char tempBuffer[4096];
372	gezelter	2
373	tim	1024	Vector3d pos;
374			Vector3d vel;
375			pos = integrableObject->getPos();
376	gezelter	1437
377			if (isinf(pos[0]) \|\| isnan(pos[0]) \|\|
378			isinf(pos[1]) \|\| isnan(pos[1]) \|\|
379			isinf(pos[2]) \|\| isnan(pos[2]) ) {
380			sprintf( painCave.errMsg,
381			"DumpWriter detected a numerical error writing the position"
382			" for object %d", index);
383			painCave.isFatal = 1;
384			simError();
385			}
386
387	tim	1024	vel = integrableObject->getVel();
388	gezelter	1437
389			if (isinf(vel[0]) \|\| isnan(vel[0]) \|\|
390			isinf(vel[1]) \|\| isnan(vel[1]) \|\|
391			isinf(vel[2]) \|\| isnan(vel[2]) ) {
392			sprintf( painCave.errMsg,
393			"DumpWriter detected a numerical error writing the velocity"
394			" for object %d", index);
395			painCave.isFatal = 1;
396			simError();
397			}
398
399	gezelter	1025	sprintf(tempBuffer, "%18.10g %18.10g %18.10g %13e %13e %13e",
400	tim	1024	pos[0], pos[1], pos[2],
401			vel[0], vel[1], vel[2]);
402			line += tempBuffer;
403	gezelter	2
404	tim	1024	if (integrableObject->isDirectional()) {
405			type += "qj";
406			Quat4d q;
407			Vector3d ji;
408			q = integrableObject->getQ();
409	gezelter	1437
410			if (isinf(q[0]) \|\| isnan(q[0]) \|\|
411			isinf(q[1]) \|\| isnan(q[1]) \|\|
412			isinf(q[2]) \|\| isnan(q[2]) \|\|
413			isinf(q[3]) \|\| isnan(q[3]) ) {
414			sprintf( painCave.errMsg,
415			"DumpWriter detected a numerical error writing the quaternion"
416			" for object %d", index);
417			painCave.isFatal = 1;
418			simError();
419			}
420
421	tim	1024	ji = integrableObject->getJ();
422	gezelter	1437
423			if (isinf(ji[0]) \|\| isnan(ji[0]) \|\|
424			isinf(ji[1]) \|\| isnan(ji[1]) \|\|
425			isinf(ji[2]) \|\| isnan(ji[2]) ) {
426			sprintf( painCave.errMsg,
427			"DumpWriter detected a numerical error writing the angular"
428			" momentum for object %d", index);
429			painCave.isFatal = 1;
430			simError();
431			}
432
433	gezelter	1025	sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e %13e",
434	tim	1024	q[0], q[1], q[2], q[3],
435			ji[0], ji[1], ji[2]);
436			line += tempBuffer;
437			}
438	gezelter	2
439	tim	1024	if (needForceVector_) {
440	gezelter	1437	type += "f";
441	gezelter	1714	Vector3d frc = integrableObject->getFrc();
442	gezelter	1437	if (isinf(frc[0]) \|\| isnan(frc[0]) \|\|
443			isinf(frc[1]) \|\| isnan(frc[1]) \|\|
444			isinf(frc[2]) \|\| isnan(frc[2]) ) {
445			sprintf( painCave.errMsg,
446			"DumpWriter detected a numerical error writing the force"
447			" for object %d", index);
448			painCave.isFatal = 1;
449			simError();
450			}
451			sprintf(tempBuffer, " %13e %13e %13e",
452			frc[0], frc[1], frc[2]);
453	tim	1024	line += tempBuffer;
454	gezelter	1437
455			if (integrableObject->isDirectional()) {
456			type += "t";
457	gezelter	1714	Vector3d trq = integrableObject->getTrq();
458	gezelter	1437	if (isinf(trq[0]) \|\| isnan(trq[0]) \|\|
459			isinf(trq[1]) \|\| isnan(trq[1]) \|\|
460			isinf(trq[2]) \|\| isnan(trq[2]) ) {
461			sprintf( painCave.errMsg,
462			"DumpWriter detected a numerical error writing the torque"
463			" for object %d", index);
464			painCave.isFatal = 1;
465			simError();
466	gezelter	1714	}
467	gezelter	1437	sprintf(tempBuffer, " %13e %13e %13e",
468			trq[0], trq[1], trq[2]);
469			line += tempBuffer;
470			}
471	gezelter	246	}
472	gezelter	1714
473	gezelter	1629	if (needParticlePot_) {
474			type += "u";
475	gezelter	1714	RealType particlePot = integrableObject->getParticlePot();
476	gezelter	1629	if (isinf(particlePot) \|\| isnan(particlePot)) {
477			sprintf( painCave.errMsg,
478			"DumpWriter detected a numerical error writing the particle "
479			" potential for object %d", index);
480			painCave.isFatal = 1;
481			simError();
482			}
483			sprintf(tempBuffer, " %13e", particlePot);
484			line += tempBuffer;
485			}
486	gezelter	1437
487	gezelter	1714	if (needFlucQ_) {
488			type += "cw";
489			RealType fqPos = integrableObject->getFlucQPos();
490			if (isinf(fqPos) \|\| isnan(fqPos) ) {
491			sprintf( painCave.errMsg,
492			"DumpWriter detected a numerical error writing the"
493			" fluctuating charge for object %d", index);
494			painCave.isFatal = 1;
495			simError();
496			}
497			sprintf(tempBuffer, " %13e ", fqPos);
498			line += tempBuffer;
499
500			RealType fqVel = integrableObject->getFlucQVel();
501			if (isinf(fqVel) \|\| isnan(fqVel) ) {
502			sprintf( painCave.errMsg,
503			"DumpWriter detected a numerical error writing the"
504			" fluctuating charge velocity for object %d", index);
505			painCave.isFatal = 1;
506			simError();
507			}
508			sprintf(tempBuffer, " %13e ", fqVel);
509			line += tempBuffer;
510
511			if (needForceVector_) {
512			type += "g";
513			RealType fqFrc = integrableObject->getFlucQFrc();
514			if (isinf(fqFrc) \|\| isnan(fqFrc) ) {
515			sprintf( painCave.errMsg,
516			"DumpWriter detected a numerical error writing the"
517			" fluctuating charge force for object %d", index);
518			painCave.isFatal = 1;
519			simError();
520			}
521			sprintf(tempBuffer, " %13e ", fqFrc);
522			line += tempBuffer;
523			}
524			}
525
526			if (needElectricField_) {
527			type += "e";
528			Vector3d eField= integrableObject->getElectricField();
529			if (isinf(eField[0]) \|\| isnan(eField[0]) \|\|
530			isinf(eField[1]) \|\| isnan(eField[1]) \|\|
531			isinf(eField[2]) \|\| isnan(eField[2]) ) {
532			sprintf( painCave.errMsg,
533			"DumpWriter detected a numerical error writing the electric"
534			" field for object %d", index);
535			painCave.isFatal = 1;
536			simError();
537			}
538			sprintf(tempBuffer, " %13e %13e %13e",
539			eField[0], eField[1], eField[2]);
540			line += tempBuffer;
541			}
542
543	gezelter	1025	sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
544	tim	1024	return std::string(tempBuffer);
545	gezelter	507	}
546	gezelter	2
547	gezelter	507	void DumpWriter::writeDump() {
548	tim	615	writeFrame(*dumpFile_);
549	gezelter	507	}
550	tim	376
551	gezelter	507	void DumpWriter::writeEor() {
552	tim	615	std::ostream* eorStream;
553	tim	376
554			#ifdef IS_MPI
555			if (worldRank == 0) {
556			#endif // is_mpi
557
558	tim	615	eorStream = createOStream(eorFilename_);
559	tim	376
560			#ifdef IS_MPI
561			}
562			#endif // is_mpi
563
564	tim	615	writeFrame(*eorStream);
565
566			#ifdef IS_MPI
567			if (worldRank == 0) {
568			#endif // is_mpi
569	tim	1024	writeClosing(*eorStream);
570			delete eorStream;
571	tim	615	#ifdef IS_MPI
572			}
573			#endif // is_mpi
574
575	gezelter	507	}
576	tim	376
577
578	gezelter	507	void DumpWriter::writeDumpAndEor() {
579	tim	376	std::vector<std::streambuf*> buffers;
580	tim	615	std::ostream* eorStream;
581	tim	376	#ifdef IS_MPI
582			if (worldRank == 0) {
583			#endif // is_mpi
584
585	tim	615	buffers.push_back(dumpFile_->rdbuf());
586	tim	376
587	tim	615	eorStream = createOStream(eorFilename_);
588	tim	376
589	tim	615	buffers.push_back(eorStream->rdbuf());
590	tim	376
591			#ifdef IS_MPI
592			}
593			#endif // is_mpi
594
595			TeeBuf tbuf(buffers.begin(), buffers.end());
596			std::ostream os(&tbuf);
597
598			writeFrame(os);
599	tim	615
600			#ifdef IS_MPI
601			if (worldRank == 0) {
602			#endif // is_mpi
603	tim	1024	writeClosing(*eorStream);
604			delete eorStream;
605	tim	615	#ifdef IS_MPI
606			}
607			#endif // is_mpi
608	tim	376
609	gezelter	507	}
610	tim	376
611	tim	1024	std::ostream* DumpWriter::createOStream(const std::string& filename) {
612	tim	619
613	tim	615	std::ostream* newOStream;
614	tim	619	#ifdef HAVE_LIBZ
615	tim	615	if (needCompression_) {
616	tim	1024	newOStream = new ogzstream(filename.c_str());
617	tim	615	} else {
618	tim	1024	newOStream = new std::ofstream(filename.c_str());
619	tim	615	}
620	tim	619	#else
621			newOStream = new std::ofstream(filename.c_str());
622			#endif
623	tim	1024	//write out MetaData first
624	gezelter	1746	(*newOStream) << "<OpenMD version=2>" << std::endl;
625	tim	1024	(*newOStream) << " <MetaData>" << std::endl;
626			(*newOStream) << info_->getRawMetaData();
627			(*newOStream) << " </MetaData>" << std::endl;
628	tim	615	return newOStream;
629	tim	1024	}
630	tim	376
631	tim	1024	void DumpWriter::writeClosing(std::ostream& os) {
632
633	gezelter	1390	os << "</OpenMD>\n";
634	tim	1024	os.flush();
635			}
636
637	gezelter	1390	}//end namespace OpenMD