OpenMD 3.1
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
DumpWriter.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include "io/DumpWriter.hpp"
46
47#ifdef IS_MPI
48#include <mpi.h>
49#endif
50
51#include <config.h>
52
53#include "io/Globals.hpp"
54#include "io/basic_teebuf.hpp"
55#include "primitives/Molecule.hpp"
56#include "utils/simError.h"
57#ifdef HAVE_ZLIB
58#include "io/gzstream.hpp"
59#endif
60
61using namespace std;
62namespace OpenMD {
63
64 DumpWriter::DumpWriter(SimInfo* info) :
65 info_(info), filename_(info->getDumpFileName()),
66 eorFilename_(info->getFinalConfigFileName()) {
67 Globals* simParams = info->getSimParams();
68 needCompression_ = simParams->getCompressDumpFile();
69 needForceVector_ = simParams->getOutputForceVector();
70 needParticlePot_ = simParams->getOutputParticlePotential();
71 needFlucQ_ = simParams->getOutputFluctuatingCharges();
72 needElectricField_ = simParams->getOutputElectricField();
73 needSitePotential_ = simParams->getOutputSitePotential();
74 needDensity_ = simParams->getOutputDensity();
75
76 if (needParticlePot_ || needFlucQ_ || needElectricField_ ||
77 needSitePotential_ || needDensity_) {
78 doSiteData_ = true;
79 } else {
80 doSiteData_ = false;
81 }
82
83 createDumpFile_ = true;
84#ifdef HAVE_LIBZ
85 if (needCompression_) {
86 filename_ += ".gz";
87 eorFilename_ += ".gz";
88 }
89#endif
90
91#ifdef IS_MPI
92
93 if (worldRank == 0) {
94#endif // is_mpi
95
96 dumpFile_ = createOStream(filename_);
97
98 if (!dumpFile_) {
99 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
100 "Could not open \"%s\" for dump output.\n", filename_.c_str());
101 painCave.isFatal = 1;
102 simError();
103 }
104
105#ifdef IS_MPI
106 }
107
108#endif // is_mpi
109 }
110
111 DumpWriter::DumpWriter(SimInfo* info, const std::string& filename) :
112 info_(info), filename_(filename) {
113 Globals* simParams = info->getSimParams();
114 eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
115
116 needCompression_ = simParams->getCompressDumpFile();
117 needForceVector_ = simParams->getOutputForceVector();
118 needParticlePot_ = simParams->getOutputParticlePotential();
119 needFlucQ_ = simParams->getOutputFluctuatingCharges();
120 needElectricField_ = simParams->getOutputElectricField();
121 needSitePotential_ = simParams->getOutputSitePotential();
122 needDensity_ = simParams->getOutputDensity();
123
124 if (needParticlePot_ || needFlucQ_ || needElectricField_ ||
125 needSitePotential_ || needDensity_) {
126 doSiteData_ = true;
127 } else {
128 doSiteData_ = false;
129 }
130
131 createDumpFile_ = true;
132#ifdef HAVE_LIBZ
133 if (needCompression_) {
134 filename_ += ".gz";
135 eorFilename_ += ".gz";
136 }
137#endif
138
139#ifdef IS_MPI
140
141 if (worldRank == 0) {
142#endif // is_mpi
143
144 dumpFile_ = createOStream(filename_);
145
146 if (!dumpFile_) {
147 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
148 "Could not open \"%s\" for dump output.\n", filename_.c_str());
149 painCave.isFatal = 1;
150 simError();
151 }
152
153#ifdef IS_MPI
154 }
155
156#endif // is_mpi
157 }
158
159 DumpWriter::DumpWriter(SimInfo* info, const std::string& filename,
160 bool writeDumpFile) :
161 info_(info),
162 filename_(filename) {
163 Globals* simParams = info->getSimParams();
164 eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
165
166 needCompression_ = simParams->getCompressDumpFile();
167 needForceVector_ = simParams->getOutputForceVector();
168 needParticlePot_ = simParams->getOutputParticlePotential();
169 needFlucQ_ = simParams->getOutputFluctuatingCharges();
170 needElectricField_ = simParams->getOutputElectricField();
171 needSitePotential_ = simParams->getOutputSitePotential();
172 needDensity_ = simParams->getOutputDensity();
173
174 if (needParticlePot_ || needFlucQ_ || needElectricField_ ||
175 needSitePotential_ || needDensity_) {
176 doSiteData_ = true;
177 } else {
178 doSiteData_ = false;
179 }
180
181#ifdef HAVE_LIBZ
182 if (needCompression_) {
183 filename_ += ".gz";
184 eorFilename_ += ".gz";
185 }
186#endif
187
188#ifdef IS_MPI
189
190 if (worldRank == 0) {
191#endif // is_mpi
192
193 createDumpFile_ = writeDumpFile;
194 if (createDumpFile_) {
195 dumpFile_ = createOStream(filename_);
196
197 if (!dumpFile_) {
198 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
199 "Could not open \"%s\" for dump output.\n",
200 filename_.c_str());
201 painCave.isFatal = 1;
202 simError();
203 }
204 }
205#ifdef IS_MPI
206 }
207
208#endif // is_mpi
209 }
210
211 DumpWriter::~DumpWriter() {
212#ifdef IS_MPI
213
214 if (worldRank == 0) {
215#endif // is_mpi
216 if (createDumpFile_) {
217 writeClosing(*dumpFile_);
218 delete dumpFile_;
219 }
220#ifdef IS_MPI
221 }
222
223#endif // is_mpi
224 }
225
226 void DumpWriter::writeFrameProperties(std::ostream& os, Snapshot* s) {
227 char buffer[1024];
228
229 os << " <FrameData>\n";
230
231 RealType currentTime = s->getTime();
232
233 if (std::isinf(currentTime) || std::isnan(currentTime)) {
234 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
235 "DumpWriter detected a numerical error writing the time");
236 painCave.isFatal = 1;
237 simError();
238 }
239
240 snprintf(buffer, 1024, " Time: %.10g\n", currentTime);
241 os << buffer;
242
243 Mat3x3d hmat;
244 hmat = s->getHmat();
245
246 for (unsigned int i = 0; i < 3; i++) {
247 for (unsigned int j = 0; j < 3; j++) {
248 if (std::isinf(hmat(i, j)) || std::isnan(hmat(i, j))) {
249 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
250 "DumpWriter detected a numerical error writing the box");
251 painCave.isFatal = 1;
252 simError();
253 }
254 }
255 }
256
257 snprintf(
258 buffer, 1024,
259 " Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { "
260 "%.10g, "
261 "%.10g, %.10g }}\n",
262 hmat(0, 0), hmat(1, 0), hmat(2, 0), hmat(0, 1), hmat(1, 1), hmat(2, 1),
263 hmat(0, 2), hmat(1, 2), hmat(2, 2));
264 os << buffer;
265
266 pair<RealType, RealType> thermostat = s->getThermostat();
267
268 if (std::isinf(thermostat.first) || std::isnan(thermostat.first) ||
269 std::isinf(thermostat.second) || std::isnan(thermostat.second)) {
270 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
271 "DumpWriter detected a numerical error writing the thermostat");
272 painCave.isFatal = 1;
273 simError();
274 }
275 snprintf(buffer, 1024, " Thermostat: %.10g , %.10g\n", thermostat.first,
276 thermostat.second);
277 os << buffer;
278
279 Mat3x3d eta;
280 eta = s->getBarostat();
281
282 for (unsigned int i = 0; i < 3; i++) {
283 for (unsigned int j = 0; j < 3; j++) {
284 if (std::isinf(eta(i, j)) || std::isnan(eta(i, j))) {
285 snprintf(
286 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
287 "DumpWriter detected a numerical error writing the barostat");
288 painCave.isFatal = 1;
289 simError();
290 }
291 }
292 }
293
294 snprintf(
295 buffer, 1024,
296 " Barostat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { "
297 "%.10g, "
298 "%.10g, %.10g }}\n",
299 eta(0, 0), eta(1, 0), eta(2, 0), eta(0, 1), eta(1, 1), eta(2, 1),
300 eta(0, 2), eta(1, 2), eta(2, 2));
301 os << buffer;
302
303 // SPF Data
304 std::shared_ptr<SPFData> spfData = s->getSPFData();
305
306 if (std::isinf(spfData->pos[0]) || std::isnan(spfData->pos[0]) ||
307 std::isinf(spfData->pos[1]) || std::isnan(spfData->pos[1]) ||
308 std::isinf(spfData->pos[2]) || std::isnan(spfData->pos[2]) ||
309 std::isinf(spfData->lambda) || std::isnan(spfData->lambda) ||
310 std::isinf(spfData->globalID) || std::isnan(spfData->globalID)) {
311 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
312 "DumpWriter detected a numerical error writing the spf data "
313 "structure");
314 painCave.isFatal = 1;
315 simError();
316 }
317 snprintf(buffer, 1024,
318 " SPFData: {{ %.10g, %.10g, %.10g }, %.10g, %d }\n",
319 spfData->pos[0], spfData->pos[1], spfData->pos[2], spfData->lambda,
320 spfData->globalID);
321 os << buffer;
322
323 os << " </FrameData>\n";
324 }
325
326 void DumpWriter::writeFrame(std::ostream& os) {
327#ifdef IS_MPI
328 MPI_Status istatus;
329#endif
330
331 Molecule* mol;
332 StuntDouble* sd;
333 SimInfo::MoleculeIterator mi;
334 Molecule::IntegrableObjectIterator ii;
335 RigidBody::AtomIterator ai;
336
337#ifndef IS_MPI
338 os << " <Snapshot>\n";
339
340 writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
341
342 os << " <StuntDoubles>\n";
343 for (mol = info_->beginMolecule(mi); mol != NULL;
344 mol = info_->nextMolecule(mi)) {
345 for (sd = mol->beginIntegrableObject(ii); sd != NULL;
346 sd = mol->nextIntegrableObject(ii)) {
347 os << prepareDumpLine(sd);
348 }
349 }
350 os << " </StuntDoubles>\n";
351
352 if (doSiteData_) {
353 os << " <SiteData>\n";
354 for (mol = info_->beginMolecule(mi); mol != NULL;
355 mol = info_->nextMolecule(mi)) {
356 for (sd = mol->beginIntegrableObject(ii); sd != NULL;
357 sd = mol->nextIntegrableObject(ii)) {
358 int ioIndex = sd->getGlobalIntegrableObjectIndex();
359 // do one for the IO itself
360 os << prepareSiteLine(sd, ioIndex, 0);
361
362 if (sd->isRigidBody()) {
363 RigidBody* rb = static_cast<RigidBody*>(sd);
364 int siteIndex = 0;
365 for (Atom* atom = rb->beginAtom(ai); atom != NULL;
366 atom = rb->nextAtom(ai)) {
367 os << prepareSiteLine(atom, ioIndex, siteIndex);
368 siteIndex++;
369 }
370 }
371 }
372 }
373 os << " </SiteData>\n";
374 }
375 os << " </Snapshot>\n";
376
377 os.flush();
378 os.rdbuf()->pubsync();
379#else
380
381 const int primaryNode = 0;
382 int worldRank;
383 int nProc;
384
385 MPI_Comm_size(MPI_COMM_WORLD, &nProc);
386 MPI_Comm_rank(MPI_COMM_WORLD, &worldRank);
387
388 if (worldRank == primaryNode) {
389 os << " <Snapshot>\n";
390 writeFrameProperties(os,
391 info_->getSnapshotManager()->getCurrentSnapshot());
392 os << " <StuntDoubles>\n";
393 }
394
395 // every node prepares the dump lines for integrable objects belong to
396 // itself
397 std::string buffer;
398 for (mol = info_->beginMolecule(mi); mol != NULL;
399 mol = info_->nextMolecule(mi)) {
400 for (sd = mol->beginIntegrableObject(ii); sd != NULL;
401 sd = mol->nextIntegrableObject(ii)) {
402 buffer += prepareDumpLine(sd);
403 }
404 }
405
406 if (worldRank == primaryNode) {
407 os << buffer;
408
409 for (int i = 1; i < nProc; ++i) {
410 // tell processor i to start sending us data:
411
412 MPI_Bcast(&i, 1, MPI_INT, primaryNode, MPI_COMM_WORLD);
413
414 // receive the length of the string buffer that was
415 // prepared by processor i:
416 int recvLength;
417 MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD,
418 &istatus);
419
420 // create a buffer to receive the data
421 char* recvBuffer = new char[recvLength];
422 if (recvBuffer == NULL) {
423 } else {
424 // receive the data:
425 MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, MPI_ANY_TAG,
426 MPI_COMM_WORLD, &istatus);
427 // send it to the file:
428 os << recvBuffer;
429 // get rid of the receive buffer:
430 delete[] recvBuffer;
431 }
432 }
433 } else {
434 int sendBufferLength = buffer.size() + 1;
435 int myturn = 0;
436 for (int i = 1; i < nProc; ++i) {
437 // wait for the primary node to call our number:
438 MPI_Bcast(&myturn, 1, MPI_INT, primaryNode, MPI_COMM_WORLD);
439 if (myturn == worldRank) {
440 // send the length of our buffer:
441
442 MPI_Send(&sendBufferLength, 1, MPI_INT, primaryNode, 0,
443 MPI_COMM_WORLD);
444
445 // send our buffer:
446 MPI_Send((void*)buffer.c_str(), sendBufferLength, MPI_CHAR,
447 primaryNode, 0, MPI_COMM_WORLD);
448 }
449 }
450 }
451
452 if (worldRank == primaryNode) { os << " </StuntDoubles>\n"; }
453
454 if (doSiteData_) {
455 if (worldRank == primaryNode) { os << " <SiteData>\n"; }
456 buffer.clear();
457 for (mol = info_->beginMolecule(mi); mol != NULL;
458 mol = info_->nextMolecule(mi)) {
459 for (sd = mol->beginIntegrableObject(ii); sd != NULL;
460 sd = mol->nextIntegrableObject(ii)) {
461 int ioIndex = sd->getGlobalIntegrableObjectIndex();
462 // do one for the IO itself
463 buffer += prepareSiteLine(sd, ioIndex, 0);
464
465 if (sd->isRigidBody()) {
466 RigidBody* rb = static_cast<RigidBody*>(sd);
467 int siteIndex = 0;
468 for (Atom* atom = rb->beginAtom(ai); atom != NULL;
469 atom = rb->nextAtom(ai)) {
470 buffer += prepareSiteLine(atom, ioIndex, siteIndex);
471 siteIndex++;
472 }
473 }
474 }
475 }
476
477 if (worldRank == primaryNode) {
478 os << buffer;
479
480 for (int i = 1; i < nProc; ++i) {
481 // tell processor i to start sending us data:
482 MPI_Bcast(&i, 1, MPI_INT, primaryNode, MPI_COMM_WORLD);
483
484 // receive the length of the string buffer that was
485 // prepared by processor i:
486 int recvLength;
487 MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD,
488 &istatus);
489
490 // create a buffer to receive the data
491 char* recvBuffer = new char[recvLength];
492 if (recvBuffer == NULL) {
493 } else {
494 // receive the data:
495 MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, MPI_ANY_TAG,
496 MPI_COMM_WORLD, &istatus);
497 // send it to the file:
498 os << recvBuffer;
499 // get rid of the receive buffer:
500 delete[] recvBuffer;
501 }
502 }
503 } else {
504 int sendBufferLength = buffer.size() + 1;
505 int myturn = 0;
506 for (int i = 1; i < nProc; ++i) {
507 // wait for the primary node to call our number:
508 MPI_Bcast(&myturn, 1, MPI_INT, primaryNode, MPI_COMM_WORLD);
509 if (myturn == worldRank) {
510 // send the length of our buffer:
511 MPI_Send(&sendBufferLength, 1, MPI_INT, primaryNode, 0,
512 MPI_COMM_WORLD);
513 // send our buffer:
514 MPI_Send((void*)buffer.c_str(), sendBufferLength, MPI_CHAR,
515 primaryNode, 0, MPI_COMM_WORLD);
516 }
517 }
518 }
519
520 if (worldRank == primaryNode) { os << " </SiteData>\n"; }
521 }
522
523 if (worldRank == primaryNode) {
524 os << " </Snapshot>\n";
525 os.flush();
526 os.rdbuf()->pubsync();
527 }
528
529#endif // is_mpi
530 }
531
532 std::string DumpWriter::prepareDumpLine(StuntDouble* sd) {
533 int index = sd->getGlobalIntegrableObjectIndex();
534 std::string type("pv");
535 std::string line;
536 char tempBuffer[4096];
537
538 Vector3d pos;
539 Vector3d vel;
540 pos = sd->getPos();
541
542 if (std::isinf(pos[0]) || std::isnan(pos[0]) || std::isinf(pos[1]) ||
543 std::isnan(pos[1]) || std::isinf(pos[2]) || std::isnan(pos[2])) {
544 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
545 "DumpWriter detected a numerical error writing the position"
546 " for object %d",
547 index);
548 painCave.isFatal = 1;
549 simError();
550 }
551
552 vel = sd->getVel();
553
554 if (std::isinf(vel[0]) || std::isnan(vel[0]) || std::isinf(vel[1]) ||
555 std::isnan(vel[1]) || std::isinf(vel[2]) || std::isnan(vel[2])) {
556 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
557 "DumpWriter detected a numerical error writing the velocity"
558 " for object %d",
559 index);
560 painCave.isFatal = 1;
561 simError();
562 }
563
564 snprintf(tempBuffer, 4096, "%18.10g %18.10g %18.10g %13e %13e %13e", pos[0],
565 pos[1], pos[2], vel[0], vel[1], vel[2]);
566 line += tempBuffer;
567
568 if (sd->isDirectional()) {
569 type += "qj";
570 Quat4d q;
571 Vector3d ji;
572 q = sd->getQ();
573
574 if (std::isinf(q[0]) || std::isnan(q[0]) || std::isinf(q[1]) ||
575 std::isnan(q[1]) || std::isinf(q[2]) || std::isnan(q[2]) ||
576 std::isinf(q[3]) || std::isnan(q[3])) {
577 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
578 "DumpWriter detected a numerical error writing the quaternion"
579 " for object %d",
580 index);
581 painCave.isFatal = 1;
582 simError();
583 }
584
585 ji = sd->getJ();
586
587 if (std::isinf(ji[0]) || std::isnan(ji[0]) || std::isinf(ji[1]) ||
588 std::isnan(ji[1]) || std::isinf(ji[2]) || std::isnan(ji[2])) {
589 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
590 "DumpWriter detected a numerical error writing the angular"
591 " momentum for object %d",
592 index);
593 painCave.isFatal = 1;
594 simError();
595 }
596
597 snprintf(tempBuffer, 4096, " %13e %13e %13e %13e %13e %13e %13e", q[0],
598 q[1], q[2], q[3], ji[0], ji[1], ji[2]);
599 line += tempBuffer;
600 }
601
602 if (needForceVector_) {
603 type += "f";
604 Vector3d frc = sd->getFrc();
605 if (std::isinf(frc[0]) || std::isnan(frc[0]) || std::isinf(frc[1]) ||
606 std::isnan(frc[1]) || std::isinf(frc[2]) || std::isnan(frc[2])) {
607 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
608 "DumpWriter detected a numerical error writing the force"
609 " for object %d",
610 index);
611 painCave.isFatal = 1;
612 simError();
613 }
614 snprintf(tempBuffer, 4096, " %13e %13e %13e", frc[0], frc[1], frc[2]);
615 line += tempBuffer;
616
617 if (sd->isDirectional()) {
618 type += "t";
619 Vector3d trq = sd->getTrq();
620 if (std::isinf(trq[0]) || std::isnan(trq[0]) || std::isinf(trq[1]) ||
621 std::isnan(trq[1]) || std::isinf(trq[2]) || std::isnan(trq[2])) {
622 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
623 "DumpWriter detected a numerical error writing the torque"
624 " for object %d",
625 index);
626 painCave.isFatal = 1;
627 simError();
628 }
629 snprintf(tempBuffer, 4096, " %13e %13e %13e", trq[0], trq[1], trq[2]);
630 line += tempBuffer;
631 }
632 }
633
634 snprintf(tempBuffer, 4096, "%10d %7s %s\n", index, type.c_str(),
635 line.c_str());
636 return std::string(tempBuffer);
637 }
638
639 std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex,
640 int siteIndex) {
641 int asl = info_->getSnapshotManager()->getAtomStorageLayout();
642 int rbsl = info_->getSnapshotManager()->getRigidBodyStorageLayout();
643 int sl {};
644
645 std::string id;
646 std::string type;
647 std::string line;
648 char tempBuffer[4096];
649
650 if (sd->isRigidBody()) {
651 sl = rbsl;
652 snprintf(tempBuffer, 4096, "%10d ", ioIndex);
653 id = std::string(tempBuffer);
654 } else {
655 sl = asl;
656 snprintf(tempBuffer, 4096, "%10d %10d", ioIndex, siteIndex);
657 id = std::string(tempBuffer);
658 }
659
660 if (needFlucQ_) {
661 if (sl & DataStorage::dslFlucQPosition) {
662 type += "c";
663 RealType fqPos = sd->getFlucQPos();
664 if (std::isinf(fqPos) || std::isnan(fqPos)) {
665 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
666 "DumpWriter detected a numerical error writing the"
667 " fluctuating charge for object %s",
668 id.c_str());
669 painCave.isFatal = 1;
670 simError();
671 }
672 snprintf(tempBuffer, 4096, " %13e ", fqPos);
673 line += tempBuffer;
674 }
675
676 if (sl & DataStorage::dslFlucQVelocity) {
677 type += "w";
678 RealType fqVel = sd->getFlucQVel();
679 if (std::isinf(fqVel) || std::isnan(fqVel)) {
680 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
681 "DumpWriter detected a numerical error writing the"
682 " fluctuating charge velocity for object %s",
683 id.c_str());
684 painCave.isFatal = 1;
685 simError();
686 }
687 snprintf(tempBuffer, 4096, " %13e ", fqVel);
688 line += tempBuffer;
689 }
690
691 if (needForceVector_) {
692 if (sl & DataStorage::dslFlucQForce) {
693 type += "g";
694 RealType fqFrc = sd->getFlucQFrc();
695 if (std::isinf(fqFrc) || std::isnan(fqFrc)) {
696 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
697 "DumpWriter detected a numerical error writing the"
698 " fluctuating charge force for object %s",
699 id.c_str());
700 painCave.isFatal = 1;
701 simError();
702 }
703 snprintf(tempBuffer, 4096, " %13e ", fqFrc);
704 line += tempBuffer;
705 }
706 }
707 }
708
709 if (needElectricField_) {
710 if (sl & DataStorage::dslElectricField) {
711 type += "e";
712 Vector3d eField = sd->getElectricField();
713 if (std::isinf(eField[0]) || std::isnan(eField[0]) ||
714 std::isinf(eField[1]) || std::isnan(eField[1]) ||
715 std::isinf(eField[2]) || std::isnan(eField[2])) {
716 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
717 "DumpWriter detected a numerical error writing the electric"
718 " field for object %s",
719 id.c_str());
720 painCave.isFatal = 1;
721 simError();
722 }
723 snprintf(tempBuffer, 4096, " %13e %13e %13e", eField[0], eField[1],
724 eField[2]);
725 line += tempBuffer;
726 }
727 }
728
729 if (needSitePotential_) {
730 if (sl & DataStorage::dslSitePotential) {
731 type += "s";
732 RealType sPot = sd->getSitePotential();
733 if (std::isinf(sPot) || std::isnan(sPot)) {
734 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
735 "DumpWriter detected a numerical error writing the"
736 " site potential for object %s",
737 id.c_str());
738 painCave.isFatal = 1;
739 simError();
740 }
741 snprintf(tempBuffer, 4096, " %13e ", sPot);
742 line += tempBuffer;
743 }
744 }
745
746 if (needParticlePot_) {
747 if (sl & DataStorage::dslParticlePot) {
748 type += "u";
749 RealType particlePot = sd->getParticlePot();
750 if (std::isinf(particlePot) || std::isnan(particlePot)) {
751 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
752 "DumpWriter detected a numerical error writing the particle "
753 " potential for object %s",
754 id.c_str());
755 painCave.isFatal = 1;
756 simError();
757 }
758 snprintf(tempBuffer, 4096, " %13e", particlePot);
759 line += tempBuffer;
760 }
761 }
762
763 if (needDensity_) {
764 if (sl & DataStorage::dslDensity) {
765 type += "d";
766 RealType density = sd->getDensity();
767 if (std::isinf(density) || std::isnan(density)) {
768 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
769 "DumpWriter detected a numerical error writing the density "
770 " for object %s",
771 id.c_str());
772 painCave.isFatal = 1;
773 simError();
774 }
775 snprintf(tempBuffer, 4096, " %13e", density);
776 line += tempBuffer;
777 }
778 }
779
780 snprintf(tempBuffer, 4096, "%s %7s %s\n", id.c_str(), type.c_str(),
781 line.c_str());
782 return std::string(tempBuffer);
783 }
784
785 void DumpWriter::writeDump() { writeFrame(*dumpFile_); }
786
787 void DumpWriter::writeEor() {
788 std::ostream* eorStream = NULL;
789
790#ifdef IS_MPI
791 if (worldRank == 0) {
792#endif // is_mpi
793
794 eorStream = createOStream(eorFilename_);
795
796#ifdef IS_MPI
797 }
798#endif
799
800 writeFrame(*eorStream);
801
802#ifdef IS_MPI
803 if (worldRank == 0) {
804#endif
805
806 writeClosing(*eorStream);
807 delete eorStream;
808
809#ifdef IS_MPI
810 }
811#endif // is_mpi
812 }
813
814 void DumpWriter::writeDumpAndEor() {
815 std::vector<std::streambuf*> buffers;
816 std::ostream* eorStream = NULL;
817#ifdef IS_MPI
818 if (worldRank == 0) {
819#endif // is_mpi
820 buffers.push_back(dumpFile_->rdbuf());
821 eorStream = createOStream(eorFilename_);
822 buffers.push_back(eorStream->rdbuf());
823#ifdef IS_MPI
824 }
825#endif // is_mpi
826
827 TeeBuf tbuf(buffers.begin(), buffers.end());
828 std::ostream os(&tbuf);
829 writeFrame(os);
830
831#ifdef IS_MPI
832 if (worldRank == 0) {
833#endif // is_mpi
834 writeClosing(*eorStream);
835 delete eorStream;
836#ifdef IS_MPI
837 }
838#endif // is_mpi
839 }
840
841 std::ostream* DumpWriter::createOStream(const std::string& filename) {
842 std::ostream* newOStream;
843#ifdef HAVE_ZLIB
844 if (needCompression_) {
845 newOStream = new ogzstream(filename.c_str());
846 } else {
847 newOStream = new std::ofstream(filename.c_str());
848 }
849#else
850 newOStream = new std::ofstream(filename.c_str());
851#endif
852 // write out MetaData first
853 (*newOStream) << "<OpenMD version=2>" << std::endl;
854 (*newOStream) << " <MetaData>" << std::endl;
855 (*newOStream) << info_->getRawMetaData();
856 (*newOStream) << " </MetaData>" << std::endl;
857 return newOStream;
858 }
859
860 void DumpWriter::writeClosing(std::ostream& os) {
861 os << "</OpenMD>\n";
862 os.flush();
863 }
864
865} // namespace OpenMD
OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Definition ActionCorrFunc.cpp:60