ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/io/DumpWriter.cpp
Revision: 1711
Committed: Sat May 19 02:58:35 2012 UTC (12 years, 11 months ago) by gezelter
File size: 16419 byte(s)
Log Message:
Some fixes for DataStorage issues.  Removed outdated zangle stuff that
has been replaced by the more modern restraints.

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 cerr << "DW npp = " << needParticlePot_ << "\n";
66 createDumpFile_ = true;
67 #ifdef HAVE_LIBZ
68 if (needCompression_) {
69 filename_ += ".gz";
70 eorFilename_ += ".gz";
71 }
72 #endif
73
74 #ifdef IS_MPI
75
76 if (worldRank == 0) {
77 #endif // is_mpi
78
79 dumpFile_ = createOStream(filename_);
80
81 if (!dumpFile_) {
82 sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
83 filename_.c_str());
84 painCave.isFatal = 1;
85 simError();
86 }
87
88 #ifdef IS_MPI
89
90 }
91
92 #endif // is_mpi
93
94 }
95
96
97 DumpWriter::DumpWriter(SimInfo* info, const std::string& filename)
98 : info_(info), filename_(filename){
99
100 Globals* simParams = info->getSimParams();
101 eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
102
103 needCompression_ = simParams->getCompressDumpFile();
104 needForceVector_ = simParams->getOutputForceVector();
105 createDumpFile_ = true;
106 #ifdef HAVE_LIBZ
107 if (needCompression_) {
108 filename_ += ".gz";
109 eorFilename_ += ".gz";
110 }
111 #endif
112
113 #ifdef IS_MPI
114
115 if (worldRank == 0) {
116 #endif // is_mpi
117
118
119 dumpFile_ = createOStream(filename_);
120
121 if (!dumpFile_) {
122 sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
123 filename_.c_str());
124 painCave.isFatal = 1;
125 simError();
126 }
127
128 #ifdef IS_MPI
129
130 }
131
132 #endif // is_mpi
133
134 }
135
136 DumpWriter::DumpWriter(SimInfo* info, const std::string& filename, bool writeDumpFile)
137 : info_(info), filename_(filename){
138
139 Globals* simParams = info->getSimParams();
140 eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
141
142 needCompression_ = simParams->getCompressDumpFile();
143 needForceVector_ = simParams->getOutputForceVector();
144 needParticlePot_ = simParams->getOutputParticlePotential();
145
146 #ifdef HAVE_LIBZ
147 if (needCompression_) {
148 filename_ += ".gz";
149 eorFilename_ += ".gz";
150 }
151 #endif
152
153 #ifdef IS_MPI
154
155 if (worldRank == 0) {
156 #endif // is_mpi
157
158 createDumpFile_ = writeDumpFile;
159 if (createDumpFile_) {
160 dumpFile_ = createOStream(filename_);
161
162 if (!dumpFile_) {
163 sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
164 filename_.c_str());
165 painCave.isFatal = 1;
166 simError();
167 }
168 }
169 #ifdef IS_MPI
170
171 }
172
173
174 #endif // is_mpi
175
176 }
177
178 DumpWriter::~DumpWriter() {
179
180 #ifdef IS_MPI
181
182 if (worldRank == 0) {
183 #endif // is_mpi
184 if (createDumpFile_){
185 writeClosing(*dumpFile_);
186 delete dumpFile_;
187 }
188 #ifdef IS_MPI
189
190 }
191
192 #endif // is_mpi
193
194 }
195
196 void DumpWriter::writeFrameProperties(std::ostream& os, Snapshot* s) {
197
198 char buffer[1024];
199
200 os << " <FrameData>\n";
201
202 RealType currentTime = s->getTime();
203
204 if (isinf(currentTime) || isnan(currentTime)) {
205 sprintf( painCave.errMsg,
206 "DumpWriter detected a numerical error writing the time");
207 painCave.isFatal = 1;
208 simError();
209 }
210
211 sprintf(buffer, " Time: %.10g\n", currentTime);
212 os << buffer;
213
214 Mat3x3d hmat;
215 hmat = s->getHmat();
216
217 for (unsigned int i = 0; i < 3; i++) {
218 for (unsigned int j = 0; j < 3; j++) {
219 if (isinf(hmat(i,j)) || isnan(hmat(i,j))) {
220 sprintf( painCave.errMsg,
221 "DumpWriter detected a numerical error writing the box");
222 painCave.isFatal = 1;
223 simError();
224 }
225 }
226 }
227
228 sprintf(buffer, " Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
229 hmat(0, 0), hmat(1, 0), hmat(2, 0),
230 hmat(0, 1), hmat(1, 1), hmat(2, 1),
231 hmat(0, 2), hmat(1, 2), hmat(2, 2));
232 os << buffer;
233
234 RealType chi = s->getChi();
235 RealType integralOfChiDt = s->getIntegralOfChiDt();
236 if (isinf(chi) || isnan(chi) ||
237 isinf(integralOfChiDt) || isnan(integralOfChiDt)) {
238 sprintf( painCave.errMsg,
239 "DumpWriter detected a numerical error writing the thermostat");
240 painCave.isFatal = 1;
241 simError();
242 }
243 sprintf(buffer, " Thermostat: %.10g , %.10g\n", chi, integralOfChiDt);
244 os << buffer;
245
246 Mat3x3d eta;
247 eta = s->getEta();
248
249 for (unsigned int i = 0; i < 3; i++) {
250 for (unsigned int j = 0; j < 3; j++) {
251 if (isinf(eta(i,j)) || isnan(eta(i,j))) {
252 sprintf( painCave.errMsg,
253 "DumpWriter detected a numerical error writing the barostat");
254 painCave.isFatal = 1;
255 simError();
256 }
257 }
258 }
259
260 sprintf(buffer, " Barostat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
261 eta(0, 0), eta(1, 0), eta(2, 0),
262 eta(0, 1), eta(1, 1), eta(2, 1),
263 eta(0, 2), eta(1, 2), eta(2, 2));
264 os << buffer;
265
266 os << " </FrameData>\n";
267 }
268
269 void DumpWriter::writeFrame(std::ostream& os) {
270
271 #ifdef IS_MPI
272 MPI_Status istatus;
273 #endif
274
275 Molecule* mol;
276 StuntDouble* integrableObject;
277 SimInfo::MoleculeIterator mi;
278 Molecule::IntegrableObjectIterator ii;
279
280 #ifndef IS_MPI
281 os << " <Snapshot>\n";
282
283 writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
284
285 os << " <StuntDoubles>\n";
286 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
287
288
289 for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
290 integrableObject = mol->nextIntegrableObject(ii)) {
291 os << prepareDumpLine(integrableObject);
292
293 }
294 }
295 os << " </StuntDoubles>\n";
296
297 os << " </Snapshot>\n";
298
299 os.flush();
300 #else
301 //every node prepares the dump lines for integrable objects belong to itself
302 std::string buffer;
303 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
304
305
306 for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
307 integrableObject = mol->nextIntegrableObject(ii)) {
308 buffer += prepareDumpLine(integrableObject);
309 }
310 }
311
312 const int masterNode = 0;
313 int nProc;
314 MPI_Comm_size(MPI_COMM_WORLD, &nProc);
315 if (worldRank == masterNode) {
316 os << " <Snapshot>\n";
317 writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
318 os << " <StuntDoubles>\n";
319
320 os << buffer;
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