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]  Vardeman & Gezelter, in progress (2009).                        
 */
 
#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();
    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();
    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();
    
#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;

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

      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;
    }
    
    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=1>" << 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:	1610
Committed:	Fri Aug 12 14:37:25 2011 UTC (13 years, 8 months ago) by gezelter
File size:	16377 byte(s)
Log Message:	Fixed a clang compilation problem, added the ability to output particle potential in the dump files.
#	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			* [4] Vardeman & Gezelter, in progress (2009).
40	gezelter	246	*/
41
42			#include "io/DumpWriter.hpp"
43			#include "primitives/Molecule.hpp"
44			#include "utils/simError.h"
45	tim	376	#include "io/basic_teebuf.hpp"
46	tim	615	#include "io/gzstream.hpp"
47			#include "io/Globals.hpp"
48
49	gezelter	1437
50	gezelter	2	#ifdef IS_MPI
51			#include <mpi.h>
52	gezelter	246	#endif //is_mpi
53	gezelter	2
54	gezelter	1437	using namespace std;
55	gezelter	1390	namespace OpenMD {
56	gezelter	2
57	gezelter	507	DumpWriter::DumpWriter(SimInfo* info)
58			: info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
59	tim	615
60			Globals* simParams = info->getSimParams();
61			needCompression_ = simParams->getCompressDumpFile();
62	chrisfen	726	needForceVector_ = simParams->getOutputForceVector();
63	gezelter	1610	needParticlePot_ = simParams->getOutputParticlePotential();
64	chuckv	791	createDumpFile_ = true;
65	tim	619	#ifdef HAVE_LIBZ
66	tim	615	if (needCompression_) {
67	tim	1024	filename_ += ".gz";
68			eorFilename_ += ".gz";
69	tim	615	}
70	tim	619	#endif
71	tim	615
72	tim	376	#ifdef IS_MPI
73
74	tim	1024	if (worldRank == 0) {
75	tim	376	#endif // is_mpi
76	chuckv	791
77	tim	1024	dumpFile_ = createOStream(filename_);
78	tim	615
79	tim	1024	if (!dumpFile_) {
80			sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
81			filename_.c_str());
82			painCave.isFatal = 1;
83			simError();
84			}
85	tim	376
86			#ifdef IS_MPI
87
88	tim	1024	}
89	tim	376
90			#endif // is_mpi
91
92	tim	1024	}
93	tim	376
94
95	gezelter	507	DumpWriter::DumpWriter(SimInfo* info, const std::string& filename)
96			: info_(info), filename_(filename){
97	tim	615
98			Globals* simParams = info->getSimParams();
99			eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
100
101			needCompression_ = simParams->getCompressDumpFile();
102	chrisfen	726	needForceVector_ = simParams->getOutputForceVector();
103	gezelter	1610	needParticlePot_ = simParams->getOutputParticlePotential();
104	chuckv	791	createDumpFile_ = true;
105	tim	619	#ifdef HAVE_LIBZ
106	tim	615	if (needCompression_) {
107	tim	1024	filename_ += ".gz";
108			eorFilename_ += ".gz";
109	tim	615	}
110	tim	619	#endif
111	tim	615
112	gezelter	246	#ifdef IS_MPI
113	gezelter	2
114	tim	1024	if (worldRank == 0) {
115	gezelter	246	#endif // is_mpi
116	gezelter	2
117	chuckv	791
118	tim	1024	dumpFile_ = createOStream(filename_);
119	tim	615
120	tim	1024	if (!dumpFile_) {
121			sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
122			filename_.c_str());
123			painCave.isFatal = 1;
124			simError();
125			}
126	gezelter	2
127			#ifdef IS_MPI
128
129	tim	1024	}
130	gezelter	2
131			#endif // is_mpi
132	gezelter	246
133	tim	1024	}
134	chuckv	791
135			DumpWriter::DumpWriter(SimInfo* info, const std::string& filename, bool writeDumpFile)
136	tim	1024	: info_(info), filename_(filename){
137	chuckv	791
138			Globals* simParams = info->getSimParams();
139			eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
140
141			needCompression_ = simParams->getCompressDumpFile();
142			needForceVector_ = simParams->getOutputForceVector();
143	gezelter	1610	needParticlePot_ = simParams->getOutputParticlePotential();
144	chuckv	791
145			#ifdef HAVE_LIBZ
146			if (needCompression_) {
147			filename_ += ".gz";
148			eorFilename_ += ".gz";
149			}
150			#endif
151
152			#ifdef IS_MPI
153
154			if (worldRank == 0) {
155			#endif // is_mpi
156
157			createDumpFile_ = writeDumpFile;
158			if (createDumpFile_) {
159			dumpFile_ = createOStream(filename_);
160
161			if (!dumpFile_) {
162			sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
163			filename_.c_str());
164			painCave.isFatal = 1;
165			simError();
166			}
167			}
168			#ifdef IS_MPI
169
170			}
171	tim	1024
172	chuckv	791
173			#endif // is_mpi
174
175			}
176	gezelter	2
177	gezelter	507	DumpWriter::~DumpWriter() {
178	gezelter	2
179			#ifdef IS_MPI
180	gezelter	246
181			if (worldRank == 0) {
182	gezelter	2	#endif // is_mpi
183	chuckv	791	if (createDumpFile_){
184	tim	1024	writeClosing(*dumpFile_);
185	chuckv	791	delete dumpFile_;
186			}
187	gezelter	2	#ifdef IS_MPI
188	gezelter	246
189			}
190
191	gezelter	2	#endif // is_mpi
192	gezelter	246
193	gezelter	507	}
194	gezelter	2
195	tim	1024	void DumpWriter::writeFrameProperties(std::ostream& os, Snapshot* s) {
196	gezelter	2
197	tim	1024	char buffer[1024];
198
199			os << " <FrameData>\n";
200
201			RealType currentTime = s->getTime();
202	gezelter	1437
203			if (isinf(currentTime) \|\| isnan(currentTime)) {
204			sprintf( painCave.errMsg,
205			"DumpWriter detected a numerical error writing the time");
206			painCave.isFatal = 1;
207			simError();
208			}
209
210	gezelter	1025	sprintf(buffer, " Time: %.10g\n", currentTime);
211	tim	1024	os << buffer;
212
213	gezelter	246	Mat3x3d hmat;
214			hmat = s->getHmat();
215	gezelter	1437
216			for (unsigned int i = 0; i < 3; i++) {
217			for (unsigned int j = 0; j < 3; j++) {
218			if (isinf(hmat(i,j)) \|\| isnan(hmat(i,j))) {
219			sprintf( painCave.errMsg,
220			"DumpWriter detected a numerical error writing the box");
221			painCave.isFatal = 1;
222			simError();
223			}
224			}
225			}
226
227	tim	1024	sprintf(buffer, " Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
228			hmat(0, 0), hmat(1, 0), hmat(2, 0),
229			hmat(0, 1), hmat(1, 1), hmat(2, 1),
230			hmat(0, 2), hmat(1, 2), hmat(2, 2));
231			os << buffer;
232	gezelter	2
233	tim	1024	RealType chi = s->getChi();
234			RealType integralOfChiDt = s->getIntegralOfChiDt();
235	gezelter	1437	if (isinf(chi) \|\| isnan(chi) \|\|
236			isinf(integralOfChiDt) \|\| isnan(integralOfChiDt)) {
237			sprintf( painCave.errMsg,
238			"DumpWriter detected a numerical error writing the thermostat");
239			painCave.isFatal = 1;
240			simError();
241			}
242	tim	1024	sprintf(buffer, " Thermostat: %.10g , %.10g\n", chi, integralOfChiDt);
243			os << buffer;
244	gezelter	246
245	tim	1024	Mat3x3d eta;
246			eta = s->getEta();
247	gezelter	1437
248			for (unsigned int i = 0; i < 3; i++) {
249			for (unsigned int j = 0; j < 3; j++) {
250			if (isinf(eta(i,j)) \|\| isnan(eta(i,j))) {
251			sprintf( painCave.errMsg,
252			"DumpWriter detected a numerical error writing the barostat");
253			painCave.isFatal = 1;
254			simError();
255			}
256			}
257			}
258
259	tim	1024	sprintf(buffer, " Barostat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
260			eta(0, 0), eta(1, 0), eta(2, 0),
261			eta(0, 1), eta(1, 1), eta(2, 1),
262			eta(0, 2), eta(1, 2), eta(2, 2));
263			os << buffer;
264	gezelter	246
265	tim	1024	os << " </FrameData>\n";
266	gezelter	507	}
267	gezelter	2
268	gezelter	507	void DumpWriter::writeFrame(std::ostream& os) {
269	gezelter	246
270	tim	1024	#ifdef IS_MPI
271			MPI_Status istatus;
272			#endif
273	gezelter	246
274			Molecule* mol;
275			StuntDouble* integrableObject;
276			SimInfo::MoleculeIterator mi;
277			Molecule::IntegrableObjectIterator ii;
278	gezelter	2
279	gezelter	246	#ifndef IS_MPI
280	tim	1024	os << " <Snapshot>\n";
281	gezelter	1025
282	tim	1024	writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
283	gezelter	2
284	tim	1024	os << " <StuntDoubles>\n";
285	gezelter	246	for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
286	gezelter	2
287	xsun	1217
288			for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
289	tim	1024	integrableObject = mol->nextIntegrableObject(ii)) {
290	xsun	1217	os << prepareDumpLine(integrableObject);
291
292	tim	1024	}
293			}
294			os << " </StuntDoubles>\n";
295	xsun	1217
296	tim	1024	os << " </Snapshot>\n";
297	gezelter	2
298	tim	1024	os.flush();
299			#else
300			//every node prepares the dump lines for integrable objects belong to itself
301			std::string buffer;
302			for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
303	xsun	1217
304
305	tim	1024	for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
306			integrableObject = mol->nextIntegrableObject(ii)) {
307	xsun	1217	buffer += prepareDumpLine(integrableObject);
308	gezelter	507	}
309	gezelter	2	}
310	tim	1024
311	gezelter	246	const int masterNode = 0;
312	chuckv	1564	int nProc;
313			MPI_Comm_size(MPI_COMM_WORLD, &nProc);
314	tim	1024	if (worldRank == masterNode) {
315			os << " <Snapshot>\n";
316			writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
317			os << " <StuntDoubles>\n";
318
319	gezelter	1025	os << buffer;
320
321	chuckv	1564
322	tim	1024	for (int i = 1; i < nProc; ++i) {
323	gezelter	2
324	tim	1024	// receive the length of the string buffer that was
325			// prepared by processor i
326	chuckv	1564
327			MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
328	tim	1024	int recvLength;
329			MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
330			char* recvBuffer = new char[recvLength];
331			if (recvBuffer == NULL) {
332			} else {
333	gezelter	1025	MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
334	tim	1024	os << recvBuffer;
335	gezelter	1313	delete [] recvBuffer;
336	tim	1024	}
337			}
338			os << " </StuntDoubles>\n";
339
340			os << " </Snapshot>\n";
341	gezelter	507	os.flush();
342	gezelter	246	} else {
343	gezelter	1025	int sendBufferLength = buffer.size() + 1;
344	chuckv	1564	int myturn = 0;
345			for (int i = 1; i < nProc; ++i){
346			MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
347			if (myturn == worldRank){
348			MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
349			MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
350			}
351			}
352	tim	1024	}
353	gezelter	2
354	tim	1024	#endif // is_mpi
355	gezelter	2
356	tim	1024	}
357	gezelter	2
358	tim	1024	std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
359
360			int index = integrableObject->getGlobalIntegrableObjectIndex();
361			std::string type("pv");
362			std::string line;
363			char tempBuffer[4096];
364	gezelter	2
365	tim	1024	Vector3d pos;
366			Vector3d vel;
367			pos = integrableObject->getPos();
368	gezelter	1437
369			if (isinf(pos[0]) \|\| isnan(pos[0]) \|\|
370			isinf(pos[1]) \|\| isnan(pos[1]) \|\|
371			isinf(pos[2]) \|\| isnan(pos[2]) ) {
372			sprintf( painCave.errMsg,
373			"DumpWriter detected a numerical error writing the position"
374			" for object %d", index);
375			painCave.isFatal = 1;
376			simError();
377			}
378
379	tim	1024	vel = integrableObject->getVel();
380	gezelter	1437
381			if (isinf(vel[0]) \|\| isnan(vel[0]) \|\|
382			isinf(vel[1]) \|\| isnan(vel[1]) \|\|
383			isinf(vel[2]) \|\| isnan(vel[2]) ) {
384			sprintf( painCave.errMsg,
385			"DumpWriter detected a numerical error writing the velocity"
386			" for object %d", index);
387			painCave.isFatal = 1;
388			simError();
389			}
390
391	gezelter	1025	sprintf(tempBuffer, "%18.10g %18.10g %18.10g %13e %13e %13e",
392	tim	1024	pos[0], pos[1], pos[2],
393			vel[0], vel[1], vel[2]);
394			line += tempBuffer;
395	gezelter	2
396	tim	1024	if (integrableObject->isDirectional()) {
397			type += "qj";
398			Quat4d q;
399			Vector3d ji;
400			q = integrableObject->getQ();
401	gezelter	1437
402			if (isinf(q[0]) \|\| isnan(q[0]) \|\|
403			isinf(q[1]) \|\| isnan(q[1]) \|\|
404			isinf(q[2]) \|\| isnan(q[2]) \|\|
405			isinf(q[3]) \|\| isnan(q[3]) ) {
406			sprintf( painCave.errMsg,
407			"DumpWriter detected a numerical error writing the quaternion"
408			" for object %d", index);
409			painCave.isFatal = 1;
410			simError();
411			}
412
413	tim	1024	ji = integrableObject->getJ();
414	gezelter	1437
415			if (isinf(ji[0]) \|\| isnan(ji[0]) \|\|
416			isinf(ji[1]) \|\| isnan(ji[1]) \|\|
417			isinf(ji[2]) \|\| isnan(ji[2]) ) {
418			sprintf( painCave.errMsg,
419			"DumpWriter detected a numerical error writing the angular"
420			" momentum for object %d", index);
421			painCave.isFatal = 1;
422			simError();
423			}
424
425	gezelter	1025	sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e %13e",
426	tim	1024	q[0], q[1], q[2], q[3],
427			ji[0], ji[1], ji[2]);
428			line += tempBuffer;
429			}
430	gezelter	2
431	tim	1024	if (needForceVector_) {
432	gezelter	1437	type += "f";
433	tim	1024	Vector3d frc;
434	gezelter	1437
435	tim	1024	frc = integrableObject->getFrc();
436	gezelter	1437
437			if (isinf(frc[0]) \|\| isnan(frc[0]) \|\|
438			isinf(frc[1]) \|\| isnan(frc[1]) \|\|
439			isinf(frc[2]) \|\| isnan(frc[2]) ) {
440			sprintf( painCave.errMsg,
441			"DumpWriter detected a numerical error writing the force"
442			" for object %d", index);
443			painCave.isFatal = 1;
444			simError();
445			}
446			sprintf(tempBuffer, " %13e %13e %13e",
447			frc[0], frc[1], frc[2]);
448	tim	1024	line += tempBuffer;
449	gezelter	1437
450			if (integrableObject->isDirectional()) {
451			type += "t";
452			Vector3d trq;
453
454			trq = integrableObject->getTrq();
455
456			if (isinf(trq[0]) \|\| isnan(trq[0]) \|\|
457			isinf(trq[1]) \|\| isnan(trq[1]) \|\|
458			isinf(trq[2]) \|\| isnan(trq[2]) ) {
459			sprintf( painCave.errMsg,
460			"DumpWriter detected a numerical error writing the torque"
461			" for object %d", index);
462			painCave.isFatal = 1;
463			simError();
464			}
465
466			sprintf(tempBuffer, " %13e %13e %13e",
467			trq[0], trq[1], trq[2]);
468			line += tempBuffer;
469			}
470	gezelter	246	}
471	gezelter	1610	if (needParticlePot_) {
472			type += "u";
473			RealType particlePot;
474
475			particlePot = integrableObject->getParticlePot();
476
477			if (isinf(particlePot) \|\| isnan(particlePot)) {
478			sprintf( painCave.errMsg,
479			"DumpWriter detected a numerical error writing the particle "
480			" potential for object %d", index);
481			painCave.isFatal = 1;
482			simError();
483			}
484			sprintf(tempBuffer, " %13e", particlePot);
485			line += tempBuffer;
486			}
487	gezelter	1437
488	gezelter	1025	sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
489	tim	1024	return std::string(tempBuffer);
490	gezelter	507	}
491	gezelter	2
492	gezelter	507	void DumpWriter::writeDump() {
493	tim	615	writeFrame(*dumpFile_);
494	gezelter	507	}
495	tim	376
496	gezelter	507	void DumpWriter::writeEor() {
497	tim	615	std::ostream* eorStream;
498	tim	376
499			#ifdef IS_MPI
500			if (worldRank == 0) {
501			#endif // is_mpi
502
503	tim	615	eorStream = createOStream(eorFilename_);
504	tim	376
505			#ifdef IS_MPI
506			}
507			#endif // is_mpi
508
509	tim	615	writeFrame(*eorStream);
510
511			#ifdef IS_MPI
512			if (worldRank == 0) {
513			#endif // is_mpi
514	tim	1024	writeClosing(*eorStream);
515			delete eorStream;
516	tim	615	#ifdef IS_MPI
517			}
518			#endif // is_mpi
519
520	gezelter	507	}
521	tim	376
522
523	gezelter	507	void DumpWriter::writeDumpAndEor() {
524	tim	376	std::vector<std::streambuf*> buffers;
525	tim	615	std::ostream* eorStream;
526	tim	376	#ifdef IS_MPI
527			if (worldRank == 0) {
528			#endif // is_mpi
529
530	tim	615	buffers.push_back(dumpFile_->rdbuf());
531	tim	376
532	tim	615	eorStream = createOStream(eorFilename_);
533	tim	376
534	tim	615	buffers.push_back(eorStream->rdbuf());
535	tim	376
536			#ifdef IS_MPI
537			}
538			#endif // is_mpi
539
540			TeeBuf tbuf(buffers.begin(), buffers.end());
541			std::ostream os(&tbuf);
542
543			writeFrame(os);
544	tim	615
545			#ifdef IS_MPI
546			if (worldRank == 0) {
547			#endif // is_mpi
548	tim	1024	writeClosing(*eorStream);
549			delete eorStream;
550	tim	615	#ifdef IS_MPI
551			}
552			#endif // is_mpi
553	tim	376
554	gezelter	507	}
555	tim	376
556	tim	1024	std::ostream* DumpWriter::createOStream(const std::string& filename) {
557	tim	619
558	tim	615	std::ostream* newOStream;
559	tim	619	#ifdef HAVE_LIBZ
560	tim	615	if (needCompression_) {
561	tim	1024	newOStream = new ogzstream(filename.c_str());
562	tim	615	} else {
563	tim	1024	newOStream = new std::ofstream(filename.c_str());
564	tim	615	}
565	tim	619	#else
566			newOStream = new std::ofstream(filename.c_str());
567			#endif
568	tim	1024	//write out MetaData first
569	gezelter	1390	(*newOStream) << "<OpenMD version=1>" << std::endl;
570	tim	1024	(*newOStream) << " <MetaData>" << std::endl;
571			(*newOStream) << info_->getRawMetaData();
572			(*newOStream) << " </MetaData>" << std::endl;
573	tim	615	return newOStream;
574	tim	1024	}
575	tim	376
576	tim	1024	void DumpWriter::writeClosing(std::ostream& os) {
577
578	gezelter	1390	os << "</OpenMD>\n";
579	tim	1024	os.flush();
580			}
581
582	gezelter	1390	}//end namespace OpenMD