ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/io/DumpWriter.cpp
Revision: 1714
Committed: Sat May 19 18:12:46 2012 UTC (12 years, 11 months ago) by gezelter
File size: 18773 byte(s)
Log Message:
Added ability to read / write dump files with fluctuating charges and
electric fields.

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

Properties

Name Value
svn:keywords Author Id Revision Date