ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/io/DumpReader.cpp
(Generate patch)

Comparing:
trunk/src/io/DumpReader.cpp (file contents), Revision 398 by tim, Mon Mar 7 22:39:33 2005 UTC vs.
branches/development/src/io/DumpReader.cpp (file contents), Revision 1787 by gezelter, Wed Aug 29 18:13:11 2012 UTC

# Line 1 | Line 1
1 < /*
2 < * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
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. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
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.
# Line 37 | Line 28
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 < #define _LARGEFILE_SOURCE64
47 < #define _FILE_OFFSET_BITS 64
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 ||
230 >        storageLayout & DataStorage::dslDipole ||
231 >        storageLayout & DataStorage::dslQuadrupole) {
232 >      needQuaternion_ = true;
233 >    } else {
234 >      needQuaternion_ = false;
235 >    }
236 >    
237 >    if (storageLayout & DataStorage::dslAngularMomentum) {
238 >      needAngMom_ = true;
239 >    } else {
240 >      needAngMom_ = false;    
241 >    }
242 >    
243 >    readSet(whichFrame);
244  
245 < #include <sys/types.h>
246 < #include <sys/stat.h>
245 >    if (needCOMprops_) {
246 >      Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
247 >      Thermo thermo(info_);
248 >      Vector3d com;
249  
250 < #include <iostream>
251 < #include <math.h>
250 >      if (needPos_ && needVel_) {
251 >        Vector3d comvel;
252 >        Vector3d comw;
253 >        thermo.getComAll(com, comvel);
254 >        comw = thermo.getAngularMomentum();
255 >      } else {
256 >        com = thermo.getCom();
257 >      }                    
258 >    }
259 >  }
260 >  
261 >  void DumpReader::readSet(int whichFrame) {    
262 >    std::string line;
263  
264 < #include <stdio.h>
265 < #include <stdlib.h>
266 < #include <string.h>
264 > #ifndef IS_MPI
265 >    inFile_->clear();  
266 >    inFile_->seekg(framePos_[whichFrame]);
267  
268 < #include "io/DumpReader.hpp"
56 < #include "primitives/Molecule.hpp"
57 < #include "utils/simError.h"
58 < #include "utils/MemoryUtils.hpp"
59 < #include "utils/StringTokenizer.hpp"
268 >    std::istream& inputStream = *inFile_;    
269  
270 < #ifdef IS_MPI
270 > #else
271 >    int masterNode = 0;
272 >    std::stringstream sstream;
273 >    if (worldRank == masterNode) {
274 >      std::string sendBuffer;
275  
276 < #include <mpi.h>
277 < #define TAKE_THIS_TAG_CHAR 0
278 < #define TAKE_THIS_TAG_INT 1
276 >      inFile_->clear();  
277 >      inFile_->seekg(framePos_[whichFrame]);
278 >      
279 >      while (inFile_->getline(buffer, bufferSize)) {
280  
281 < #endif // is_mpi
281 >        line = buffer;
282 >        sendBuffer += line;
283 >        sendBuffer += '\n';
284 >        if (line.find("</Snapshot>") != std::string::npos) {
285 >          break;
286 >        }        
287 >      }
288  
289 +      int sendBufferSize = sendBuffer.size();
290 +      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);    
291 +      MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);    
292 +      
293 +      sstream.str(sendBuffer);
294 +    } else {
295 +      int sendBufferSize;
296 +      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);    
297 +      char * recvBuffer = new char[sendBufferSize+1];
298 +      assert(recvBuffer);
299 +      recvBuffer[sendBufferSize] = '\0';
300 +      MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);    
301 +      sstream.str(recvBuffer);
302 +      delete [] recvBuffer;
303 +    }      
304  
305 < namespace oopse {
71 <
72 < DumpReader::DumpReader(SimInfo* info, const std::string& filename)
73 <                     : info_(info), filename_(filename), isScanned_(false), nframes_(0) {
74 <
75 < #ifdef IS_MPI
76 <
77 <    if (worldRank == 0) {
305 >    std::istream& inputStream = sstream;  
306   #endif
307  
308 <        inFile_ = fopen(filename_.c_str(), "r");
308 >    inputStream.getline(buffer, bufferSize);
309  
310 <        if (inFile_ == NULL) {
311 <            sprintf(painCave.errMsg, "DumpReader: Cannot open file: %s\n", filename_.c_str());
312 <            painCave.isFatal = 1;
313 <            simError();
314 <        }
315 <
316 < #ifdef IS_MPI
310 >    line = buffer;
311 >    if (line.find("<Snapshot>") == std::string::npos) {
312 >      sprintf(painCave.errMsg,
313 >              "DumpReader Error: can not find <Snapshot>\n");
314 >      painCave.isFatal = 1;
315 >      simError();
316 >    }
317 >    
318 >    //read frameData
319 >    readFrameProperties(inputStream);
320  
321 <    }
321 >    //read StuntDoubles
322 >    readStuntDoubles(inputStream);    
323  
324 <    strcpy(checkPointMsg, "Dump file opened for reading successfully.");
325 <    MPIcheckPoint();
324 >    inputStream.getline(buffer, bufferSize);
325 >    line = buffer;
326  
327 < #endif
327 >    if (line.find("<SiteData>") != std::string::npos) {
328 >      //read SiteData
329 >      readSiteData(inputStream);        
330 >    } else {
331 >      if (line.find("</Snapshot>") == std::string::npos) {
332 >        sprintf(painCave.errMsg,
333 >                "DumpReader Error: can not find </Snapshot>\n");
334 >        painCave.isFatal = 1;
335 >        simError();
336 >      }        
337 >    }
338 >  }
339 >  
340 >  void DumpReader::parseDumpLine(const std::string& line) {
341  
342 <    return;
343 < }
342 >      
343 >    StringTokenizer tokenizer(line);
344 >    int nTokens;
345 >    
346 >    nTokens = tokenizer.countTokens();
347 >    
348 >    if (nTokens < 2) {  
349 >      sprintf(painCave.errMsg,
350 >              "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
351 >      painCave.isFatal = 1;
352 >      simError();
353 >    }
354  
355 < DumpReader::~DumpReader() {
355 >    int index = tokenizer.nextTokenAsInt();
356 >
357 >    StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
358  
359 < #ifdef IS_MPI
359 >    if (sd == NULL) {
360 >      return;
361 >    }
362 >    std::string type = tokenizer.nextToken();
363 >    int size = type.size();
364  
365 <    if (worldRank == 0) {
366 < #endif
367 <
368 <        int error;
369 <        error = fclose(inFile_);
370 <
371 <        if (error) {
372 <            sprintf(painCave.errMsg, "DumpReader Error: Error closing %s\n", filename_.c_str());
373 <            painCave.isFatal = 1;            
374 <            simError();
365 >    size_t found;
366 >    
367 >    if (needPos_) {
368 >      found = type.find("p");      
369 >      if (found == std::string::npos) {
370 >        sprintf(painCave.errMsg,
371 >                "DumpReader Error: StuntDouble %d has no Position\n"
372 >                "\tField (\"p\") specified.\n%s\n", index,
373 >                line.c_str());  
374 >        painCave.isFatal = 1;
375 >        simError();
376 >      }
377 >    }
378 >    
379 >    if (sd->isDirectional()) {
380 >      if (needQuaternion_) {
381 >        found = type.find("q");      
382 >        if (found == std::string::npos) {
383 >          sprintf(painCave.errMsg,
384 >                  "DumpReader Error: Directional StuntDouble %d has no\n"
385 >                  "\tQuaternion Field (\"q\") specified.\n%s\n", index,
386 >                  line.c_str());  
387 >          painCave.isFatal = 1;
388 >          simError();
389          }
390 <
116 <        MemoryUtils::deletePointers(framePos_);
117 <
118 < #ifdef IS_MPI
119 <
390 >      }      
391      }
392  
393 <    strcpy(checkPointMsg, "Dump file closed successfully.");
394 <    MPIcheckPoint();
393 >    for(int i = 0; i < size; ++i) {
394 >      switch(type[i]) {
395 >        
396 >        case 'p': {
397 >            Vector3d pos;
398 >            pos[0] = tokenizer.nextTokenAsDouble();
399 >            pos[1] = tokenizer.nextTokenAsDouble();
400 >            pos[2] = tokenizer.nextTokenAsDouble();
401 >            if (needPos_) {
402 >              sd->setPos(pos);
403 >            }            
404 >            break;
405 >        }
406 >        case 'v' : {
407 >            Vector3d vel;
408 >            vel[0] = tokenizer.nextTokenAsDouble();
409 >            vel[1] = tokenizer.nextTokenAsDouble();
410 >            vel[2] = tokenizer.nextTokenAsDouble();
411 >            if (needVel_) {
412 >              sd->setVel(vel);
413 >            }
414 >            break;
415 >        }
416  
417 < #endif
417 >        case 'q' : {
418 >           Quat4d q;
419 >           if (sd->isDirectional()) {
420 >              
421 >             q[0] = tokenizer.nextTokenAsDouble();
422 >             q[1] = tokenizer.nextTokenAsDouble();
423 >             q[2] = tokenizer.nextTokenAsDouble();
424 >             q[3] = tokenizer.nextTokenAsDouble();
425 >              
426 >             RealType qlen = q.length();
427 >             if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
428 >                
429 >               sprintf(painCave.errMsg,
430 >                       "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
431 >               painCave.isFatal = 1;
432 >               simError();
433 >                
434 >             }  
435 >              
436 >             q.normalize();
437 >             if (needQuaternion_) {            
438 >               sd->setQ(q);
439 >             }              
440 >           }            
441 >           break;
442 >        }  
443 >        case 'j' : {
444 >          Vector3d ji;
445 >          if (sd->isDirectional()) {
446 >             ji[0] = tokenizer.nextTokenAsDouble();
447 >             ji[1] = tokenizer.nextTokenAsDouble();
448 >             ji[2] = tokenizer.nextTokenAsDouble();
449 >             if (needAngMom_) {
450 >               sd->setJ(ji);
451 >             }
452 >          }
453 >          break;
454 >        }  
455 >        case 'f': {
456  
457 <    return;
458 < }
457 >          Vector3d force;
458 >          force[0] = tokenizer.nextTokenAsDouble();
459 >          force[1] = tokenizer.nextTokenAsDouble();
460 >          force[2] = tokenizer.nextTokenAsDouble();          
461 >          sd->setFrc(force);
462 >          break;
463 >        }
464 >        case 't' : {
465  
466 < int DumpReader::getNFrames(void) {
466 >           Vector3d torque;
467 >           torque[0] = tokenizer.nextTokenAsDouble();
468 >           torque[1] = tokenizer.nextTokenAsDouble();
469 >           torque[2] = tokenizer.nextTokenAsDouble();          
470 >           sd->setTrq(torque);          
471 >           break;
472 >        }
473 >        case 'u' : {
474  
475 <    if (!isScanned_)
476 <        scanFile();
475 >           RealType particlePot;
476 >           particlePot = tokenizer.nextTokenAsDouble();
477 >           sd->setParticlePot(particlePot);          
478 >           break;
479 >        }
480 >        case 'c' : {
481  
482 <    return nframes_;
483 < }
482 >           RealType flucQPos;
483 >           flucQPos = tokenizer.nextTokenAsDouble();
484 >           sd->setFlucQPos(flucQPos);          
485 >           break;
486 >        }
487 >        case 'w' : {
488  
489 < void DumpReader::scanFile(void) {
490 <  int i, j;
491 <  int lineNum = 0;
492 <  char readBuffer[maxBufferSize];
493 <  fpos_t * currPos;
489 >           RealType flucQVel;
490 >           flucQVel = tokenizer.nextTokenAsDouble();
491 >           sd->setFlucQVel(flucQVel);          
492 >           break;
493 >        }
494 >        case 'g' : {
495  
496 < #ifdef IS_MPI
497 <
498 <    if (worldRank == 0) {
499 < #endif // is_mpi
148 <
149 <        rewind(inFile_);
150 <
151 <        currPos = new fpos_t;
152 <        fgetpos(inFile_, currPos);
153 <        fgets(readBuffer, sizeof(readBuffer), inFile_);
154 <        lineNum++;
155 <
156 <        if (feof(inFile_)) {
157 <            sprintf(painCave.errMsg,
158 <                    "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
159 <                    filename_.c_str(),
160 <                    lineNum);
161 <            painCave.isFatal = 1;
162 <            simError();
496 >           RealType flucQFrc;
497 >           flucQFrc = tokenizer.nextTokenAsDouble();
498 >           sd->setFlucQFrc(flucQFrc);          
499 >           break;
500          }
501 +        case 'e' : {
502  
503 <        while (!feof(inFile_)) {
504 <            framePos_.push_back(currPos);
505 <
506 <            i = atoi(readBuffer);
507 <
508 <            fgets(readBuffer, sizeof(readBuffer), inFile_);
171 <            lineNum++;
172 <
173 <            if (feof(inFile_)) {
174 <                sprintf(painCave.errMsg,
175 <                        "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
176 <                        filename_.c_str(),
177 <                        lineNum);
178 <                painCave.isFatal = 1;
179 <                simError();
180 <            }
181 <
182 <            for(j = 0; j < i; j++) {
183 <                fgets(readBuffer, sizeof(readBuffer), inFile_);
184 <                lineNum++;
185 <
186 <                if (feof(inFile_)) {
187 <                    sprintf(painCave.errMsg,
188 <                            "DumpReader Error: File \"%s\" ended unexpectedly at line %d,"
189 <                                " with atom %d\n", filename_.c_str(),
190 <                            lineNum,
191 <                            j);
192 <
193 <                    painCave.isFatal = 1;
194 <                    simError();
195 <                }
196 <            }
197 <
198 <            currPos = new fpos_t;
199 <            fgetpos(inFile_, currPos);
200 <            fgets(readBuffer, sizeof(readBuffer), inFile_);
201 <            lineNum++;
503 >           Vector3d eField;
504 >           eField[0] = tokenizer.nextTokenAsDouble();
505 >           eField[1] = tokenizer.nextTokenAsDouble();
506 >           eField[2] = tokenizer.nextTokenAsDouble();          
507 >           sd->setElectricField(eField);          
508 >           break;
509          }
510 +        default: {
511 +               sprintf(painCave.errMsg,
512 +                       "DumpReader Error: %s is an unrecognized type\n", type.c_str());
513 +               painCave.isFatal = 1;
514 +               simError();
515 +          break;  
516 +        }
517  
518 <        delete currPos;
205 <        rewind(inFile_);
206 <        
207 <        nframes_ = framePos_.size();
208 < #ifdef IS_MPI
518 >      }
519      }
520 +    
521 +  }
522 +  
523  
524 <    MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
524 >  void DumpReader::parseSiteLine(const std::string& line) {
525  
526 <    strcpy(checkPointMsg, "Successfully scanned DumpFile\n");
527 <    MPIcheckPoint();
528 <
529 < #endif // is_mpi
530 <
531 <    isScanned_ = true;
532 < }
533 <
534 < void DumpReader::readFrame(int whichFrame) {
535 <    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
526 >    StringTokenizer tokenizer(line);
527 >    int nTokens;
528 >    
529 >    nTokens = tokenizer.countTokens();
530 >    
531 >    if (nTokens < 2) {  
532 >      sprintf(painCave.errMsg,
533 >              "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
534 >      painCave.isFatal = 1;
535 >      simError();
536 >    }
537  
538 <    if (storageLayout & DataStorage::dslPosition) {
539 <        needPos_ = true;
540 <    } else {
227 <        needPos_ = false;
228 <    }
538 >    /**
539 >     * The first token is the global integrable object index.
540 >     */
541  
542 <    if (storageLayout & DataStorage::dslVelocity) {
543 <        needVel_ = true;
544 <    } else {
545 <        needVel_ = false;
542 >    int index = tokenizer.nextTokenAsInt();
543 >    StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
544 >    if (sd == NULL) {
545 >      return;
546      }
547  
548 <    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
549 <        needQuaternion_ = true;
550 <    } else {
551 <        needQuaternion_ = false;
548 >    /**
549 >     * Test to see if the next token is an integer or not.  If not,
550 >     * we've got data on the integrable object itself.  If there is an
551 >     * integer, we're parsing data for a site on a rigid body.
552 >     */
553 >
554 >    std::string indexTest = tokenizer.peekNextToken();
555 >    std::istringstream i(indexTest);
556 >    int siteIndex;
557 >    if (i >> siteIndex) {
558 >      // chew up this token and parse as an int:
559 >      siteIndex = tokenizer.nextTokenAsInt();
560 >      RigidBody* rb = static_cast<RigidBody*>(sd);
561 >      sd = rb->getAtoms()[siteIndex];
562      }
563  
564 <    if (storageLayout & DataStorage::dslAngularMomentum) {
565 <        needAngMom_ = true;
566 <    } else {
567 <        needAngMom_ = false;    
564 >    /**
565 >     * The next token contains information on what follows.
566 >     */
567 >    std::string type = tokenizer.nextToken();
568 >    int size = type.size();
569 >    
570 >    for(int i = 0; i < size; ++i) {
571 >      switch(type[i]) {
572 >        
573 >      case 'u' : {
574 >        
575 >        RealType particlePot;
576 >        particlePot = tokenizer.nextTokenAsDouble();
577 >        sd->setParticlePot(particlePot);
578 >        break;
579 >      }
580 >      case 'c' : {
581 >        
582 >        RealType flucQPos;
583 >        flucQPos = tokenizer.nextTokenAsDouble();
584 >        sd->setFlucQPos(flucQPos);
585 >        break;
586 >      }
587 >      case 'w' : {
588 >        
589 >        RealType flucQVel;
590 >        flucQVel = tokenizer.nextTokenAsDouble();
591 >        sd->setFlucQVel(flucQVel);
592 >        break;
593 >      }
594 >      case 'g' : {
595 >        
596 >        RealType flucQFrc;
597 >        flucQFrc = tokenizer.nextTokenAsDouble();
598 >        sd->setFlucQFrc(flucQFrc);
599 >        break;
600 >      }
601 >      case 'e' : {
602 >        
603 >        Vector3d eField;
604 >        eField[0] = tokenizer.nextTokenAsDouble();
605 >        eField[1] = tokenizer.nextTokenAsDouble();
606 >        eField[2] = tokenizer.nextTokenAsDouble();  
607 >        sd->setElectricField(eField);          
608 >        break;
609 >      }
610 >      default: {
611 >        sprintf(painCave.errMsg,
612 >                "DumpReader Error: %s is an unrecognized type\n", type.c_str());
613 >        painCave.isFatal = 1;
614 >        simError();
615 >        break;  
616 >      }
617 >      }
618 >    }    
619 >  }
620 >  
621 >  
622 >  void  DumpReader::readStuntDoubles(std::istream& inputStream) {
623 >    
624 >    inputStream.getline(buffer, bufferSize);
625 >    std::string line(buffer);
626 >    
627 >    if (line.find("<StuntDoubles>") == std::string::npos) {
628 >      sprintf(painCave.errMsg,
629 >              "DumpReader Error: Missing <StuntDoubles>\n");
630 >      painCave.isFatal = 1;
631 >      simError();
632      }
633  
634 <    readSet(whichFrame);
635 < }
634 >    while(inputStream.getline(buffer, bufferSize)) {
635 >      line = buffer;
636 >      
637 >      if(line.find("</StuntDoubles>") != std::string::npos) {
638 >        break;
639 >      }
640  
641 < void DumpReader::readSet(int whichFrame) {
252 <  int i;
253 <  int nTotObjs;                  // the number of atoms
254 <  char read_buffer[maxBufferSize];  //the line buffer for reading
255 <  char * eof_test;               // ptr to see when we reach the end of the file
256 <
257 <  Molecule* mol;
258 <  StuntDouble* integrableObject;
259 <  SimInfo::MoleculeIterator mi;
260 <  Molecule::IntegrableObjectIterator ii;
261 <
262 < #ifndef IS_MPI
263 <
264 <    fsetpos(inFile_, framePos_[whichFrame]);
265 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
266 <
267 <    if (eof_test == NULL) {
268 <        sprintf(painCave.errMsg,
269 <                "DumpReader error: error reading 1st line of \"%s\"\n",
270 <                filename_.c_str());
271 <        painCave.isFatal = 1;
272 <        simError();
641 >      parseDumpLine(line);
642      }
643 +  
644 +  }
645  
646 <    nTotObjs = atoi(read_buffer);
646 >  void  DumpReader::readSiteData(std::istream& inputStream) {
647  
648 <    if (nTotObjs != info_->getNGlobalIntegrableObjects()) {
649 <        sprintf(painCave.errMsg,
650 <                "DumpReader error. %s nIntegrable, %d, "
651 <                    "does not match the meta-data file's nIntegrable, %d.\n",
652 <                filename_.c_str(),
653 <                nTotObjs,
283 <                info_->getNGlobalIntegrableObjects());
284 <
285 <        painCave.isFatal = 1;
286 <        simError();
648 >    inputStream.getline(buffer, bufferSize);
649 >    std::string line(buffer);
650 >    
651 >    if (line.find("<SiteData>") == std::string::npos) {
652 >      // site data isn't required for a simulation, so skip
653 >      return;
654      }
655  
656 <    //read the box mat from the comment line
656 >    while(inputStream.getline(buffer, bufferSize)) {
657 >      line = buffer;
658 >      
659 >      if(line.find("</SiteData>") != std::string::npos) {
660 >        break;
661 >      }
662  
663 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
292 <
293 <    if (eof_test == NULL) {
294 <        sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
295 <                filename_.c_str());
296 <        painCave.isFatal = 1;
297 <        simError();
663 >      parseSiteLine(line);
664      }
665 +  
666 +  }
667  
668 <    parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
668 >  void DumpReader::readFrameProperties(std::istream& inputStream) {
669  
670 <    //parse dump lines
670 >    Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
671 >    inputStream.getline(buffer, bufferSize);
672 >    std::string line(buffer);
673  
674 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
675 <
676 <        for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
677 <            integrableObject = mol->nextIntegrableObject(ii)) {          
678 <
309 <            eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
310 <
311 <            if (eof_test == NULL) {
312 <                sprintf(painCave.errMsg,
313 <                        "DumpReader Error: error in reading file %s\n"
314 <                            "natoms  = %d; index = %d\n"
315 <                            "error reading the line from the file.\n",
316 <                        filename_.c_str(),
317 <                        nTotObjs,
318 <                        i);
319 <
320 <                painCave.isFatal = 1;
321 <                simError();
322 <            }
323 <
324 <            parseDumpLine(read_buffer, integrableObject);
325 <            
326 <            }
674 >    if (line.find("<FrameData>") == std::string::npos) {
675 >      sprintf(painCave.errMsg,
676 >              "DumpReader Error: Missing <FrameData>\n");
677 >      painCave.isFatal = 1;
678 >      simError();
679      }
680  
681 <    // MPI Section of code..........
681 >    while(inputStream.getline(buffer, bufferSize)) {
682 >      line = buffer;
683 >      
684 >      if(line.find("</FrameData>") != std::string::npos) {
685 >        break;
686 >      }
687 >      
688 >      StringTokenizer tokenizer(line, " ;\t\n\r{}:,");
689 >      if (!tokenizer.hasMoreTokens()) {
690 >        sprintf(painCave.errMsg,
691 >                "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
692 >        painCave.isFatal = 1;
693 >        simError();      
694 >      }
695  
696 < #else //IS_MPI
697 <
698 <    // first thing first, suspend fatalities.
699 <    int masterNode = 0;
700 <    int nCurObj;
701 <    painCave.isEventLoop = 1;
702 <
703 <    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
704 <    int haveError;
705 <
706 <    MPI_Status istatus;
707 <    int nitems;
708 <
709 <    nTotObjs = info_->getNGlobalIntegrableObjects();
710 <    haveError = 0;
711 <
712 <    if (worldRank == masterNode) {
713 <        fsetpos(inFile_, framePos_[whichFrame]);
714 <
715 <        eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
696 >      std::string propertyName = tokenizer.nextToken();
697 >      if (propertyName == "Time") {
698 >        RealType currTime = tokenizer.nextTokenAsDouble();
699 >        s->setTime(currTime);
700 >      } else if (propertyName == "Hmat"){
701 >        Mat3x3d hmat;
702 >        hmat(0, 0) = tokenizer.nextTokenAsDouble();
703 >        hmat(0, 1) = tokenizer.nextTokenAsDouble();
704 >        hmat(0, 2) = tokenizer.nextTokenAsDouble();
705 >        hmat(1, 0) = tokenizer.nextTokenAsDouble();
706 >        hmat(1, 1) = tokenizer.nextTokenAsDouble();
707 >        hmat(1, 2) = tokenizer.nextTokenAsDouble();
708 >        hmat(2, 0) = tokenizer.nextTokenAsDouble();
709 >        hmat(2, 1) = tokenizer.nextTokenAsDouble();
710 >        hmat(2, 2) = tokenizer.nextTokenAsDouble();
711 >        s->setHmat(hmat);      
712 >      } else if (propertyName == "Thermostat") {
713 >        pair<RealType, RealType> thermostat;
714 >        thermostat.first = tokenizer.nextTokenAsDouble();
715 >        thermostat.second = tokenizer.nextTokenAsDouble();
716 >        s->setThermostat(thermostat);
717 >     } else if (propertyName == "Barostat") {
718 >        Mat3x3d eta;
719 >        eta(0, 0) = tokenizer.nextTokenAsDouble();
720 >        eta(0, 1) = tokenizer.nextTokenAsDouble();
721 >        eta(0, 2) = tokenizer.nextTokenAsDouble();
722 >        eta(1, 0) = tokenizer.nextTokenAsDouble();
723 >        eta(1, 1) = tokenizer.nextTokenAsDouble();
724 >        eta(1, 2) = tokenizer.nextTokenAsDouble();
725 >        eta(2, 0) = tokenizer.nextTokenAsDouble();
726 >        eta(2, 1) = tokenizer.nextTokenAsDouble();
727 >        eta(2, 2) = tokenizer.nextTokenAsDouble();
728 >        s->setBarostat(eta);
729 >      } else {
730 >        sprintf(painCave.errMsg,
731 >                "DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str());
732 >        painCave.isFatal = 0;
733 >        simError();        
734 >      }
735 >      
736 >    }
737  
738 <        if (eof_test == NULL) {
353 <            sprintf(painCave.errMsg, "DumpReader Error: Error reading 1st line of %s \n ",
354 <                    filename_.c_str());
355 <            painCave.isFatal = 1;
356 <            simError();
357 <        }
738 >  }
739  
740 <        nitems = atoi(read_buffer);
741 <
361 <        // Check to see that the number of integrable objects in the
362 <        // intial configuration file is the same as derived from the
363 <        // meta-data file.
364 <
365 <        if (nTotObjs != nitems) {
366 <            sprintf(painCave.errMsg,
367 <                    "DumpReader Error. %s nIntegrable, %d, "
368 <                        "does not match the meta-data file's nIntegrable, %d.\n",
369 <                    filename_.c_str(),
370 <                    nTotObjs,
371 <                    info_->getNGlobalIntegrableObjects());
372 <
373 <            painCave.isFatal = 1;
374 <            simError();
375 <        }
376 <
377 <        //read the boxMat from the comment line
378 <
379 <        eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
380 <
381 <        if (eof_test == NULL) {
382 <            sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
383 <                    filename_.c_str());
384 <            painCave.isFatal = 1;
385 <            simError();
386 <        }
387 <
388 <        //Every single processor will parse the comment line by itself
389 <        //By using this way, we might lose some efficiency, but if we want to add
390 <        //more parameters into comment line, we only need to modify function
391 <        //parseCommentLine
392 <
393 <        MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
394 <        parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
395 <
396 <        for(i = 0; i < info_->getNGlobalMolecules(); i++) {
397 <            int which_node = info_->getMolToProc(i);
398 <
399 <            if (which_node == masterNode) {
400 <                //molecules belong to master node
401 <
402 <                mol = info_->getMoleculeByGlobalIndex(i);
403 <
404 <                if (mol == NULL) {
405 <                    sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
406 <                        painCave.isFatal = 1;
407 <                    simError();
408 <                }
409 <
410 <                for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
411 <                       integrableObject = mol->nextIntegrableObject(ii)){
412 <                        
413 <                    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
414 <
415 <                    if (eof_test == NULL) {
416 <                        sprintf(painCave.errMsg,
417 <                                "DumpReader Error: error in reading file %s\n"
418 <                                    "natoms  = %d; index = %d\n"
419 <                                    "error reading the line from the file.\n",
420 <                                filename_.c_str(),
421 <                                nTotObjs,
422 <                                i);
423 <
424 <                        painCave.isFatal = 1;
425 <                        simError();
426 <                    }
427 <
428 <                    parseDumpLine(read_buffer, integrableObject);
429 <                }
430 <            } else {
431 <                //molecule belongs to slave nodes
432 <
433 <                MPI_Recv(&nCurObj, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
434 <                         MPI_COMM_WORLD, &istatus);
435 <
436 <                for(int j = 0; j < nCurObj; j++) {
437 <                    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
438 <
439 <                    if (eof_test == NULL) {
440 <                        sprintf(painCave.errMsg,
441 <                                "DumpReader Error: error in reading file %s\n"
442 <                                    "natoms  = %d; index = %d\n"
443 <                                    "error reading the line from the file.\n",
444 <                                filename_.c_str(),
445 <                                nTotObjs,
446 <                                i);
447 <
448 <                        painCave.isFatal = 1;
449 <                        simError();
450 <                    }
451 <                    
452 <                    MPI_Send(read_buffer, maxBufferSize, MPI_CHAR, which_node,
453 <                             TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
454 <                }
455 <            }
456 <        }
457 <    } else {
458 <        //actions taken at slave nodes
459 <        MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
460 <
461 <        /**@todo*/
462 <        parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
463 <
464 <        for(i = 0; i < info_->getNGlobalMolecules(); i++) {
465 <            int which_node = info_->getMolToProc(i);
466 <
467 <            if (which_node == worldRank) {
468 <                //molecule with global index i belongs to this processor
469 <                
470 <                mol = info_->getMoleculeByGlobalIndex(i);
471 <                if (mol == NULL) {
472 <                    sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
473 <                    painCave.isFatal = 1;
474 <                    simError();
475 <                }
476 <                
477 <                nCurObj = mol->getNIntegrableObjects();
478 <
479 <                MPI_Send(&nCurObj, 1, MPI_INT, masterNode, TAKE_THIS_TAG_INT,
480 <                         MPI_COMM_WORLD);
481 <
482 <                for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
483 <                       integrableObject = mol->nextIntegrableObject(ii)){
484 <                        
485 <                    MPI_Recv(read_buffer, maxBufferSize, MPI_CHAR, masterNode,
486 <                             TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
487 <
488 <                    parseDumpLine(read_buffer, integrableObject);
489 <                }
490 <                      
491 <            }
492 <            
493 <        }
494 <        
495 <    }
496 <
497 < #endif
498 <
499 < }
500 <
501 < void DumpReader::parseDumpLine(char *line, StuntDouble *integrableObject) {
502 <
503 <    Vector3d pos;  // position place holders
504 <    Vector3d vel;  // velocity placeholders
505 <    Quat4d q;    // the quaternions
506 <    Vector3d ji;   // angular velocity placeholders;
507 <    StringTokenizer tokenizer(line);
508 <    int nTokens;
509 <    
510 <    nTokens = tokenizer.countTokens();
511 <
512 <    if (nTokens < 14) {
513 <            sprintf(painCave.errMsg,
514 <                    "DumpReader Error: Not enough Tokens.\n%s\n", line);
515 <            painCave.isFatal = 1;
516 <            simError();
517 <    }
518 <
519 <    std::string name = tokenizer.nextToken();
520 <
521 <    if (name != integrableObject->getType()) {
522 <
523 <            sprintf(painCave.errMsg,
524 <                    "DumpReader Error: Atom type [%s] in %s does not match Atom Type [%s] in .md file.\n",
525 <                    name.c_str(), filename_.c_str(), integrableObject->getType().c_str());
526 <            painCave.isFatal = 1;
527 <            simError();        
528 <    }
529 <
530 <    pos[0] = tokenizer.nextTokenAsDouble();
531 <    pos[1] = tokenizer.nextTokenAsDouble();
532 <    pos[2] = tokenizer.nextTokenAsDouble();
533 <    if (needPos_) {
534 <        integrableObject->setPos(pos);
535 <    }
536 <    
537 <    vel[0] = tokenizer.nextTokenAsDouble();
538 <    vel[1] = tokenizer.nextTokenAsDouble();
539 <    vel[2] = tokenizer.nextTokenAsDouble();
540 <    if (needVel_) {
541 <        integrableObject->setVel(vel);
542 <    }
543 <    
544 <    if (integrableObject->isDirectional()) {
545 <        
546 <        q[0] = tokenizer.nextTokenAsDouble();
547 <        q[1] = tokenizer.nextTokenAsDouble();
548 <        q[2] = tokenizer.nextTokenAsDouble();
549 <        q[3] = tokenizer.nextTokenAsDouble();
550 <
551 <        double qlen = q.length();
552 <        if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
553 <            
554 <            sprintf(painCave.errMsg,
555 <                    "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0).\n");
556 <            painCave.isFatal = 1;
557 <            simError();
558 <            
559 <        }
560 <
561 <        q.normalize();
562 <        if (needQuaternion_) {          
563 <            integrableObject->setQ(q);
564 <        }
565 <
566 <        ji[0] = tokenizer.nextTokenAsDouble();
567 <        ji[1] = tokenizer.nextTokenAsDouble();
568 <        ji[2] = tokenizer.nextTokenAsDouble();
569 <        if (needAngMom_) {
570 <            integrableObject->setJ(ji);
571 <        }
572 <    }
573 <
574 < }
575 <
576 <
577 < void DumpReader::parseCommentLine(char* line, Snapshot* s) {
578 <    double currTime;
579 <    Mat3x3d hmat;
580 <    double chi;
581 <    double integralOfChiDt;
582 <    Mat3x3d eta;
583 <
584 <    StringTokenizer tokenizer(line);
585 <    int nTokens;
586 <
587 <    nTokens = tokenizer.countTokens();
588 <
589 <    //comment line should at least contain 10 tokens: current time(1 token) and  h-matrix(9 tokens)
590 <    if (nTokens < 10) {
591 <            sprintf(painCave.errMsg,
592 <                    "DumpReader Error: Not enough tokens in comment line: %s", line);
593 <            painCave.isFatal = 1;
594 <            simError();  
595 <    }
596 <
597 <    //read current time
598 <    currTime = tokenizer.nextTokenAsDouble();
599 <    s->setTime(currTime);
600 <    
601 <    //read h-matrix
602 <    hmat(0, 0) = tokenizer.nextTokenAsDouble();
603 <    hmat(0, 1) = tokenizer.nextTokenAsDouble();
604 <    hmat(0, 2) = tokenizer.nextTokenAsDouble();
605 <    hmat(1, 0) = tokenizer.nextTokenAsDouble();
606 <    hmat(1, 1) = tokenizer.nextTokenAsDouble();
607 <    hmat(1, 2) = tokenizer.nextTokenAsDouble();
608 <    hmat(2, 0) = tokenizer.nextTokenAsDouble();
609 <    hmat(2, 1) = tokenizer.nextTokenAsDouble();
610 <    hmat(2, 2) = tokenizer.nextTokenAsDouble();
611 <    s->setHmat(hmat);
612 <    
613 <    //read chi and integralOfChidt, they should apprear in pair
614 <    if (tokenizer.countTokens() >= 2) {
615 <        chi = tokenizer.nextTokenAsDouble();
616 <        integralOfChiDt = tokenizer.nextTokenAsDouble();            
617 <
618 <        s->setChi(chi);
619 <        s->setIntegralOfChiDt(integralOfChiDt);
620 <    }
621 <    
622 <    //read eta (eta is 3x3 matrix)
623 <    if (tokenizer.countTokens() >= 9) {
624 <        eta(0, 0) = tokenizer.nextTokenAsDouble();
625 <        eta(0, 1) = tokenizer.nextTokenAsDouble();
626 <        eta(0, 2) = tokenizer.nextTokenAsDouble();
627 <        eta(1, 0) = tokenizer.nextTokenAsDouble();
628 <        eta(1, 1) = tokenizer.nextTokenAsDouble();
629 <        eta(1, 2) = tokenizer.nextTokenAsDouble();
630 <        eta(2, 0) = tokenizer.nextTokenAsDouble();
631 <        eta(2, 1) = tokenizer.nextTokenAsDouble();
632 <        eta(2, 2) = tokenizer.nextTokenAsDouble();      
633 <
634 <        s->setEta(eta);
635 <    }
636 <
637 <    
638 < }
639 <
640 < }//end namespace oopse
740 >  
741 > }//end namespace OpenMD

Comparing:
trunk/src/io/DumpReader.cpp (property svn:keywords), Revision 398 by tim, Mon Mar 7 22:39:33 2005 UTC vs.
branches/development/src/io/DumpReader.cpp (property svn:keywords), Revision 1787 by gezelter, Wed Aug 29 18:13:11 2012 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines