ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/io/DumpWriter.cpp
Revision: 1610
Committed: Fri Aug 12 14:37:25 2011 UTC (13 years, 8 months ago) by gezelter
File size: 16377 byte(s)
Log Message:
Fixed a clang compilation problem, added the ability to output
particle potential in the dump files.

File Contents

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

Properties

Name Value
svn:keywords Author Id Revision Date