ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/io/DumpReader.cpp
Revision: 1764
Committed: Tue Jul 3 18:32:27 2012 UTC (12 years, 10 months ago) by gezelter
File size: 21617 byte(s)
Log Message:
Refactored Snapshot and Stats to use the Accumulator classes.  Collected
a number of methods into Thermo that belonged there.

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
80 if (inFile_->fail()) {
81 sprintf(painCave.errMsg,
82 "DumpReader: Cannot open file: %s\n",
83 filename_.c_str());
84 painCave.isFatal = 1;
85 simError();
86 }
87
88 #ifdef IS_MPI
89
90 }
91
92 strcpy(checkPointMsg, "Dump file opened for reading successfully.");
93 errorCheckPoint();
94
95 #endif
96
97 return;
98 }
99
100 DumpReader::~DumpReader() {
101
102 #ifdef IS_MPI
103
104 if (worldRank == 0) {
105 #endif
106
107 delete inFile_;
108
109 #ifdef IS_MPI
110
111 }
112
113 strcpy(checkPointMsg, "Dump file closed successfully.");
114 errorCheckPoint();
115
116 #endif
117
118 return;
119 }
120
121 int DumpReader::getNFrames(void) {
122
123 if (!isScanned_)
124 scanFile();
125
126 return nframes_;
127 }
128
129 void DumpReader::scanFile(void) {
130 int lineNo = 0;
131 std::streampos prevPos;
132 std::streampos currPos;
133
134 #ifdef IS_MPI
135
136 if (worldRank == 0) {
137 #endif // is_mpi
138
139 currPos = inFile_->tellg();
140 prevPos = currPos;
141 bool foundOpenSnapshotTag = false;
142 bool foundClosedSnapshotTag = false;
143 bool foundOpenSiteDataTag = false;
144 while(inFile_->getline(buffer, bufferSize)) {
145 ++lineNo;
146
147 std::string line = buffer;
148 currPos = inFile_->tellg();
149 if (line.find("<Snapshot>")!= std::string::npos) {
150 if (foundOpenSnapshotTag) {
151 sprintf(painCave.errMsg,
152 "DumpReader:<Snapshot> is multiply nested at line %d in %s \n", lineNo,
153 filename_.c_str());
154 painCave.isFatal = 1;
155 simError();
156 }
157 foundOpenSnapshotTag = true;
158 foundClosedSnapshotTag = false;
159 framePos_.push_back(prevPos);
160
161 } else if (line.find("</Snapshot>") != std::string::npos){
162 if (!foundOpenSnapshotTag) {
163 sprintf(painCave.errMsg,
164 "DumpReader:</Snapshot> appears before <Snapshot> at line %d in %s \n", lineNo,
165 filename_.c_str());
166 painCave.isFatal = 1;
167 simError();
168 }
169
170 if (foundClosedSnapshotTag) {
171 sprintf(painCave.errMsg,
172 "DumpReader:</Snapshot> appears multiply nested at line %d in %s \n", lineNo,
173 filename_.c_str());
174 painCave.isFatal = 1;
175 simError();
176 }
177 foundClosedSnapshotTag = true;
178 foundOpenSnapshotTag = false;
179 }
180 prevPos = currPos;
181 }
182
183 // only found <Snapshot> for the last frame means the file is corrupted, we should discard
184 // it and give a warning message
185 if (foundOpenSnapshotTag) {
186 sprintf(painCave.errMsg,
187 "DumpReader: last frame in %s is invalid\n", filename_.c_str());
188 painCave.isFatal = 0;
189 simError();
190 framePos_.pop_back();
191 }
192
193 nframes_ = framePos_.size();
194
195 if (nframes_ == 0) {
196 sprintf(painCave.errMsg,
197 "DumpReader: %s does not contain a valid frame\n", filename_.c_str());
198 painCave.isFatal = 1;
199 simError();
200 }
201 #ifdef IS_MPI
202 }
203
204 MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
205
206 #endif // is_mpi
207
208 isScanned_ = true;
209 }
210
211 void DumpReader::readFrame(int whichFrame) {
212 if (!isScanned_)
213 scanFile();
214
215 int storageLayout = info_->getSnapshotManager()->getStorageLayout();
216
217 if (storageLayout & DataStorage::dslPosition) {
218 needPos_ = true;
219 } else {
220 needPos_ = false;
221 }
222
223 if (storageLayout & DataStorage::dslVelocity) {
224 needVel_ = true;
225 } else {
226 needVel_ = false;
227 }
228
229 if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
230 needQuaternion_ = true;
231 } else {
232 needQuaternion_ = false;
233 }
234
235 if (storageLayout & DataStorage::dslAngularMomentum) {
236 needAngMom_ = true;
237 } else {
238 needAngMom_ = false;
239 }
240
241 readSet(whichFrame);
242
243 if (needCOMprops_) {
244 Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
245 Thermo thermo(info_);
246 Vector3d com;
247
248 if (needPos_ && needVel_) {
249 Vector3d comvel;
250 Vector3d comw;
251 thermo.getComAll(com, comvel);
252 comw = thermo.getAngularMomentum();
253 } else {
254 com = thermo.getCom();
255 }
256 }
257 }
258
259 void DumpReader::readSet(int whichFrame) {
260 std::string line;
261
262 #ifndef IS_MPI
263 inFile_->clear();
264 inFile_->seekg(framePos_[whichFrame]);
265
266 std::istream& inputStream = *inFile_;
267
268 #else
269 int masterNode = 0;
270 std::stringstream sstream;
271 if (worldRank == masterNode) {
272 std::string sendBuffer;
273
274 inFile_->clear();
275 inFile_->seekg(framePos_[whichFrame]);
276
277 while (inFile_->getline(buffer, bufferSize)) {
278
279 line = buffer;
280 sendBuffer += line;
281 sendBuffer += '\n';
282 if (line.find("</Snapshot>") != std::string::npos) {
283 break;
284 }
285 }
286
287 int sendBufferSize = sendBuffer.size();
288 MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
289 MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
290
291 sstream.str(sendBuffer);
292 } else {
293 int sendBufferSize;
294 MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
295 char * recvBuffer = new char[sendBufferSize+1];
296 assert(recvBuffer);
297 recvBuffer[sendBufferSize] = '\0';
298 MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
299 sstream.str(recvBuffer);
300 delete [] recvBuffer;
301 }
302
303 std::istream& inputStream = sstream;
304 #endif
305
306 inputStream.getline(buffer, bufferSize);
307
308 line = buffer;
309 if (line.find("<Snapshot>") == std::string::npos) {
310 sprintf(painCave.errMsg,
311 "DumpReader Error: can not find <Snapshot>\n");
312 painCave.isFatal = 1;
313 simError();
314 }
315
316 //read frameData
317 readFrameProperties(inputStream);
318
319 //read StuntDoubles
320 readStuntDoubles(inputStream);
321
322 inputStream.getline(buffer, bufferSize);
323 line = buffer;
324
325 if (line.find("<SiteData>") != std::string::npos) {
326 //read SiteData
327 readSiteData(inputStream);
328 } else {
329 if (line.find("</Snapshot>") == std::string::npos) {
330 sprintf(painCave.errMsg,
331 "DumpReader Error: can not find </Snapshot>\n");
332 painCave.isFatal = 1;
333 simError();
334 }
335 }
336 }
337
338 void DumpReader::parseDumpLine(const std::string& line) {
339
340
341 StringTokenizer tokenizer(line);
342 int nTokens;
343
344 nTokens = tokenizer.countTokens();
345
346 if (nTokens < 2) {
347 sprintf(painCave.errMsg,
348 "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
349 painCave.isFatal = 1;
350 simError();
351 }
352
353 int index = tokenizer.nextTokenAsInt();
354
355 StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
356
357 if (integrableObject == NULL) {
358 return;
359 }
360 std::string type = tokenizer.nextToken();
361 int size = type.size();
362
363 size_t found;
364
365 if (needPos_) {
366 found = type.find("p");
367 if (found == std::string::npos) {
368 sprintf(painCave.errMsg,
369 "DumpReader Error: StuntDouble %d has no Position\n"
370 "\tField (\"p\") specified.\n%s\n", index,
371 line.c_str());
372 painCave.isFatal = 1;
373 simError();
374 }
375 }
376
377 if (integrableObject->isDirectional()) {
378 if (needQuaternion_) {
379 found = type.find("q");
380 if (found == std::string::npos) {
381 sprintf(painCave.errMsg,
382 "DumpReader Error: Directional StuntDouble %d has no\n"
383 "\tQuaternion Field (\"q\") specified.\n%s\n", index,
384 line.c_str());
385 painCave.isFatal = 1;
386 simError();
387 }
388 }
389 }
390
391 for(int i = 0; i < size; ++i) {
392 switch(type[i]) {
393
394 case 'p': {
395 Vector3d pos;
396 pos[0] = tokenizer.nextTokenAsDouble();
397 pos[1] = tokenizer.nextTokenAsDouble();
398 pos[2] = tokenizer.nextTokenAsDouble();
399 if (needPos_) {
400 integrableObject->setPos(pos);
401 }
402 break;
403 }
404 case 'v' : {
405 Vector3d vel;
406 vel[0] = tokenizer.nextTokenAsDouble();
407 vel[1] = tokenizer.nextTokenAsDouble();
408 vel[2] = tokenizer.nextTokenAsDouble();
409 if (needVel_) {
410 integrableObject->setVel(vel);
411 }
412 break;
413 }
414
415 case 'q' : {
416 Quat4d q;
417 if (integrableObject->isDirectional()) {
418
419 q[0] = tokenizer.nextTokenAsDouble();
420 q[1] = tokenizer.nextTokenAsDouble();
421 q[2] = tokenizer.nextTokenAsDouble();
422 q[3] = tokenizer.nextTokenAsDouble();
423
424 RealType qlen = q.length();
425 if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
426
427 sprintf(painCave.errMsg,
428 "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
429 painCave.isFatal = 1;
430 simError();
431
432 }
433
434 q.normalize();
435 if (needQuaternion_) {
436 integrableObject->setQ(q);
437 }
438 }
439 break;
440 }
441 case 'j' : {
442 Vector3d ji;
443 if (integrableObject->isDirectional()) {
444 ji[0] = tokenizer.nextTokenAsDouble();
445 ji[1] = tokenizer.nextTokenAsDouble();
446 ji[2] = tokenizer.nextTokenAsDouble();
447 if (needAngMom_) {
448 integrableObject->setJ(ji);
449 }
450 }
451 break;
452 }
453 case 'f': {
454
455 Vector3d force;
456 force[0] = tokenizer.nextTokenAsDouble();
457 force[1] = tokenizer.nextTokenAsDouble();
458 force[2] = tokenizer.nextTokenAsDouble();
459 integrableObject->setFrc(force);
460 break;
461 }
462 case 't' : {
463
464 Vector3d torque;
465 torque[0] = tokenizer.nextTokenAsDouble();
466 torque[1] = tokenizer.nextTokenAsDouble();
467 torque[2] = tokenizer.nextTokenAsDouble();
468 integrableObject->setTrq(torque);
469 break;
470 }
471 case 'u' : {
472
473 RealType particlePot;
474 particlePot = tokenizer.nextTokenAsDouble();
475 integrableObject->setParticlePot(particlePot);
476 break;
477 }
478 case 'c' : {
479
480 RealType flucQPos;
481 flucQPos = tokenizer.nextTokenAsDouble();
482 integrableObject->setFlucQPos(flucQPos);
483 break;
484 }
485 case 'w' : {
486
487 RealType flucQVel;
488 flucQVel = tokenizer.nextTokenAsDouble();
489 integrableObject->setFlucQVel(flucQVel);
490 break;
491 }
492 case 'g' : {
493
494 RealType flucQFrc;
495 flucQFrc = tokenizer.nextTokenAsDouble();
496 integrableObject->setFlucQFrc(flucQFrc);
497 break;
498 }
499 case 'e' : {
500
501 Vector3d eField;
502 eField[0] = tokenizer.nextTokenAsDouble();
503 eField[1] = tokenizer.nextTokenAsDouble();
504 eField[2] = tokenizer.nextTokenAsDouble();
505 integrableObject->setElectricField(eField);
506 break;
507 }
508 default: {
509 sprintf(painCave.errMsg,
510 "DumpReader Error: %s is an unrecognized type\n", type.c_str());
511 painCave.isFatal = 1;
512 simError();
513 break;
514 }
515
516 }
517 }
518
519 }
520
521
522 void DumpReader::parseSiteLine(const std::string& line) {
523
524 StringTokenizer tokenizer(line);
525 int nTokens;
526
527 nTokens = tokenizer.countTokens();
528
529 if (nTokens < 2) {
530 sprintf(painCave.errMsg,
531 "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
532 painCave.isFatal = 1;
533 simError();
534 }
535
536 /**
537 * The first token is the global integrable object index.
538 */
539
540 int index = tokenizer.nextTokenAsInt();
541 StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
542 if (integrableObject == NULL) {
543 return;
544 }
545 StuntDouble* sd = integrableObject;
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*>(integrableObject);
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