ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/io/DumpWriter.cpp
Revision: 2020
Committed: Mon Sep 22 19:18:35 2014 UTC (10 years, 7 months ago) by gezelter
File size: 25280 byte(s)
Log Message:
Fixes for restraints, renaming of UniformField

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

Properties

Name Value
svn:keywords Author Id Revision Date