ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/io/DumpWriter.cpp
Revision: 1712
Committed: Sat May 19 13:30:21 2012 UTC (12 years, 11 months ago) by gezelter
File size: 16483 byte(s)
Log Message:
Bugfixes (mostly related to particlePot and storageLayout).

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

Properties

Name Value
svn:keywords Author Id Revision Date