ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/io/DumpReader.cpp
Revision: 1794
Committed: Thu Sep 6 19:44:06 2012 UTC (12 years, 7 months ago) by gezelter
File size: 21439 byte(s)
Log Message:
Merging some of the trunk changes back to the development branch,
cleaning up a datastorage bug

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

Properties

Name Value
svn:keywords Author Id Revision Date