ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/io/DumpWriter.cpp
Revision: 1665
Committed: Tue Nov 22 20:38:56 2011 UTC (13 years, 5 months ago) by gezelter
File size: 16366 byte(s)
Log Message:
updated copyright notices

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] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 */
42
43 #include "io/DumpWriter.hpp"
44 #include "primitives/Molecule.hpp"
45 #include "utils/simError.h"
46 #include "io/basic_teebuf.hpp"
47 #include "io/gzstream.hpp"
48 #include "io/Globals.hpp"
49
50
51 #ifdef IS_MPI
52 #include <mpi.h>
53 #endif //is_mpi
54
55 using namespace std;
56 namespace OpenMD {
57
58 DumpWriter::DumpWriter(SimInfo* info)
59 : info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
60
61 Globals* simParams = info->getSimParams();
62 needCompression_ = simParams->getCompressDumpFile();
63 needForceVector_ = simParams->getOutputForceVector();
64 needParticlePot_ = simParams->getOutputParticlePotential();
65 createDumpFile_ = true;
66 #ifdef HAVE_LIBZ
67 if (needCompression_) {
68 filename_ += ".gz";
69 eorFilename_ += ".gz";
70 }
71 #endif
72
73 #ifdef IS_MPI
74
75 if (worldRank == 0) {
76 #endif // is_mpi
77
78 dumpFile_ = createOStream(filename_);
79
80 if (!dumpFile_) {
81 sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
82 filename_.c_str());
83 painCave.isFatal = 1;
84 simError();
85 }
86
87 #ifdef IS_MPI
88
89 }
90
91 #endif // is_mpi
92
93 }
94
95
96 DumpWriter::DumpWriter(SimInfo* info, const std::string& filename)
97 : info_(info), filename_(filename){
98
99 Globals* simParams = info->getSimParams();
100 eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
101
102 needCompression_ = simParams->getCompressDumpFile();
103 needForceVector_ = simParams->getOutputForceVector();
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 for (int i = 1; i < nProc; ++i) {
322
323 // receive the length of the string buffer that was
324 // prepared by processor i
325
326 MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
327 int recvLength;
328 MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
329 char* recvBuffer = new char[recvLength];
330 if (recvBuffer == NULL) {
331 } else {
332 MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
333 os << recvBuffer;
334 delete [] recvBuffer;
335 }
336 }
337 os << " </StuntDoubles>\n";
338
339 os << " </Snapshot>\n";
340 os.flush();
341 } else {
342 int sendBufferLength = buffer.size() + 1;
343 int myturn = 0;
344 for (int i = 1; i < nProc; ++i){
345 MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
346 if (myturn == worldRank){
347 MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
348 MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
349 }
350 }
351 }
352
353 #endif // is_mpi
354
355 }
356
357 std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
358
359 int index = integrableObject->getGlobalIntegrableObjectIndex();
360 std::string type("pv");
361 std::string line;
362 char tempBuffer[4096];
363
364 Vector3d pos;
365 Vector3d vel;
366 pos = integrableObject->getPos();
367
368 if (isinf(pos[0]) || isnan(pos[0]) ||
369 isinf(pos[1]) || isnan(pos[1]) ||
370 isinf(pos[2]) || isnan(pos[2]) ) {
371 sprintf( painCave.errMsg,
372 "DumpWriter detected a numerical error writing the position"
373 " for object %d", index);
374 painCave.isFatal = 1;
375 simError();
376 }
377
378 vel = integrableObject->getVel();
379
380 if (isinf(vel[0]) || isnan(vel[0]) ||
381 isinf(vel[1]) || isnan(vel[1]) ||
382 isinf(vel[2]) || isnan(vel[2]) ) {
383 sprintf( painCave.errMsg,
384 "DumpWriter detected a numerical error writing the velocity"
385 " for object %d", index);
386 painCave.isFatal = 1;
387 simError();
388 }
389
390 sprintf(tempBuffer, "%18.10g %18.10g %18.10g %13e %13e %13e",
391 pos[0], pos[1], pos[2],
392 vel[0], vel[1], vel[2]);
393 line += tempBuffer;
394
395 if (integrableObject->isDirectional()) {
396 type += "qj";
397 Quat4d q;
398 Vector3d ji;
399 q = integrableObject->getQ();
400
401 if (isinf(q[0]) || isnan(q[0]) ||
402 isinf(q[1]) || isnan(q[1]) ||
403 isinf(q[2]) || isnan(q[2]) ||
404 isinf(q[3]) || isnan(q[3]) ) {
405 sprintf( painCave.errMsg,
406 "DumpWriter detected a numerical error writing the quaternion"
407 " for object %d", index);
408 painCave.isFatal = 1;
409 simError();
410 }
411
412 ji = integrableObject->getJ();
413
414 if (isinf(ji[0]) || isnan(ji[0]) ||
415 isinf(ji[1]) || isnan(ji[1]) ||
416 isinf(ji[2]) || isnan(ji[2]) ) {
417 sprintf( painCave.errMsg,
418 "DumpWriter detected a numerical error writing the angular"
419 " momentum for object %d", index);
420 painCave.isFatal = 1;
421 simError();
422 }
423
424 sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e %13e",
425 q[0], q[1], q[2], q[3],
426 ji[0], ji[1], ji[2]);
427 line += tempBuffer;
428 }
429
430 if (needForceVector_) {
431 type += "f";
432 Vector3d frc;
433
434 frc = integrableObject->getFrc();
435
436 if (isinf(frc[0]) || isnan(frc[0]) ||
437 isinf(frc[1]) || isnan(frc[1]) ||
438 isinf(frc[2]) || isnan(frc[2]) ) {
439 sprintf( painCave.errMsg,
440 "DumpWriter detected a numerical error writing the force"
441 " for object %d", index);
442 painCave.isFatal = 1;
443 simError();
444 }
445 sprintf(tempBuffer, " %13e %13e %13e",
446 frc[0], frc[1], frc[2]);
447 line += tempBuffer;
448
449 if (integrableObject->isDirectional()) {
450 type += "t";
451 Vector3d trq;
452
453 trq = integrableObject->getTrq();
454
455 if (isinf(trq[0]) || isnan(trq[0]) ||
456 isinf(trq[1]) || isnan(trq[1]) ||
457 isinf(trq[2]) || isnan(trq[2]) ) {
458 sprintf( painCave.errMsg,
459 "DumpWriter detected a numerical error writing the torque"
460 " for object %d", index);
461 painCave.isFatal = 1;
462 simError();
463 }
464
465 sprintf(tempBuffer, " %13e %13e %13e",
466 trq[0], trq[1], trq[2]);
467 line += tempBuffer;
468 }
469 }
470 if (needParticlePot_) {
471 type += "u";
472 RealType particlePot;
473
474 particlePot = integrableObject->getParticlePot();
475
476 if (isinf(particlePot) || isnan(particlePot)) {
477 sprintf( painCave.errMsg,
478 "DumpWriter detected a numerical error writing the particle "
479 " potential for object %d", index);
480 painCave.isFatal = 1;
481 simError();
482 }
483 sprintf(tempBuffer, " %13e", particlePot);
484 line += tempBuffer;
485 }
486
487 sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
488 return std::string(tempBuffer);
489 }
490
491 void DumpWriter::writeDump() {
492 writeFrame(*dumpFile_);
493 }
494
495 void DumpWriter::writeEor() {
496 std::ostream* eorStream;
497
498 #ifdef IS_MPI
499 if (worldRank == 0) {
500 #endif // is_mpi
501
502 eorStream = createOStream(eorFilename_);
503
504 #ifdef IS_MPI
505 }
506 #endif // is_mpi
507
508 writeFrame(*eorStream);
509
510 #ifdef IS_MPI
511 if (worldRank == 0) {
512 #endif // is_mpi
513 writeClosing(*eorStream);
514 delete eorStream;
515 #ifdef IS_MPI
516 }
517 #endif // is_mpi
518
519 }
520
521
522 void DumpWriter::writeDumpAndEor() {
523 std::vector<std::streambuf*> buffers;
524 std::ostream* eorStream;
525 #ifdef IS_MPI
526 if (worldRank == 0) {
527 #endif // is_mpi
528
529 buffers.push_back(dumpFile_->rdbuf());
530
531 eorStream = createOStream(eorFilename_);
532
533 buffers.push_back(eorStream->rdbuf());
534
535 #ifdef IS_MPI
536 }
537 #endif // is_mpi
538
539 TeeBuf tbuf(buffers.begin(), buffers.end());
540 std::ostream os(&tbuf);
541
542 writeFrame(os);
543
544 #ifdef IS_MPI
545 if (worldRank == 0) {
546 #endif // is_mpi
547 writeClosing(*eorStream);
548 delete eorStream;
549 #ifdef IS_MPI
550 }
551 #endif // is_mpi
552
553 }
554
555 std::ostream* DumpWriter::createOStream(const std::string& filename) {
556
557 std::ostream* newOStream;
558 #ifdef HAVE_LIBZ
559 if (needCompression_) {
560 newOStream = new ogzstream(filename.c_str());
561 } else {
562 newOStream = new std::ofstream(filename.c_str());
563 }
564 #else
565 newOStream = new std::ofstream(filename.c_str());
566 #endif
567 //write out MetaData first
568 (*newOStream) << "<OpenMD version=1>" << std::endl;
569 (*newOStream) << " <MetaData>" << std::endl;
570 (*newOStream) << info_->getRawMetaData();
571 (*newOStream) << " </MetaData>" << std::endl;
572 return newOStream;
573 }
574
575 void DumpWriter::writeClosing(std::ostream& os) {
576
577 os << "</OpenMD>\n";
578 os.flush();
579 }
580
581 }//end namespace OpenMD

Properties

Name Value
svn:keywords Author Id Revision Date