ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/io/DumpWriter.cpp
Revision: 1798
Committed: Thu Sep 13 14:10:11 2012 UTC (12 years, 7 months ago) by gezelter
File size: 23835 byte(s)
Log Message:
Merged trunk changes into the development branch

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     * [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 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     Atom* atom;
314 gezelter 2
315 gezelter 246 #ifndef IS_MPI
316 tim 1024 os << " <Snapshot>\n";
317 gezelter 1025
318 tim 1024 writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
319 gezelter 2
320 tim 1024 os << " <StuntDoubles>\n";
321 gezelter 246 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
322 gezelter 2
323 xsun 1217
324 gezelter 1769 for (sd = mol->beginIntegrableObject(ii); sd != NULL;
325     sd = mol->nextIntegrableObject(ii)) {
326     os << prepareDumpLine(sd);
327 xsun 1217
328 tim 1024 }
329     }
330     os << " </StuntDoubles>\n";
331 gezelter 1752
332     if (doSiteData_) {
333     os << " <SiteData>\n";
334 gezelter 1798 for (mol = info_->beginMolecule(mi); mol != NULL;
335     mol = info_->nextMolecule(mi)) {
336 gezelter 1752
337 gezelter 1769 for (sd = mol->beginIntegrableObject(ii); sd != NULL;
338     sd = mol->nextIntegrableObject(ii)) {
339 gezelter 1752
340 gezelter 1769 int ioIndex = sd->getGlobalIntegrableObjectIndex();
341 gezelter 1752 // do one for the IO itself
342 gezelter 1769 os << prepareSiteLine(sd, ioIndex, 0);
343 gezelter 1752
344 gezelter 1769 if (sd->isRigidBody()) {
345 gezelter 1752
346 gezelter 1769 RigidBody* rb = static_cast<RigidBody*>(sd);
347 gezelter 1752 int siteIndex = 0;
348     for (atom = rb->beginAtom(ai); atom != NULL;
349     atom = rb->nextAtom(ai)) {
350     os << prepareSiteLine(atom, ioIndex, siteIndex);
351     siteIndex++;
352     }
353     }
354     }
355     }
356     os << " </SiteData>\n";
357     }
358 tim 1024 os << " </Snapshot>\n";
359 gezelter 2
360 tim 1024 os.flush();
361     #else
362 gezelter 1798
363     const int masterNode = 0;
364     int worldRank = MPI::COMM_WORLD.Get_rank();
365     int nProc = MPI::COMM_WORLD.Get_size();
366    
367     if (worldRank == masterNode) {
368     os << " <Snapshot>\n";
369     writeFrameProperties(os,
370     info_->getSnapshotManager()->getCurrentSnapshot());
371     os << " <StuntDoubles>\n";
372     }
373    
374 tim 1024 //every node prepares the dump lines for integrable objects belong to itself
375     std::string buffer;
376 gezelter 1798 for (mol = info_->beginMolecule(mi); mol != NULL;
377     mol = info_->nextMolecule(mi)) {
378 gezelter 1769 for (sd = mol->beginIntegrableObject(ii); sd != NULL;
379     sd = mol->nextIntegrableObject(ii)) {
380 gezelter 1798 buffer += prepareDumpLine(sd);
381 gezelter 507 }
382 gezelter 2 }
383 tim 1024
384     if (worldRank == masterNode) {
385 gezelter 1025 os << buffer;
386 gezelter 1798
387 tim 1024 for (int i = 1; i < nProc; ++i) {
388 gezelter 1798 // tell processor i to start sending us data:
389     MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode);
390 gezelter 2
391 tim 1024 // receive the length of the string buffer that was
392 gezelter 1798 // prepared by processor i:
393     int recvLength;
394     MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
395     istatus);
396 gezelter 2
397 gezelter 1798 // create a buffer to receive the data
398 tim 1024 char* recvBuffer = new char[recvLength];
399     if (recvBuffer == NULL) {
400     } else {
401 gezelter 1798 // receive the data:
402     MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i,
403     MPI::ANY_TAG, istatus);
404     // send it to the file:
405 tim 1024 os << recvBuffer;
406 gezelter 1798 // get rid of the receive buffer:
407 gezelter 1313 delete [] recvBuffer;
408 tim 1024 }
409     }
410 gezelter 246 } else {
411 gezelter 1025 int sendBufferLength = buffer.size() + 1;
412 gezelter 1629 int myturn = 0;
413     for (int i = 1; i < nProc; ++i){
414 gezelter 1798 // wait for the master node to call our number:
415     MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
416 gezelter 1629 if (myturn == worldRank){
417 gezelter 1798 // send the length of our buffer:
418     MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0);
419    
420     // send our buffer:
421     MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength,
422     MPI::CHAR, masterNode, 0);
423 gezelter 1629 }
424     }
425 tim 1024 }
426 gezelter 1798
427     if (worldRank == masterNode) {
428     os << " </StuntDoubles>\n";
429     }
430 gezelter 2
431 gezelter 1798 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     for (atom = rb->beginAtom(ai); atom != NULL;
451     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 1798
514 tim 1024 }
515 gezelter 2
516 gezelter 1769 std::string DumpWriter::prepareDumpLine(StuntDouble* sd) {
517 tim 1024
518 gezelter 1769 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 1769 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 1769 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 1769 if (sd->isDirectional()) {
555 tim 1024 type += "qj";
556     Quat4d q;
557     Vector3d ji;
558 gezelter 1769 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 1769 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 1769 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 1769 if (sd->isDirectional()) {
606 gezelter 1437 type += "t";
607 gezelter 1769 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 1714 }
617 gezelter 1437 sprintf(tempBuffer, " %13e %13e %13e",
618     trq[0], trq[1], trq[2]);
619     line += tempBuffer;
620     }
621 gezelter 246 }
622 gezelter 1714
623 gezelter 1752 sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
624     return std::string(tempBuffer);
625     }
626    
627 gezelter 1769 std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
628 gezelter 1752
629    
630     std::string id;
631     std::string type;
632     std::string line;
633     char tempBuffer[4096];
634    
635 gezelter 1769 if (sd->isRigidBody()) {
636 gezelter 1752 sprintf(tempBuffer, "%10d ", ioIndex);
637     id = std::string(tempBuffer);
638     } else {
639     sprintf(tempBuffer, "%10d %10d", ioIndex, siteIndex);
640     id = std::string(tempBuffer);
641 gezelter 1629 }
642 gezelter 1752
643 gezelter 1714 if (needFlucQ_) {
644     type += "cw";
645 gezelter 1769 RealType fqPos = sd->getFlucQPos();
646 gezelter 1714 if (isinf(fqPos) || isnan(fqPos) ) {
647     sprintf( painCave.errMsg,
648     "DumpWriter detected a numerical error writing the"
649 gezelter 1752 " fluctuating charge for object %s", id.c_str());
650 gezelter 1714 painCave.isFatal = 1;
651     simError();
652     }
653     sprintf(tempBuffer, " %13e ", fqPos);
654     line += tempBuffer;
655    
656 gezelter 1769 RealType fqVel = sd->getFlucQVel();
657 gezelter 1714 if (isinf(fqVel) || isnan(fqVel) ) {
658     sprintf( painCave.errMsg,
659     "DumpWriter detected a numerical error writing the"
660 gezelter 1752 " fluctuating charge velocity for object %s", id.c_str());
661 gezelter 1714 painCave.isFatal = 1;
662     simError();
663     }
664     sprintf(tempBuffer, " %13e ", fqVel);
665     line += tempBuffer;
666    
667     if (needForceVector_) {
668     type += "g";
669 gezelter 1769 RealType fqFrc = sd->getFlucQFrc();
670 gezelter 1714 if (isinf(fqFrc) || isnan(fqFrc) ) {
671     sprintf( painCave.errMsg,
672     "DumpWriter detected a numerical error writing the"
673 gezelter 1752 " fluctuating charge force for object %s", id.c_str());
674 gezelter 1714 painCave.isFatal = 1;
675     simError();
676     }
677     sprintf(tempBuffer, " %13e ", fqFrc);
678     line += tempBuffer;
679     }
680     }
681    
682     if (needElectricField_) {
683     type += "e";
684 gezelter 1769 Vector3d eField= sd->getElectricField();
685 gezelter 1714 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 gezelter 1752 " field for object %s", id.c_str());
691 gezelter 1714 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 gezelter 1752
700     if (needParticlePot_) {
701     type += "u";
702 gezelter 1769 RealType particlePot = sd->getParticlePot();
703 gezelter 1752 if (isinf(particlePot) || isnan(particlePot)) {
704     sprintf( painCave.errMsg,
705     "DumpWriter detected a numerical error writing the particle "
706     " potential for object %s", id.c_str());
707     painCave.isFatal = 1;
708     simError();
709     }
710     sprintf(tempBuffer, " %13e", particlePot);
711     line += tempBuffer;
712     }
713    
714    
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    
730 tim 615 eorStream = createOStream(eorFilename_);
731 tim 376
732     #ifdef IS_MPI
733     }
734     #endif // is_mpi
735    
736 tim 615 writeFrame(*eorStream);
737    
738     #ifdef IS_MPI
739     if (worldRank == 0) {
740     #endif // is_mpi
741 tim 1024 writeClosing(*eorStream);
742     delete eorStream;
743 tim 615 #ifdef IS_MPI
744     }
745     #endif // is_mpi
746    
747 gezelter 507 }
748 tim 376
749    
750 gezelter 507 void DumpWriter::writeDumpAndEor() {
751 tim 376 std::vector<std::streambuf*> buffers;
752 tim 615 std::ostream* eorStream;
753 tim 376 #ifdef IS_MPI
754     if (worldRank == 0) {
755     #endif // is_mpi
756    
757 tim 615 buffers.push_back(dumpFile_->rdbuf());
758 tim 376
759 tim 615 eorStream = createOStream(eorFilename_);
760 tim 376
761 tim 615 buffers.push_back(eorStream->rdbuf());
762 tim 376
763     #ifdef IS_MPI
764     }
765     #endif // is_mpi
766    
767     TeeBuf tbuf(buffers.begin(), buffers.end());
768     std::ostream os(&tbuf);
769    
770     writeFrame(os);
771 tim 615
772     #ifdef IS_MPI
773     if (worldRank == 0) {
774     #endif // is_mpi
775 tim 1024 writeClosing(*eorStream);
776     delete eorStream;
777 tim 615 #ifdef IS_MPI
778     }
779     #endif // is_mpi
780 tim 376
781 gezelter 507 }
782 tim 376
783 tim 1024 std::ostream* DumpWriter::createOStream(const std::string& filename) {
784 tim 619
785 tim 615 std::ostream* newOStream;
786 gezelter 1767 #ifdef HAVE_ZLIB
787 tim 615 if (needCompression_) {
788 tim 1024 newOStream = new ogzstream(filename.c_str());
789 tim 615 } else {
790 tim 1024 newOStream = new std::ofstream(filename.c_str());
791 tim 615 }
792 tim 619 #else
793     newOStream = new std::ofstream(filename.c_str());
794     #endif
795 tim 1024 //write out MetaData first
796 gezelter 1746 (*newOStream) << "<OpenMD version=2>" << std::endl;
797 tim 1024 (*newOStream) << " <MetaData>" << std::endl;
798     (*newOStream) << info_->getRawMetaData();
799     (*newOStream) << " </MetaData>" << std::endl;
800 tim 615 return newOStream;
801 tim 1024 }
802 tim 376
803 tim 1024 void DumpWriter::writeClosing(std::ostream& os) {
804    
805 gezelter 1390 os << "</OpenMD>\n";
806 tim 1024 os.flush();
807     }
808    
809 gezelter 1390 }//end namespace OpenMD

Properties

Name Value
svn:keywords Author Id Revision Date