ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/io/DumpReader.cpp
Revision: 1897
Committed: Wed Jul 3 20:21:00 2013 UTC (11 years, 9 months ago) by gezelter
File size: 21430 byte(s)
Log Message:
Fixing a DumpReader bug (I think) for site data on integrable objects that are not RigidBodies.

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, 234107 (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 #define _LARGEFILE_SOURCE64
44 #define _FILE_OFFSET_BITS 64
45
46 #include <sys/types.h>
47 #include <sys/stat.h>
48
49 #include <iostream>
50 #include <math.h>
51
52 #include <stdio.h>
53 #include <stdlib.h>
54 #include <string.h>
55
56 #include "io/DumpReader.hpp"
57 #include "primitives/Molecule.hpp"
58 #include "utils/simError.h"
59 #include "utils/MemoryUtils.hpp"
60 #include "utils/StringTokenizer.hpp"
61 #include "brains/Thermo.hpp"
62
63 #ifdef IS_MPI
64 #include <mpi.h>
65 #endif
66
67
68 namespace OpenMD {
69
70 DumpReader::DumpReader(SimInfo* info, const std::string& filename)
71 : info_(info), filename_(filename), isScanned_(false), nframes_(0), needCOMprops_(false) {
72
73 #ifdef IS_MPI
74
75 if (worldRank == 0) {
76 #endif
77
78 inFile_ = new std::ifstream(filename_.c_str(),
79 ifstream::in | ifstream::binary);
80
81 if (inFile_->fail()) {
82 sprintf(painCave.errMsg,
83 "DumpReader: Cannot open file: %s\n",
84 filename_.c_str());
85 painCave.isFatal = 1;
86 simError();
87 }
88
89 #ifdef IS_MPI
90
91 }
92
93 strcpy(checkPointMsg, "Dump file opened for reading successfully.");
94 errorCheckPoint();
95
96 #endif
97
98 return;
99 }
100
101 DumpReader::~DumpReader() {
102
103 #ifdef IS_MPI
104
105 if (worldRank == 0) {
106 #endif
107
108 delete inFile_;
109
110 #ifdef IS_MPI
111
112 }
113
114 strcpy(checkPointMsg, "Dump file closed successfully.");
115 errorCheckPoint();
116
117 #endif
118
119 return;
120 }
121
122 int DumpReader::getNFrames(void) {
123
124 if (!isScanned_)
125 scanFile();
126
127 return nframes_;
128 }
129
130 void DumpReader::scanFile(void) {
131
132 std::streampos prevPos;
133 std::streampos currPos;
134
135 #ifdef IS_MPI
136
137 if (worldRank == 0) {
138 #endif // is_mpi
139
140 currPos = inFile_->tellg();
141 prevPos = currPos;
142 bool foundOpenSnapshotTag = false;
143 bool foundClosedSnapshotTag = false;
144
145 int lineNo = 0;
146 while(inFile_->getline(buffer, bufferSize)) {
147 ++lineNo;
148
149 std::string line = buffer;
150 currPos = inFile_->tellg();
151 if (line.find("<Snapshot>")!= std::string::npos) {
152 if (foundOpenSnapshotTag) {
153 sprintf(painCave.errMsg,
154 "DumpReader:<Snapshot> is multiply nested at line %d in %s \n", lineNo,
155 filename_.c_str());
156 painCave.isFatal = 1;
157 simError();
158 }
159 foundOpenSnapshotTag = true;
160 foundClosedSnapshotTag = false;
161 framePos_.push_back(prevPos);
162
163 } else if (line.find("</Snapshot>") != std::string::npos){
164 if (!foundOpenSnapshotTag) {
165 sprintf(painCave.errMsg,
166 "DumpReader:</Snapshot> appears before <Snapshot> at line %d in %s \n", lineNo,
167 filename_.c_str());
168 painCave.isFatal = 1;
169 simError();
170 }
171
172 if (foundClosedSnapshotTag) {
173 sprintf(painCave.errMsg,
174 "DumpReader:</Snapshot> appears multiply nested at line %d in %s \n", lineNo,
175 filename_.c_str());
176 painCave.isFatal = 1;
177 simError();
178 }
179 foundClosedSnapshotTag = true;
180 foundOpenSnapshotTag = false;
181 }
182 prevPos = currPos;
183 }
184
185 // only found <Snapshot> for the last frame means the file is corrupted, we should discard
186 // it and give a warning message
187 if (foundOpenSnapshotTag) {
188 sprintf(painCave.errMsg,
189 "DumpReader: last frame in %s is invalid\n", filename_.c_str());
190 painCave.isFatal = 0;
191 simError();
192 framePos_.pop_back();
193 }
194
195 nframes_ = framePos_.size();
196
197 if (nframes_ == 0) {
198 sprintf(painCave.errMsg,
199 "DumpReader: %s does not contain a valid frame\n", filename_.c_str());
200 painCave.isFatal = 1;
201 simError();
202 }
203 #ifdef IS_MPI
204 }
205
206 MPI::COMM_WORLD.Bcast(&nframes_, 1, MPI::INT, 0);
207
208 #endif // is_mpi
209
210 isScanned_ = true;
211 }
212
213 void DumpReader::readFrame(int whichFrame) {
214 if (!isScanned_)
215 scanFile();
216
217 int storageLayout = info_->getSnapshotManager()->getStorageLayout();
218
219 if (storageLayout & DataStorage::dslPosition) {
220 needPos_ = true;
221 } else {
222 needPos_ = false;
223 }
224
225 if (storageLayout & DataStorage::dslVelocity) {
226 needVel_ = true;
227 } else {
228 needVel_ = false;
229 }
230
231 if (storageLayout & DataStorage::dslAmat ||
232 storageLayout & DataStorage::dslDipole ||
233 storageLayout & DataStorage::dslQuadrupole) {
234 needQuaternion_ = true;
235 } else {
236 needQuaternion_ = false;
237 }
238
239 if (storageLayout & DataStorage::dslAngularMomentum) {
240 needAngMom_ = true;
241 } else {
242 needAngMom_ = false;
243 }
244
245 readSet(whichFrame);
246
247 if (needCOMprops_) {
248 Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
249 Thermo thermo(info_);
250 Vector3d com;
251
252 if (needPos_ && needVel_) {
253 Vector3d comvel;
254 Vector3d comw;
255 thermo.getComAll(com, comvel);
256 comw = thermo.getAngularMomentum();
257 } else {
258 com = thermo.getCom();
259 }
260 }
261 }
262
263 void DumpReader::readSet(int whichFrame) {
264 std::string line;
265
266 #ifndef IS_MPI
267 inFile_->clear();
268 inFile_->seekg(framePos_[whichFrame]);
269
270 std::istream& inputStream = *inFile_;
271
272 #else
273 int masterNode = 0;
274 std::stringstream sstream;
275 if (worldRank == masterNode) {
276 std::string sendBuffer;
277
278 inFile_->clear();
279 inFile_->seekg(framePos_[whichFrame]);
280
281 while (inFile_->getline(buffer, bufferSize)) {
282
283 line = buffer;
284 sendBuffer += line;
285 sendBuffer += '\n';
286 if (line.find("</Snapshot>") != std::string::npos) {
287 break;
288 }
289 }
290
291 int sendBufferSize = sendBuffer.size();
292 MPI::COMM_WORLD.Bcast(&sendBufferSize, 1, MPI::INT, masterNode);
293 MPI::COMM_WORLD.Bcast((void *)sendBuffer.c_str(), sendBufferSize,
294 MPI::CHAR, masterNode);
295
296 sstream.str(sendBuffer);
297 } else {
298 int sendBufferSize;
299 MPI::COMM_WORLD.Bcast(&sendBufferSize, 1, MPI::INT, masterNode);
300 char * recvBuffer = new char[sendBufferSize+1];
301 assert(recvBuffer);
302 recvBuffer[sendBufferSize] = '\0';
303 MPI::COMM_WORLD.Bcast(recvBuffer, sendBufferSize, MPI::CHAR, masterNode);
304 sstream.str(recvBuffer);
305 delete [] recvBuffer;
306 }
307
308 std::istream& inputStream = sstream;
309 #endif
310
311 inputStream.getline(buffer, bufferSize);
312
313 line = buffer;
314 if (line.find("<Snapshot>") == std::string::npos) {
315 sprintf(painCave.errMsg,
316 "DumpReader Error: can not find <Snapshot>\n");
317 painCave.isFatal = 1;
318 simError();
319 }
320
321 //read frameData
322 readFrameProperties(inputStream);
323
324 //read StuntDoubles
325 readStuntDoubles(inputStream);
326
327 inputStream.getline(buffer, bufferSize);
328 line = buffer;
329
330 if (line.find("<SiteData>") != std::string::npos) {
331 //read SiteData
332 readSiteData(inputStream);
333 } else {
334 if (line.find("</Snapshot>") == std::string::npos) {
335 sprintf(painCave.errMsg,
336 "DumpReader Error: can not find </Snapshot>\n");
337 painCave.isFatal = 1;
338 simError();
339 }
340 }
341 }
342
343 void DumpReader::parseDumpLine(const std::string& line) {
344
345
346 StringTokenizer tokenizer(line);
347 int nTokens;
348
349 nTokens = tokenizer.countTokens();
350
351 if (nTokens < 2) {
352 sprintf(painCave.errMsg,
353 "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
354 painCave.isFatal = 1;
355 simError();
356 }
357
358 int index = tokenizer.nextTokenAsInt();
359
360 StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
361
362 if (sd == NULL) {
363 return;
364 }
365 std::string type = tokenizer.nextToken();
366 int size = type.size();
367
368 size_t found;
369
370 if (needPos_) {
371 found = type.find("p");
372 if (found == std::string::npos) {
373 sprintf(painCave.errMsg,
374 "DumpReader Error: StuntDouble %d has no Position\n"
375 "\tField (\"p\") specified.\n%s\n", index,
376 line.c_str());
377 painCave.isFatal = 1;
378 simError();
379 }
380 }
381
382 if (sd->isDirectional()) {
383 if (needQuaternion_) {
384 found = type.find("q");
385 if (found == std::string::npos) {
386 sprintf(painCave.errMsg,
387 "DumpReader Error: Directional StuntDouble %d has no\n"
388 "\tQuaternion Field (\"q\") specified.\n%s\n", index,
389 line.c_str());
390 painCave.isFatal = 1;
391 simError();
392 }
393 }
394 }
395
396 for(int i = 0; i < size; ++i) {
397 switch(type[i]) {
398
399 case 'p': {
400 Vector3d pos;
401 pos[0] = tokenizer.nextTokenAsDouble();
402 pos[1] = tokenizer.nextTokenAsDouble();
403 pos[2] = tokenizer.nextTokenAsDouble();
404 if (needPos_) {
405 sd->setPos(pos);
406 }
407 break;
408 }
409 case 'v' : {
410 Vector3d vel;
411 vel[0] = tokenizer.nextTokenAsDouble();
412 vel[1] = tokenizer.nextTokenAsDouble();
413 vel[2] = tokenizer.nextTokenAsDouble();
414 if (needVel_) {
415 sd->setVel(vel);
416 }
417 break;
418 }
419
420 case 'q' : {
421 Quat4d q;
422 if (sd->isDirectional()) {
423
424 q[0] = tokenizer.nextTokenAsDouble();
425 q[1] = tokenizer.nextTokenAsDouble();
426 q[2] = tokenizer.nextTokenAsDouble();
427 q[3] = tokenizer.nextTokenAsDouble();
428
429 RealType qlen = q.length();
430 if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
431
432 sprintf(painCave.errMsg,
433 "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
434 painCave.isFatal = 1;
435 simError();
436
437 }
438
439 q.normalize();
440 if (needQuaternion_) {
441 sd->setQ(q);
442 }
443 }
444 break;
445 }
446 case 'j' : {
447 Vector3d ji;
448 if (sd->isDirectional()) {
449 ji[0] = tokenizer.nextTokenAsDouble();
450 ji[1] = tokenizer.nextTokenAsDouble();
451 ji[2] = tokenizer.nextTokenAsDouble();
452 if (needAngMom_) {
453 sd->setJ(ji);
454 }
455 }
456 break;
457 }
458 case 'f': {
459
460 Vector3d force;
461 force[0] = tokenizer.nextTokenAsDouble();
462 force[1] = tokenizer.nextTokenAsDouble();
463 force[2] = tokenizer.nextTokenAsDouble();
464 sd->setFrc(force);
465 break;
466 }
467 case 't' : {
468
469 Vector3d torque;
470 torque[0] = tokenizer.nextTokenAsDouble();
471 torque[1] = tokenizer.nextTokenAsDouble();
472 torque[2] = tokenizer.nextTokenAsDouble();
473 sd->setTrq(torque);
474 break;
475 }
476 case 'u' : {
477
478 RealType particlePot;
479 particlePot = tokenizer.nextTokenAsDouble();
480 sd->setParticlePot(particlePot);
481 break;
482 }
483 case 'c' : {
484
485 RealType flucQPos;
486 flucQPos = tokenizer.nextTokenAsDouble();
487 sd->setFlucQPos(flucQPos);
488 break;
489 }
490 case 'w' : {
491
492 RealType flucQVel;
493 flucQVel = tokenizer.nextTokenAsDouble();
494 sd->setFlucQVel(flucQVel);
495 break;
496 }
497 case 'g' : {
498
499 RealType flucQFrc;
500 flucQFrc = tokenizer.nextTokenAsDouble();
501 sd->setFlucQFrc(flucQFrc);
502 break;
503 }
504 case 'e' : {
505
506 Vector3d eField;
507 eField[0] = tokenizer.nextTokenAsDouble();
508 eField[1] = tokenizer.nextTokenAsDouble();
509 eField[2] = tokenizer.nextTokenAsDouble();
510 sd->setElectricField(eField);
511 break;
512 }
513 default: {
514 sprintf(painCave.errMsg,
515 "DumpReader Error: %s is an unrecognized type\n", type.c_str());
516 painCave.isFatal = 1;
517 simError();
518 break;
519 }
520
521 }
522 }
523
524 }
525
526
527 void DumpReader::parseSiteLine(const std::string& line) {
528
529 StringTokenizer tokenizer(line);
530 int nTokens;
531
532 nTokens = tokenizer.countTokens();
533
534 if (nTokens < 2) {
535 sprintf(painCave.errMsg,
536 "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
537 painCave.isFatal = 1;
538 simError();
539 }
540
541 /**
542 * The first token is the global integrable object index.
543 */
544
545 int index = tokenizer.nextTokenAsInt();
546 StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
547 if (sd == NULL) {
548 return;
549 }
550
551 /**
552 * Test to see if the next token is an integer or not. If not,
553 * we've got data on the integrable object itself. If there is an
554 * integer, we're parsing data for a site on a rigid body.
555 */
556
557 std::string indexTest = tokenizer.peekNextToken();
558 std::istringstream i(indexTest);
559 int siteIndex;
560 if (i >> siteIndex) {
561 // chew up this token and parse as an int:
562 siteIndex = tokenizer.nextTokenAsInt();
563
564 if (sd->isRigidBody()) {
565 RigidBody* rb = static_cast<RigidBody*>(sd);
566 sd = rb->getAtoms()[siteIndex];
567 }
568 }
569
570 /**
571 * The next token contains information on what follows.
572 */
573 std::string type = tokenizer.nextToken();
574 int size = type.size();
575
576 for(int i = 0; i < size; ++i) {
577 switch(type[i]) {
578
579 case 'u' : {
580
581 RealType particlePot;
582 particlePot = tokenizer.nextTokenAsDouble();
583 sd->setParticlePot(particlePot);
584 break;
585 }
586 case 'c' : {
587
588 RealType flucQPos;
589 flucQPos = tokenizer.nextTokenAsDouble();
590 sd->setFlucQPos(flucQPos);
591 break;
592 }
593 case 'w' : {
594
595 RealType flucQVel;
596 flucQVel = tokenizer.nextTokenAsDouble();
597 sd->setFlucQVel(flucQVel);
598 break;
599 }
600 case 'g' : {
601
602 RealType flucQFrc;
603 flucQFrc = tokenizer.nextTokenAsDouble();
604 sd->setFlucQFrc(flucQFrc);
605 break;
606 }
607 case 'e' : {
608
609 Vector3d eField;
610 eField[0] = tokenizer.nextTokenAsDouble();
611 eField[1] = tokenizer.nextTokenAsDouble();
612 eField[2] = tokenizer.nextTokenAsDouble();
613 sd->setElectricField(eField);
614 break;
615 }
616 default: {
617 sprintf(painCave.errMsg,
618 "DumpReader Error: %s is an unrecognized type\n", type.c_str());
619 painCave.isFatal = 1;
620 simError();
621 break;
622 }
623 }
624 }
625 }
626
627
628 void DumpReader::readStuntDoubles(std::istream& inputStream) {
629
630 inputStream.getline(buffer, bufferSize);
631 std::string line(buffer);
632
633 if (line.find("<StuntDoubles>") == std::string::npos) {
634 sprintf(painCave.errMsg,
635 "DumpReader Error: Missing <StuntDoubles>\n");
636 painCave.isFatal = 1;
637 simError();
638 }
639
640 while(inputStream.getline(buffer, bufferSize)) {
641 line = buffer;
642
643 if(line.find("</StuntDoubles>") != std::string::npos) {
644 break;
645 }
646
647 parseDumpLine(line);
648 }
649
650 }
651
652 void DumpReader::readSiteData(std::istream& inputStream) {
653
654 std::string line(buffer);
655
656 // We already found the starting <SiteData> tag or we wouldn't be
657 // here, so just start parsing until we get to the ending
658 // </SiteData> tag:
659
660 while(inputStream.getline(buffer, bufferSize)) {
661 line = buffer;
662
663 if(line.find("</SiteData>") != std::string::npos) {
664 break;
665 }
666
667 parseSiteLine(line);
668 }
669
670 }
671
672 void DumpReader::readFrameProperties(std::istream& inputStream) {
673
674 Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
675 inputStream.getline(buffer, bufferSize);
676 std::string line(buffer);
677
678 if (line.find("<FrameData>") == std::string::npos) {
679 sprintf(painCave.errMsg,
680 "DumpReader Error: Missing <FrameData>\n");
681 painCave.isFatal = 1;
682 simError();
683 }
684
685 while(inputStream.getline(buffer, bufferSize)) {
686 line = buffer;
687
688 if(line.find("</FrameData>") != std::string::npos) {
689 break;
690 }
691
692 StringTokenizer tokenizer(line, " ;\t\n\r{}:,");
693 if (!tokenizer.hasMoreTokens()) {
694 sprintf(painCave.errMsg,
695 "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
696 painCave.isFatal = 1;
697 simError();
698 }
699
700 std::string propertyName = tokenizer.nextToken();
701 if (propertyName == "Time") {
702 RealType currTime = tokenizer.nextTokenAsDouble();
703 s->setTime(currTime);
704 } else if (propertyName == "Hmat"){
705 Mat3x3d hmat;
706 hmat(0, 0) = tokenizer.nextTokenAsDouble();
707 hmat(0, 1) = tokenizer.nextTokenAsDouble();
708 hmat(0, 2) = tokenizer.nextTokenAsDouble();
709 hmat(1, 0) = tokenizer.nextTokenAsDouble();
710 hmat(1, 1) = tokenizer.nextTokenAsDouble();
711 hmat(1, 2) = tokenizer.nextTokenAsDouble();
712 hmat(2, 0) = tokenizer.nextTokenAsDouble();
713 hmat(2, 1) = tokenizer.nextTokenAsDouble();
714 hmat(2, 2) = tokenizer.nextTokenAsDouble();
715 s->setHmat(hmat);
716 } else if (propertyName == "Thermostat") {
717 pair<RealType, RealType> thermostat;
718 thermostat.first = tokenizer.nextTokenAsDouble();
719 thermostat.second = tokenizer.nextTokenAsDouble();
720 s->setThermostat(thermostat);
721 } else if (propertyName == "Barostat") {
722 Mat3x3d eta;
723 eta(0, 0) = tokenizer.nextTokenAsDouble();
724 eta(0, 1) = tokenizer.nextTokenAsDouble();
725 eta(0, 2) = tokenizer.nextTokenAsDouble();
726 eta(1, 0) = tokenizer.nextTokenAsDouble();
727 eta(1, 1) = tokenizer.nextTokenAsDouble();
728 eta(1, 2) = tokenizer.nextTokenAsDouble();
729 eta(2, 0) = tokenizer.nextTokenAsDouble();
730 eta(2, 1) = tokenizer.nextTokenAsDouble();
731 eta(2, 2) = tokenizer.nextTokenAsDouble();
732 s->setBarostat(eta);
733 } else {
734 sprintf(painCave.errMsg,
735 "DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str());
736 painCave.isFatal = 0;
737 simError();
738 }
739
740 }
741
742 }
743
744
745 }//end namespace OpenMD

Properties

Name Value
svn:keywords Author Id Revision Date