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