ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/io/DumpWriter.cpp
Revision: 1875
Committed: Fri May 17 14:41:42 2013 UTC (11 years, 11 months ago) by gezelter
File size: 23820 byte(s)
Log Message:
Fixed a bunch of stylistic and performance issues discovered via cppcheck.

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 1850 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (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 gezelter 1767 #ifdef HAVE_ZLIB
48 tim 615 #include "io/gzstream.hpp"
49 gezelter 1767 #endif
50 tim 615 #include "io/Globals.hpp"
51    
52 gezelter 1767 #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 1764 #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 1714 needCompression_ = simParams->getCompressDumpFile();
69     needForceVector_ = simParams->getOutputForceVector();
70     needParticlePot_ = simParams->getOutputParticlePotential();
71     needFlucQ_ = simParams->getOutputFluctuatingCharges();
72     needElectricField_ = simParams->getOutputElectricField();
73    
74 gezelter 1752 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 1714 needCompression_ = simParams->getCompressDumpFile();
118     needForceVector_ = simParams->getOutputForceVector();
119     needParticlePot_ = simParams->getOutputParticlePotential();
120     needFlucQ_ = simParams->getOutputFluctuatingCharges();
121     needElectricField_ = simParams->getOutputElectricField();
122    
123 gezelter 1752 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 1714 needCompression_ = simParams->getCompressDumpFile();
167     needForceVector_ = simParams->getOutputForceVector();
168     needParticlePot_ = simParams->getOutputParticlePotential();
169     needFlucQ_ = simParams->getOutputFluctuatingCharges();
170     needElectricField_ = simParams->getOutputElectricField();
171    
172 gezelter 1752 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 1764 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 1764 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 1764 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 1798 MPI::Status istatus;
306 tim 1024 #endif
307 gezelter 246
308     Molecule* mol;
309 gezelter 1769 StuntDouble* sd;
310 gezelter 246 SimInfo::MoleculeIterator mi;
311     Molecule::IntegrableObjectIterator ii;
312 gezelter 1752 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 1875 for (mol = info_->beginMolecule(mi); mol != NULL;
321     mol = info_->nextMolecule(mi)) {
322 xsun 1217
323 gezelter 1769 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 1752
331     if (doSiteData_) {
332     os << " <SiteData>\n";
333 gezelter 1798 for (mol = info_->beginMolecule(mi); mol != NULL;
334     mol = info_->nextMolecule(mi)) {
335 gezelter 1752
336 gezelter 1769 for (sd = mol->beginIntegrableObject(ii); sd != NULL;
337     sd = mol->nextIntegrableObject(ii)) {
338 gezelter 1752
339 gezelter 1769 int ioIndex = sd->getGlobalIntegrableObjectIndex();
340 gezelter 1752 // do one for the IO itself
341 gezelter 1769 os << prepareSiteLine(sd, ioIndex, 0);
342 gezelter 1752
343 gezelter 1769 if (sd->isRigidBody()) {
344 gezelter 1752
345 gezelter 1769 RigidBody* rb = static_cast<RigidBody*>(sd);
346 gezelter 1752 int siteIndex = 0;
347 gezelter 1875 for (Atom* atom = rb->beginAtom(ai); atom != NULL;
348 gezelter 1752 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 1798
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 1798 for (mol = info_->beginMolecule(mi); mol != NULL;
376     mol = info_->nextMolecule(mi)) {
377 gezelter 1769 for (sd = mol->beginIntegrableObject(ii); sd != NULL;
378     sd = mol->nextIntegrableObject(ii)) {
379 gezelter 1798 buffer += prepareDumpLine(sd);
380 gezelter 507 }
381 gezelter 2 }
382 tim 1024
383     if (worldRank == masterNode) {
384 gezelter 1025 os << buffer;
385 gezelter 1798
386 tim 1024 for (int i = 1; i < nProc; ++i) {
387 gezelter 1798 // 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 1798 // prepared by processor i:
392     int recvLength;
393     MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
394     istatus);
395 gezelter 2
396 gezelter 1798 // create a buffer to receive the data
397 tim 1024 char* recvBuffer = new char[recvLength];
398     if (recvBuffer == NULL) {
399     } else {
400 gezelter 1798 // 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 1798 // 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 gezelter 1629 int myturn = 0;
412     for (int i = 1; i < nProc; ++i){
413 gezelter 1798 // wait for the master node to call our number:
414     MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
415 gezelter 1629 if (myturn == worldRank){
416 gezelter 1798 // 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 1629 }
423     }
424 tim 1024 }
425 gezelter 1798
426     if (worldRank == masterNode) {
427     os << " </StuntDoubles>\n";
428     }
429 gezelter 2
430 gezelter 1798 if (doSiteData_) {
431     if (worldRank == masterNode) {
432     os << " <SiteData>\n";
433     }
434     buffer.clear();
435     for (mol = info_->beginMolecule(mi); mol != NULL;
436     mol = info_->nextMolecule(mi)) {
437    
438     for (sd = mol->beginIntegrableObject(ii); sd != NULL;
439     sd = mol->nextIntegrableObject(ii)) {
440    
441     int ioIndex = sd->getGlobalIntegrableObjectIndex();
442     // do one for the IO itself
443     buffer += prepareSiteLine(sd, ioIndex, 0);
444    
445     if (sd->isRigidBody()) {
446    
447     RigidBody* rb = static_cast<RigidBody*>(sd);
448     int siteIndex = 0;
449 gezelter 1875 for (Atom* atom = rb->beginAtom(ai); atom != NULL;
450 gezelter 1798 atom = rb->nextAtom(ai)) {
451     buffer += prepareSiteLine(atom, ioIndex, siteIndex);
452     siteIndex++;
453     }
454     }
455     }
456     }
457    
458     if (worldRank == masterNode) {
459     os << buffer;
460    
461     for (int i = 1; i < nProc; ++i) {
462    
463     // tell processor i to start sending us data:
464     MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode);
465    
466     // receive the length of the string buffer that was
467     // prepared by processor i:
468     int recvLength;
469     MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
470     istatus);
471    
472     // create a buffer to receive the data
473     char* recvBuffer = new char[recvLength];
474     if (recvBuffer == NULL) {
475     } else {
476     // receive the data:
477     MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i,
478     MPI::ANY_TAG, istatus);
479     // send it to the file:
480     os << recvBuffer;
481     // get rid of the receive buffer:
482     delete [] recvBuffer;
483     }
484     }
485     } else {
486     int sendBufferLength = buffer.size() + 1;
487     int myturn = 0;
488     for (int i = 1; i < nProc; ++i){
489     // wait for the master node to call our number:
490     MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
491     if (myturn == worldRank){
492     // send the length of our buffer:
493     MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0);
494     // send our buffer:
495     MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength,
496     MPI::CHAR, masterNode, 0);
497     }
498     }
499     }
500    
501     if (worldRank == masterNode) {
502     os << " </SiteData>\n";
503     }
504     }
505    
506     if (worldRank == masterNode) {
507     os << " </Snapshot>\n";
508     os.flush();
509     }
510    
511 tim 1024 #endif // is_mpi
512 gezelter 1798
513 tim 1024 }
514 gezelter 2
515 gezelter 1769 std::string DumpWriter::prepareDumpLine(StuntDouble* sd) {
516 tim 1024
517 gezelter 1769 int index = sd->getGlobalIntegrableObjectIndex();
518 tim 1024 std::string type("pv");
519     std::string line;
520     char tempBuffer[4096];
521 gezelter 2
522 tim 1024 Vector3d pos;
523     Vector3d vel;
524 gezelter 1769 pos = sd->getPos();
525 gezelter 1437
526     if (isinf(pos[0]) || isnan(pos[0]) ||
527     isinf(pos[1]) || isnan(pos[1]) ||
528     isinf(pos[2]) || isnan(pos[2]) ) {
529     sprintf( painCave.errMsg,
530     "DumpWriter detected a numerical error writing the position"
531     " for object %d", index);
532     painCave.isFatal = 1;
533     simError();
534     }
535    
536 gezelter 1769 vel = sd->getVel();
537 gezelter 1437
538     if (isinf(vel[0]) || isnan(vel[0]) ||
539     isinf(vel[1]) || isnan(vel[1]) ||
540     isinf(vel[2]) || isnan(vel[2]) ) {
541     sprintf( painCave.errMsg,
542     "DumpWriter detected a numerical error writing the velocity"
543     " for object %d", index);
544     painCave.isFatal = 1;
545     simError();
546     }
547    
548 gezelter 1025 sprintf(tempBuffer, "%18.10g %18.10g %18.10g %13e %13e %13e",
549 tim 1024 pos[0], pos[1], pos[2],
550     vel[0], vel[1], vel[2]);
551     line += tempBuffer;
552 gezelter 2
553 gezelter 1769 if (sd->isDirectional()) {
554 tim 1024 type += "qj";
555     Quat4d q;
556     Vector3d ji;
557 gezelter 1769 q = sd->getQ();
558 gezelter 1437
559     if (isinf(q[0]) || isnan(q[0]) ||
560     isinf(q[1]) || isnan(q[1]) ||
561     isinf(q[2]) || isnan(q[2]) ||
562     isinf(q[3]) || isnan(q[3]) ) {
563     sprintf( painCave.errMsg,
564     "DumpWriter detected a numerical error writing the quaternion"
565     " for object %d", index);
566     painCave.isFatal = 1;
567     simError();
568     }
569    
570 gezelter 1769 ji = sd->getJ();
571 gezelter 1437
572     if (isinf(ji[0]) || isnan(ji[0]) ||
573     isinf(ji[1]) || isnan(ji[1]) ||
574     isinf(ji[2]) || isnan(ji[2]) ) {
575     sprintf( painCave.errMsg,
576     "DumpWriter detected a numerical error writing the angular"
577     " momentum for object %d", index);
578     painCave.isFatal = 1;
579     simError();
580     }
581    
582 gezelter 1025 sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e %13e",
583 tim 1024 q[0], q[1], q[2], q[3],
584     ji[0], ji[1], ji[2]);
585     line += tempBuffer;
586     }
587 gezelter 2
588 tim 1024 if (needForceVector_) {
589 gezelter 1437 type += "f";
590 gezelter 1769 Vector3d frc = sd->getFrc();
591 gezelter 1437 if (isinf(frc[0]) || isnan(frc[0]) ||
592     isinf(frc[1]) || isnan(frc[1]) ||
593     isinf(frc[2]) || isnan(frc[2]) ) {
594     sprintf( painCave.errMsg,
595     "DumpWriter detected a numerical error writing the force"
596     " for object %d", index);
597     painCave.isFatal = 1;
598     simError();
599     }
600     sprintf(tempBuffer, " %13e %13e %13e",
601     frc[0], frc[1], frc[2]);
602 tim 1024 line += tempBuffer;
603 gezelter 1437
604 gezelter 1769 if (sd->isDirectional()) {
605 gezelter 1437 type += "t";
606 gezelter 1769 Vector3d trq = sd->getTrq();
607 gezelter 1437 if (isinf(trq[0]) || isnan(trq[0]) ||
608     isinf(trq[1]) || isnan(trq[1]) ||
609     isinf(trq[2]) || isnan(trq[2]) ) {
610     sprintf( painCave.errMsg,
611     "DumpWriter detected a numerical error writing the torque"
612     " for object %d", index);
613     painCave.isFatal = 1;
614     simError();
615 gezelter 1714 }
616 gezelter 1437 sprintf(tempBuffer, " %13e %13e %13e",
617     trq[0], trq[1], trq[2]);
618     line += tempBuffer;
619     }
620 gezelter 246 }
621 gezelter 1714
622 gezelter 1752 sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
623     return std::string(tempBuffer);
624     }
625    
626 gezelter 1769 std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
627 gezelter 1752
628    
629     std::string id;
630     std::string type;
631     std::string line;
632     char tempBuffer[4096];
633    
634 gezelter 1769 if (sd->isRigidBody()) {
635 gezelter 1752 sprintf(tempBuffer, "%10d ", ioIndex);
636     id = std::string(tempBuffer);
637     } else {
638     sprintf(tempBuffer, "%10d %10d", ioIndex, siteIndex);
639     id = std::string(tempBuffer);
640 gezelter 1629 }
641 gezelter 1752
642 gezelter 1714 if (needFlucQ_) {
643     type += "cw";
644 gezelter 1769 RealType fqPos = sd->getFlucQPos();
645 gezelter 1714 if (isinf(fqPos) || isnan(fqPos) ) {
646     sprintf( painCave.errMsg,
647     "DumpWriter detected a numerical error writing the"
648 gezelter 1752 " fluctuating charge for object %s", id.c_str());
649 gezelter 1714 painCave.isFatal = 1;
650     simError();
651     }
652     sprintf(tempBuffer, " %13e ", fqPos);
653     line += tempBuffer;
654    
655 gezelter 1769 RealType fqVel = sd->getFlucQVel();
656 gezelter 1714 if (isinf(fqVel) || isnan(fqVel) ) {
657     sprintf( painCave.errMsg,
658     "DumpWriter detected a numerical error writing the"
659 gezelter 1752 " fluctuating charge velocity for object %s", id.c_str());
660 gezelter 1714 painCave.isFatal = 1;
661     simError();
662     }
663     sprintf(tempBuffer, " %13e ", fqVel);
664     line += tempBuffer;
665    
666     if (needForceVector_) {
667     type += "g";
668 gezelter 1769 RealType fqFrc = sd->getFlucQFrc();
669 gezelter 1714 if (isinf(fqFrc) || isnan(fqFrc) ) {
670     sprintf( painCave.errMsg,
671     "DumpWriter detected a numerical error writing the"
672 gezelter 1752 " fluctuating charge force for object %s", id.c_str());
673 gezelter 1714 painCave.isFatal = 1;
674     simError();
675     }
676     sprintf(tempBuffer, " %13e ", fqFrc);
677     line += tempBuffer;
678     }
679     }
680    
681     if (needElectricField_) {
682     type += "e";
683 gezelter 1769 Vector3d eField= sd->getElectricField();
684 gezelter 1714 if (isinf(eField[0]) || isnan(eField[0]) ||
685     isinf(eField[1]) || isnan(eField[1]) ||
686     isinf(eField[2]) || isnan(eField[2]) ) {
687     sprintf( painCave.errMsg,
688     "DumpWriter detected a numerical error writing the electric"
689 gezelter 1752 " field for object %s", id.c_str());
690 gezelter 1714 painCave.isFatal = 1;
691     simError();
692     }
693     sprintf(tempBuffer, " %13e %13e %13e",
694     eField[0], eField[1], eField[2]);
695     line += tempBuffer;
696     }
697    
698 gezelter 1752
699     if (needParticlePot_) {
700     type += "u";
701 gezelter 1769 RealType particlePot = sd->getParticlePot();
702 gezelter 1752 if (isinf(particlePot) || isnan(particlePot)) {
703     sprintf( painCave.errMsg,
704     "DumpWriter detected a numerical error writing the particle "
705     " potential for object %s", id.c_str());
706     painCave.isFatal = 1;
707     simError();
708     }
709     sprintf(tempBuffer, " %13e", particlePot);
710     line += tempBuffer;
711     }
712    
713    
714     sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
715 tim 1024 return std::string(tempBuffer);
716 gezelter 507 }
717 gezelter 2
718 gezelter 507 void DumpWriter::writeDump() {
719 tim 615 writeFrame(*dumpFile_);
720 gezelter 507 }
721 tim 376
722 gezelter 507 void DumpWriter::writeEor() {
723 tim 615 std::ostream* eorStream;
724 tim 376
725     #ifdef IS_MPI
726     if (worldRank == 0) {
727     #endif // is_mpi
728 gezelter 1874
729 tim 615 eorStream = createOStream(eorFilename_);
730 gezelter 1874 writeFrame(*eorStream);
731    
732 tim 376 #ifdef IS_MPI
733     }
734 tim 615 if (worldRank == 0) {
735     #endif // is_mpi
736 gezelter 1874
737 tim 1024 writeClosing(*eorStream);
738     delete eorStream;
739 gezelter 1874
740 tim 615 #ifdef IS_MPI
741     }
742     #endif // is_mpi
743    
744 gezelter 507 }
745 tim 376
746    
747 gezelter 507 void DumpWriter::writeDumpAndEor() {
748 tim 376 std::vector<std::streambuf*> buffers;
749 tim 615 std::ostream* eorStream;
750 tim 376 #ifdef IS_MPI
751     if (worldRank == 0) {
752     #endif // is_mpi
753    
754 tim 615 buffers.push_back(dumpFile_->rdbuf());
755 tim 376
756 tim 615 eorStream = createOStream(eorFilename_);
757 tim 376
758 tim 615 buffers.push_back(eorStream->rdbuf());
759 tim 376
760     #ifdef IS_MPI
761     }
762     #endif // is_mpi
763    
764     TeeBuf tbuf(buffers.begin(), buffers.end());
765     std::ostream os(&tbuf);
766    
767     writeFrame(os);
768 tim 615
769     #ifdef IS_MPI
770     if (worldRank == 0) {
771     #endif // is_mpi
772 tim 1024 writeClosing(*eorStream);
773     delete eorStream;
774 tim 615 #ifdef IS_MPI
775     }
776     #endif // is_mpi
777 tim 376
778 gezelter 507 }
779 tim 376
780 tim 1024 std::ostream* DumpWriter::createOStream(const std::string& filename) {
781 tim 619
782 tim 615 std::ostream* newOStream;
783 gezelter 1767 #ifdef HAVE_ZLIB
784 tim 615 if (needCompression_) {
785 tim 1024 newOStream = new ogzstream(filename.c_str());
786 tim 615 } else {
787 tim 1024 newOStream = new std::ofstream(filename.c_str());
788 tim 615 }
789 tim 619 #else
790     newOStream = new std::ofstream(filename.c_str());
791     #endif
792 tim 1024 //write out MetaData first
793 gezelter 1746 (*newOStream) << "<OpenMD version=2>" << std::endl;
794 tim 1024 (*newOStream) << " <MetaData>" << std::endl;
795     (*newOStream) << info_->getRawMetaData();
796     (*newOStream) << " </MetaData>" << std::endl;
797 tim 615 return newOStream;
798 tim 1024 }
799 tim 376
800 tim 1024 void DumpWriter::writeClosing(std::ostream& os) {
801    
802 gezelter 1390 os << "</OpenMD>\n";
803 tim 1024 os.flush();
804     }
805    
806 gezelter 1390 }//end namespace OpenMD

Properties

Name Value
svn:keywords Author Id Revision Date