ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/io/DumpReader.cpp
Revision: 1793
Committed: Fri Aug 31 21:16:10 2012 UTC (12 years, 8 months ago) by gezelter
File size: 21341 byte(s)
Log Message:
Cleaning up some warning messages on linux

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

Properties

Name Value
svn:keywords Author Id Revision Date