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 |