ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/io/DumpWriter.cpp
Revision: 1880
Committed: Mon Jun 17 18:28:30 2013 UTC (11 years, 10 months ago) by gezelter
File size: 23848 byte(s)
Log Message:
Preparing for official 2.1 release (clean-up)

File Contents

# 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 gezelter 1879 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39 gezelter 1782 * [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 gezelter 1782 #ifdef HAVE_ZLIB
48 tim 615 #include "io/gzstream.hpp"
49 gezelter 1782 #endif
50 tim 615 #include "io/Globals.hpp"
51    
52 gezelter 1782 #ifdef _MSC_VER
53     #define isnan(x) _isnan((x))
54     #define isinf(x) (!_finite(x) && !_isnan(x))
55     #endif
56 gezelter 1437
57 gezelter 2 #ifdef IS_MPI
58     #include <mpi.h>
59 gezelter 1782 #endif
60 gezelter 2
61 gezelter 1437 using namespace std;
62 gezelter 1390 namespace OpenMD {
63 gezelter 2
64 gezelter 507 DumpWriter::DumpWriter(SimInfo* info)
65     : info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
66 tim 615
67     Globals* simParams = info->getSimParams();
68 gezelter 1782 needCompression_ = simParams->getCompressDumpFile();
69     needForceVector_ = simParams->getOutputForceVector();
70     needParticlePot_ = simParams->getOutputParticlePotential();
71     needFlucQ_ = simParams->getOutputFluctuatingCharges();
72     needElectricField_ = simParams->getOutputElectricField();
73    
74     if (needParticlePot_ || needFlucQ_ || needElectricField_) {
75     doSiteData_ = true;
76     } else {
77     doSiteData_ = false;
78     }
79    
80 chuckv 791 createDumpFile_ = true;
81 tim 619 #ifdef HAVE_LIBZ
82 tim 615 if (needCompression_) {
83 tim 1024 filename_ += ".gz";
84     eorFilename_ += ".gz";
85 tim 615 }
86 tim 619 #endif
87 tim 615
88 tim 376 #ifdef IS_MPI
89    
90 tim 1024 if (worldRank == 0) {
91 tim 376 #endif // is_mpi
92 chuckv 791
93 tim 1024 dumpFile_ = createOStream(filename_);
94 tim 615
95 tim 1024 if (!dumpFile_) {
96     sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
97     filename_.c_str());
98     painCave.isFatal = 1;
99     simError();
100     }
101 tim 376
102     #ifdef IS_MPI
103    
104 tim 1024 }
105 tim 376
106     #endif // is_mpi
107    
108 tim 1024 }
109 tim 376
110    
111 gezelter 507 DumpWriter::DumpWriter(SimInfo* info, const std::string& filename)
112     : info_(info), filename_(filename){
113 tim 615
114     Globals* simParams = info->getSimParams();
115     eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
116    
117 gezelter 1782 needCompression_ = simParams->getCompressDumpFile();
118     needForceVector_ = simParams->getOutputForceVector();
119     needParticlePot_ = simParams->getOutputParticlePotential();
120     needFlucQ_ = simParams->getOutputFluctuatingCharges();
121     needElectricField_ = simParams->getOutputElectricField();
122    
123     if (needParticlePot_ || needFlucQ_ || needElectricField_) {
124     doSiteData_ = true;
125     } else {
126     doSiteData_ = false;
127     }
128    
129 chuckv 791 createDumpFile_ = true;
130 tim 619 #ifdef HAVE_LIBZ
131 tim 615 if (needCompression_) {
132 tim 1024 filename_ += ".gz";
133     eorFilename_ += ".gz";
134 tim 615 }
135 tim 619 #endif
136 tim 615
137 gezelter 246 #ifdef IS_MPI
138 gezelter 2
139 tim 1024 if (worldRank == 0) {
140 gezelter 246 #endif // is_mpi
141 gezelter 2
142 chuckv 791
143 tim 1024 dumpFile_ = createOStream(filename_);
144 tim 615
145 tim 1024 if (!dumpFile_) {
146     sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
147     filename_.c_str());
148     painCave.isFatal = 1;
149     simError();
150     }
151 gezelter 2
152     #ifdef IS_MPI
153    
154 tim 1024 }
155 gezelter 2
156     #endif // is_mpi
157 gezelter 246
158 tim 1024 }
159 chuckv 791
160     DumpWriter::DumpWriter(SimInfo* info, const std::string& filename, bool writeDumpFile)
161 tim 1024 : info_(info), filename_(filename){
162 chuckv 791
163     Globals* simParams = info->getSimParams();
164     eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
165    
166 gezelter 1782 needCompression_ = simParams->getCompressDumpFile();
167     needForceVector_ = simParams->getOutputForceVector();
168     needParticlePot_ = simParams->getOutputParticlePotential();
169     needFlucQ_ = simParams->getOutputFluctuatingCharges();
170     needElectricField_ = simParams->getOutputElectricField();
171    
172     if (needParticlePot_ || needFlucQ_ || needElectricField_) {
173     doSiteData_ = true;
174     } else {
175     doSiteData_ = false;
176     }
177    
178 chuckv 791 #ifdef HAVE_LIBZ
179     if (needCompression_) {
180     filename_ += ".gz";
181     eorFilename_ += ".gz";
182     }
183     #endif
184    
185     #ifdef IS_MPI
186    
187     if (worldRank == 0) {
188     #endif // is_mpi
189    
190     createDumpFile_ = writeDumpFile;
191     if (createDumpFile_) {
192     dumpFile_ = createOStream(filename_);
193    
194     if (!dumpFile_) {
195     sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
196     filename_.c_str());
197     painCave.isFatal = 1;
198     simError();
199     }
200     }
201     #ifdef IS_MPI
202    
203     }
204 tim 1024
205 chuckv 791
206     #endif // is_mpi
207    
208     }
209 gezelter 2
210 gezelter 507 DumpWriter::~DumpWriter() {
211 gezelter 2
212     #ifdef IS_MPI
213 gezelter 246
214     if (worldRank == 0) {
215 gezelter 2 #endif // is_mpi
216 chuckv 791 if (createDumpFile_){
217 tim 1024 writeClosing(*dumpFile_);
218 chuckv 791 delete dumpFile_;
219     }
220 gezelter 2 #ifdef IS_MPI
221 gezelter 246
222     }
223    
224 gezelter 2 #endif // is_mpi
225 gezelter 246
226 gezelter 507 }
227 gezelter 2
228 tim 1024 void DumpWriter::writeFrameProperties(std::ostream& os, Snapshot* s) {
229 gezelter 2
230 tim 1024 char buffer[1024];
231    
232     os << " <FrameData>\n";
233    
234     RealType currentTime = s->getTime();
235 gezelter 1437
236     if (isinf(currentTime) || isnan(currentTime)) {
237     sprintf( painCave.errMsg,
238     "DumpWriter detected a numerical error writing the time");
239     painCave.isFatal = 1;
240     simError();
241     }
242    
243 gezelter 1025 sprintf(buffer, " Time: %.10g\n", currentTime);
244 tim 1024 os << buffer;
245    
246 gezelter 246 Mat3x3d hmat;
247     hmat = s->getHmat();
248 gezelter 1437
249     for (unsigned int i = 0; i < 3; i++) {
250     for (unsigned int j = 0; j < 3; j++) {
251     if (isinf(hmat(i,j)) || isnan(hmat(i,j))) {
252     sprintf( painCave.errMsg,
253     "DumpWriter detected a numerical error writing the box");
254     painCave.isFatal = 1;
255     simError();
256     }
257     }
258     }
259    
260 tim 1024 sprintf(buffer, " Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
261     hmat(0, 0), hmat(1, 0), hmat(2, 0),
262     hmat(0, 1), hmat(1, 1), hmat(2, 1),
263     hmat(0, 2), hmat(1, 2), hmat(2, 2));
264     os << buffer;
265 gezelter 2
266 gezelter 1782 pair<RealType, RealType> thermostat = s->getThermostat();
267    
268     if (isinf(thermostat.first) || isnan(thermostat.first) ||
269     isinf(thermostat.second) || isnan(thermostat.second)) {
270 gezelter 1437 sprintf( painCave.errMsg,
271     "DumpWriter detected a numerical error writing the thermostat");
272     painCave.isFatal = 1;
273     simError();
274     }
275 gezelter 1782 sprintf(buffer, " Thermostat: %.10g , %.10g\n", thermostat.first,
276     thermostat.second);
277 tim 1024 os << buffer;
278 gezelter 246
279 tim 1024 Mat3x3d eta;
280 gezelter 1782 eta = s->getBarostat();
281 gezelter 1437
282     for (unsigned int i = 0; i < 3; i++) {
283     for (unsigned int j = 0; j < 3; j++) {
284     if (isinf(eta(i,j)) || isnan(eta(i,j))) {
285     sprintf( painCave.errMsg,
286     "DumpWriter detected a numerical error writing the barostat");
287     painCave.isFatal = 1;
288     simError();
289     }
290     }
291     }
292    
293 tim 1024 sprintf(buffer, " Barostat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
294     eta(0, 0), eta(1, 0), eta(2, 0),
295     eta(0, 1), eta(1, 1), eta(2, 1),
296     eta(0, 2), eta(1, 2), eta(2, 2));
297     os << buffer;
298 gezelter 246
299 tim 1024 os << " </FrameData>\n";
300 gezelter 507 }
301 gezelter 2
302 gezelter 507 void DumpWriter::writeFrame(std::ostream& os) {
303 gezelter 246
304 tim 1024 #ifdef IS_MPI
305 gezelter 1796 MPI::Status istatus;
306 tim 1024 #endif
307 gezelter 246
308     Molecule* mol;
309 gezelter 1782 StuntDouble* sd;
310 gezelter 246 SimInfo::MoleculeIterator mi;
311     Molecule::IntegrableObjectIterator ii;
312 gezelter 1782 RigidBody::AtomIterator ai;
313 gezelter 2
314 gezelter 246 #ifndef IS_MPI
315 tim 1024 os << " <Snapshot>\n";
316 gezelter 1025
317 tim 1024 writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
318 gezelter 2
319 tim 1024 os << " <StuntDoubles>\n";
320 gezelter 1879 for (mol = info_->beginMolecule(mi); mol != NULL;
321     mol = info_->nextMolecule(mi)) {
322 xsun 1217
323 gezelter 1782 for (sd = mol->beginIntegrableObject(ii); sd != NULL;
324     sd = mol->nextIntegrableObject(ii)) {
325     os << prepareDumpLine(sd);
326 xsun 1217
327 tim 1024 }
328     }
329     os << " </StuntDoubles>\n";
330 gezelter 1782
331     if (doSiteData_) {
332     os << " <SiteData>\n";
333 gezelter 1796 for (mol = info_->beginMolecule(mi); mol != NULL;
334     mol = info_->nextMolecule(mi)) {
335 gezelter 1782
336     for (sd = mol->beginIntegrableObject(ii); sd != NULL;
337 gezelter 1880 sd = mol->nextIntegrableObject(ii)) {
338    
339 gezelter 1782 int ioIndex = sd->getGlobalIntegrableObjectIndex();
340     // do one for the IO itself
341     os << prepareSiteLine(sd, ioIndex, 0);
342    
343     if (sd->isRigidBody()) {
344    
345     RigidBody* rb = static_cast<RigidBody*>(sd);
346     int siteIndex = 0;
347 gezelter 1879 for (Atom* atom = rb->beginAtom(ai); atom != NULL;
348 gezelter 1782 atom = rb->nextAtom(ai)) {
349     os << prepareSiteLine(atom, ioIndex, siteIndex);
350     siteIndex++;
351     }
352     }
353     }
354     }
355     os << " </SiteData>\n";
356     }
357 tim 1024 os << " </Snapshot>\n";
358 gezelter 2
359 tim 1024 os.flush();
360     #else
361 gezelter 1796
362     const int masterNode = 0;
363     int worldRank = MPI::COMM_WORLD.Get_rank();
364     int nProc = MPI::COMM_WORLD.Get_size();
365    
366     if (worldRank == masterNode) {
367     os << " <Snapshot>\n";
368     writeFrameProperties(os,
369     info_->getSnapshotManager()->getCurrentSnapshot());
370     os << " <StuntDoubles>\n";
371     }
372    
373 tim 1024 //every node prepares the dump lines for integrable objects belong to itself
374     std::string buffer;
375 gezelter 1796 for (mol = info_->beginMolecule(mi); mol != NULL;
376     mol = info_->nextMolecule(mi)) {
377 gezelter 1782 for (sd = mol->beginIntegrableObject(ii); sd != NULL;
378     sd = mol->nextIntegrableObject(ii)) {
379 gezelter 1796 buffer += prepareDumpLine(sd);
380 gezelter 507 }
381 gezelter 2 }
382 tim 1024
383     if (worldRank == masterNode) {
384 gezelter 1025 os << buffer;
385 gezelter 1796
386 tim 1024 for (int i = 1; i < nProc; ++i) {
387 gezelter 1796 // tell processor i to start sending us data:
388     MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode);
389 gezelter 2
390 tim 1024 // receive the length of the string buffer that was
391 gezelter 1796 // prepared by processor i:
392     int recvLength;
393     MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
394     istatus);
395 gezelter 1782
396 gezelter 1796 // create a buffer to receive the data
397 tim 1024 char* recvBuffer = new char[recvLength];
398     if (recvBuffer == NULL) {
399     } else {
400 gezelter 1796 // receive the data:
401     MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i,
402     MPI::ANY_TAG, istatus);
403     // send it to the file:
404 tim 1024 os << recvBuffer;
405 gezelter 1796 // get rid of the receive buffer:
406 gezelter 1313 delete [] recvBuffer;
407 tim 1024 }
408     }
409 gezelter 246 } else {
410 gezelter 1025 int sendBufferLength = buffer.size() + 1;
411 chuckv 1564 int myturn = 0;
412     for (int i = 1; i < nProc; ++i){
413 gezelter 1796 // wait for the master node to call our number:
414     MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
415 chuckv 1564 if (myturn == worldRank){
416 gezelter 1796 // send the length of our buffer:
417     MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0);
418    
419     // send our buffer:
420     MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength,
421     MPI::CHAR, masterNode, 0);
422 gezelter 1879
423 chuckv 1564 }
424     }
425 tim 1024 }
426 gezelter 1796
427     if (worldRank == masterNode) {
428     os << " </StuntDoubles>\n";
429     }
430 gezelter 2
431 gezelter 1796 if (doSiteData_) {
432     if (worldRank == masterNode) {
433     os << " <SiteData>\n";
434     }
435     buffer.clear();
436     for (mol = info_->beginMolecule(mi); mol != NULL;
437     mol = info_->nextMolecule(mi)) {
438    
439     for (sd = mol->beginIntegrableObject(ii); sd != NULL;
440     sd = mol->nextIntegrableObject(ii)) {
441    
442     int ioIndex = sd->getGlobalIntegrableObjectIndex();
443     // do one for the IO itself
444     buffer += prepareSiteLine(sd, ioIndex, 0);
445    
446     if (sd->isRigidBody()) {
447    
448     RigidBody* rb = static_cast<RigidBody*>(sd);
449     int siteIndex = 0;
450 gezelter 1879 for (Atom* atom = rb->beginAtom(ai); atom != NULL;
451 gezelter 1796 atom = rb->nextAtom(ai)) {
452     buffer += prepareSiteLine(atom, ioIndex, siteIndex);
453     siteIndex++;
454     }
455     }
456     }
457     }
458    
459     if (worldRank == masterNode) {
460     os << buffer;
461    
462     for (int i = 1; i < nProc; ++i) {
463    
464     // tell processor i to start sending us data:
465     MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode);
466    
467     // receive the length of the string buffer that was
468     // prepared by processor i:
469     int recvLength;
470     MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
471     istatus);
472    
473     // create a buffer to receive the data
474     char* recvBuffer = new char[recvLength];
475     if (recvBuffer == NULL) {
476     } else {
477     // receive the data:
478     MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i,
479     MPI::ANY_TAG, istatus);
480     // send it to the file:
481     os << recvBuffer;
482     // get rid of the receive buffer:
483     delete [] recvBuffer;
484     }
485     }
486     } else {
487     int sendBufferLength = buffer.size() + 1;
488     int myturn = 0;
489     for (int i = 1; i < nProc; ++i){
490     // wait for the master node to call our number:
491     MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
492     if (myturn == worldRank){
493     // send the length of our buffer:
494     MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0);
495     // send our buffer:
496     MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength,
497     MPI::CHAR, masterNode, 0);
498     }
499     }
500     }
501    
502     if (worldRank == masterNode) {
503     os << " </SiteData>\n";
504     }
505     }
506    
507     if (worldRank == masterNode) {
508     os << " </Snapshot>\n";
509     os.flush();
510     }
511    
512 tim 1024 #endif // is_mpi
513 gezelter 1796
514 tim 1024 }
515 gezelter 2
516 gezelter 1782 std::string DumpWriter::prepareDumpLine(StuntDouble* sd) {
517 tim 1024
518 gezelter 1782 int index = sd->getGlobalIntegrableObjectIndex();
519 tim 1024 std::string type("pv");
520     std::string line;
521     char tempBuffer[4096];
522 gezelter 2
523 tim 1024 Vector3d pos;
524     Vector3d vel;
525 gezelter 1782 pos = sd->getPos();
526 gezelter 1437
527     if (isinf(pos[0]) || isnan(pos[0]) ||
528     isinf(pos[1]) || isnan(pos[1]) ||
529     isinf(pos[2]) || isnan(pos[2]) ) {
530     sprintf( painCave.errMsg,
531     "DumpWriter detected a numerical error writing the position"
532     " for object %d", index);
533     painCave.isFatal = 1;
534     simError();
535     }
536    
537 gezelter 1782 vel = sd->getVel();
538 gezelter 1437
539     if (isinf(vel[0]) || isnan(vel[0]) ||
540     isinf(vel[1]) || isnan(vel[1]) ||
541     isinf(vel[2]) || isnan(vel[2]) ) {
542     sprintf( painCave.errMsg,
543     "DumpWriter detected a numerical error writing the velocity"
544     " for object %d", index);
545     painCave.isFatal = 1;
546     simError();
547     }
548    
549 gezelter 1025 sprintf(tempBuffer, "%18.10g %18.10g %18.10g %13e %13e %13e",
550 tim 1024 pos[0], pos[1], pos[2],
551     vel[0], vel[1], vel[2]);
552     line += tempBuffer;
553 gezelter 2
554 gezelter 1782 if (sd->isDirectional()) {
555 tim 1024 type += "qj";
556     Quat4d q;
557     Vector3d ji;
558 gezelter 1782 q = sd->getQ();
559 gezelter 1437
560     if (isinf(q[0]) || isnan(q[0]) ||
561     isinf(q[1]) || isnan(q[1]) ||
562     isinf(q[2]) || isnan(q[2]) ||
563     isinf(q[3]) || isnan(q[3]) ) {
564     sprintf( painCave.errMsg,
565     "DumpWriter detected a numerical error writing the quaternion"
566     " for object %d", index);
567     painCave.isFatal = 1;
568     simError();
569     }
570    
571 gezelter 1782 ji = sd->getJ();
572 gezelter 1437
573     if (isinf(ji[0]) || isnan(ji[0]) ||
574     isinf(ji[1]) || isnan(ji[1]) ||
575     isinf(ji[2]) || isnan(ji[2]) ) {
576     sprintf( painCave.errMsg,
577     "DumpWriter detected a numerical error writing the angular"
578     " momentum for object %d", index);
579     painCave.isFatal = 1;
580     simError();
581     }
582    
583 gezelter 1025 sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e %13e",
584 tim 1024 q[0], q[1], q[2], q[3],
585     ji[0], ji[1], ji[2]);
586     line += tempBuffer;
587     }
588 gezelter 2
589 tim 1024 if (needForceVector_) {
590 gezelter 1437 type += "f";
591 gezelter 1782 Vector3d frc = sd->getFrc();
592 gezelter 1437 if (isinf(frc[0]) || isnan(frc[0]) ||
593     isinf(frc[1]) || isnan(frc[1]) ||
594     isinf(frc[2]) || isnan(frc[2]) ) {
595     sprintf( painCave.errMsg,
596     "DumpWriter detected a numerical error writing the force"
597     " for object %d", index);
598     painCave.isFatal = 1;
599     simError();
600     }
601     sprintf(tempBuffer, " %13e %13e %13e",
602     frc[0], frc[1], frc[2]);
603 tim 1024 line += tempBuffer;
604 gezelter 1437
605 gezelter 1782 if (sd->isDirectional()) {
606 gezelter 1437 type += "t";
607 gezelter 1782 Vector3d trq = sd->getTrq();
608 gezelter 1437 if (isinf(trq[0]) || isnan(trq[0]) ||
609     isinf(trq[1]) || isnan(trq[1]) ||
610     isinf(trq[2]) || isnan(trq[2]) ) {
611     sprintf( painCave.errMsg,
612     "DumpWriter detected a numerical error writing the torque"
613     " for object %d", index);
614     painCave.isFatal = 1;
615     simError();
616 gezelter 1782 }
617 gezelter 1437 sprintf(tempBuffer, " %13e %13e %13e",
618     trq[0], trq[1], trq[2]);
619     line += tempBuffer;
620     }
621 gezelter 246 }
622 gezelter 1782
623     sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
624     return std::string(tempBuffer);
625     }
626    
627     std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
628    
629    
630     std::string id;
631     std::string type;
632     std::string line;
633     char tempBuffer[4096];
634    
635     if (sd->isRigidBody()) {
636     sprintf(tempBuffer, "%10d ", ioIndex);
637     id = std::string(tempBuffer);
638     } else {
639     sprintf(tempBuffer, "%10d %10d", ioIndex, siteIndex);
640     id = std::string(tempBuffer);
641     }
642    
643     if (needFlucQ_) {
644     type += "cw";
645     RealType fqPos = sd->getFlucQPos();
646     if (isinf(fqPos) || isnan(fqPos) ) {
647     sprintf( painCave.errMsg,
648     "DumpWriter detected a numerical error writing the"
649     " fluctuating charge for object %s", id.c_str());
650     painCave.isFatal = 1;
651     simError();
652     }
653     sprintf(tempBuffer, " %13e ", fqPos);
654     line += tempBuffer;
655    
656     RealType fqVel = sd->getFlucQVel();
657     if (isinf(fqVel) || isnan(fqVel) ) {
658     sprintf( painCave.errMsg,
659     "DumpWriter detected a numerical error writing the"
660     " fluctuating charge velocity for object %s", id.c_str());
661     painCave.isFatal = 1;
662     simError();
663     }
664     sprintf(tempBuffer, " %13e ", fqVel);
665     line += tempBuffer;
666    
667     if (needForceVector_) {
668     type += "g";
669     RealType fqFrc = sd->getFlucQFrc();
670     if (isinf(fqFrc) || isnan(fqFrc) ) {
671     sprintf( painCave.errMsg,
672     "DumpWriter detected a numerical error writing the"
673     " fluctuating charge force for object %s", id.c_str());
674     painCave.isFatal = 1;
675     simError();
676     }
677     sprintf(tempBuffer, " %13e ", fqFrc);
678     line += tempBuffer;
679     }
680     }
681    
682     if (needElectricField_) {
683     type += "e";
684     Vector3d eField= sd->getElectricField();
685     if (isinf(eField[0]) || isnan(eField[0]) ||
686     isinf(eField[1]) || isnan(eField[1]) ||
687     isinf(eField[2]) || isnan(eField[2]) ) {
688     sprintf( painCave.errMsg,
689     "DumpWriter detected a numerical error writing the electric"
690     " field for object %s", id.c_str());
691     painCave.isFatal = 1;
692     simError();
693     }
694     sprintf(tempBuffer, " %13e %13e %13e",
695     eField[0], eField[1], eField[2]);
696     line += tempBuffer;
697     }
698    
699    
700 gezelter 1610 if (needParticlePot_) {
701     type += "u";
702 gezelter 1782 RealType particlePot = sd->getParticlePot();
703 gezelter 1610 if (isinf(particlePot) || isnan(particlePot)) {
704     sprintf( painCave.errMsg,
705     "DumpWriter detected a numerical error writing the particle "
706 gezelter 1782 " potential for object %s", id.c_str());
707 gezelter 1610 painCave.isFatal = 1;
708     simError();
709     }
710     sprintf(tempBuffer, " %13e", particlePot);
711     line += tempBuffer;
712     }
713 gezelter 1437
714 gezelter 1782
715     sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
716 tim 1024 return std::string(tempBuffer);
717 gezelter 507 }
718 gezelter 2
719 gezelter 507 void DumpWriter::writeDump() {
720 tim 615 writeFrame(*dumpFile_);
721 gezelter 507 }
722 tim 376
723 gezelter 507 void DumpWriter::writeEor() {
724 tim 615 std::ostream* eorStream;
725 tim 376
726     #ifdef IS_MPI
727     if (worldRank == 0) {
728     #endif // is_mpi
729 gezelter 1879
730 tim 615 eorStream = createOStream(eorFilename_);
731 tim 376
732     #ifdef IS_MPI
733     }
734 gezelter 1879 #endif
735    
736 tim 615 writeFrame(*eorStream);
737 gezelter 1879
738 tim 615 #ifdef IS_MPI
739     if (worldRank == 0) {
740 gezelter 1879 #endif
741    
742 tim 1024 writeClosing(*eorStream);
743     delete eorStream;
744 gezelter 1879
745 tim 615 #ifdef IS_MPI
746     }
747     #endif // is_mpi
748    
749 gezelter 507 }
750 tim 376
751    
752 gezelter 507 void DumpWriter::writeDumpAndEor() {
753 tim 376 std::vector<std::streambuf*> buffers;
754 tim 615 std::ostream* eorStream;
755 tim 376 #ifdef IS_MPI
756     if (worldRank == 0) {
757     #endif // is_mpi
758    
759 tim 615 buffers.push_back(dumpFile_->rdbuf());
760 tim 376
761 tim 615 eorStream = createOStream(eorFilename_);
762 tim 376
763 tim 615 buffers.push_back(eorStream->rdbuf());
764 tim 376
765     #ifdef IS_MPI
766     }
767     #endif // is_mpi
768    
769     TeeBuf tbuf(buffers.begin(), buffers.end());
770     std::ostream os(&tbuf);
771    
772     writeFrame(os);
773 tim 615
774     #ifdef IS_MPI
775     if (worldRank == 0) {
776     #endif // is_mpi
777 tim 1024 writeClosing(*eorStream);
778     delete eorStream;
779 tim 615 #ifdef IS_MPI
780     }
781     #endif // is_mpi
782 tim 376
783 gezelter 507 }
784 tim 376
785 tim 1024 std::ostream* DumpWriter::createOStream(const std::string& filename) {
786 tim 619
787 tim 615 std::ostream* newOStream;
788 gezelter 1782 #ifdef HAVE_ZLIB
789 tim 615 if (needCompression_) {
790 tim 1024 newOStream = new ogzstream(filename.c_str());
791 tim 615 } else {
792 tim 1024 newOStream = new std::ofstream(filename.c_str());
793 tim 615 }
794 tim 619 #else
795     newOStream = new std::ofstream(filename.c_str());
796     #endif
797 tim 1024 //write out MetaData first
798 gezelter 1782 (*newOStream) << "<OpenMD version=2>" << std::endl;
799 tim 1024 (*newOStream) << " <MetaData>" << std::endl;
800     (*newOStream) << info_->getRawMetaData();
801     (*newOStream) << " </MetaData>" << std::endl;
802 tim 615 return newOStream;
803 tim 1024 }
804 tim 376
805 tim 1024 void DumpWriter::writeClosing(std::ostream& os) {
806    
807 gezelter 1390 os << "</OpenMD>\n";
808 tim 1024 os.flush();
809     }
810    
811 gezelter 1390 }//end namespace OpenMD

Properties

Name Value
svn:keywords Author Id Revision Date