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 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
branches/development/src/io/DumpReader.cpp (file contents), Revision 1764 by gezelter, Tue Jul 3 18:32:27 2012 UTC

# Line 1 | Line 1
1   /*
2 < * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
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   */
41
42 #define _LARGEFILE_SOURCE64
43 #define _FILE_OFFSET_BITS 64
44
45 #include <sys/types.h>
46 #include <sys/stat.h>
47
48 #include <iostream>
49 #include <math.h>
50
51 #include <stdio.h>
52 #include <stdlib.h>
53 #include <string.h>
54
55 #include "io/DumpReader.hpp"
56 #include "primitives/Molecule.hpp"
57 #include "utils/simError.h"
58 #include "utils/MemoryUtils.hpp"
59 #include "utils/StringTokenizer.hpp"
60
61 #ifdef IS_MPI
62
63 #include <mpi.h>
64 #define TAKE_THIS_TAG_CHAR 0
65 #define TAKE_THIS_TAG_INT 1
66
67 #endif // is_mpi
68
69
70 namespace oopse {
42    
43 <  DumpReader::DumpReader(SimInfo* info, const std::string& filename)
44 <    : info_(info), filename_(filename), isScanned_(false), nframes_(0) {
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
73 > #ifdef IS_MPI
74      
75 <      if (worldRank == 0) {
76 < #endif
75 >    if (worldRank == 0) {
76 > #endif
77        
78 <        inFile_ = fopen(filename_.c_str(), "r");
78 >      inFile_ = new std::ifstream(filename_.c_str());
79        
80 <        if (inFile_ == NULL) {
81 <          sprintf(painCave.errMsg, "DumpReader: Cannot open file: %s\n", filename_.c_str());
82 <          painCave.isFatal = 1;
83 <          simError();
84 <        }
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
88 > #ifdef IS_MPI
89        
90 <      }
90 >    }
91      
92 <      strcpy(checkPointMsg, "Dump file opened for reading successfully.");
93 <      MPIcheckPoint();
92 >    strcpy(checkPointMsg, "Dump file opened for reading successfully.");
93 >    errorCheckPoint();
94      
95 < #endif
95 > #endif
96      
97 <      return;
98 <    }
97 >    return;
98 >  }
99    
100 <  DumpReader::~DumpReader() {
100 >  DumpReader::~DumpReader() {
101      
102 < #ifdef IS_MPI
102 > #ifdef IS_MPI
103      
104 <    if (worldRank == 0) {
105 < #endif
104 >    if (worldRank == 0) {
105 > #endif
106        
107 <      int error;
108 <      error = fclose(inFile_);
107 >      delete inFile_;
108        
109 <      if (error) {
111 <        sprintf(painCave.errMsg, "DumpReader Error: Error closing %s\n", filename_.c_str());
112 <        painCave.isFatal = 1;            
113 <        simError();
114 <      }
109 > #ifdef IS_MPI
110        
111 <      MemoryUtils::deletePointers(framePos_);
117 <      
118 < #ifdef IS_MPI
119 <      
120 <    }
111 >    }
112      
113 <    strcpy(checkPointMsg, "Dump file closed successfully.");
114 <    MPIcheckPoint();
113 >    strcpy(checkPointMsg, "Dump file closed successfully.");
114 >    errorCheckPoint();
115      
116 < #endif
116 > #endif
117      
118 <    return;
119 <  }
118 >    return;
119 >  }
120    
121 <  int DumpReader::getNFrames(void) {
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 <    if (!isScanned_)
133 <      scanFile();
134 > #ifdef IS_MPI
135      
136 <    return nframes_;
137 <  }
137 <  
138 <  void DumpReader::scanFile(void) {
139 <    int i, j;
140 <    int lineNum = 0;
141 <    char readBuffer[maxBufferSize];
142 <    fpos_t * currPos;
143 <    
144 < #ifdef IS_MPI
145 <    
146 <    if (worldRank == 0) {
147 < #endif // is_mpi
136 >    if (worldRank == 0) {
137 > #endif // is_mpi
138        
139 <      rewind(inFile_);
140 <      
141 <      currPos = new fpos_t;
142 <      fgetpos(inFile_, currPos);
143 <      fgets(readBuffer, sizeof(readBuffer), inFile_);
144 <      lineNum++;
145 <      
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();
163 <      }
164 <      
165 <      while (!feof(inFile_)) {
166 <        framePos_.push_back(currPos);
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 <        i = atoi(readBuffer);
148 <        
149 <        fgets(readBuffer, sizeof(readBuffer), inFile_);
150 <        lineNum++;
151 <        
152 <        if (feof(inFile_)) {
153 <          sprintf(painCave.errMsg,
154 <                  "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
155 <                  filename_.c_str(),
156 <                  lineNum);
157 <          painCave.isFatal = 1;
158 <          simError();
159 <        }
181 <        
182 <        for(j = 0; j < i; j++) {
183 <          fgets(readBuffer, sizeof(readBuffer), inFile_);
184 <          lineNum++;
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 <          if (feof(inFile_)) {
162 <            sprintf(painCave.errMsg,
163 <                    "DumpReader Error: File \"%s\" ended unexpectedly at line %d,"
164 <                    " with atom %d\n", filename_.c_str(),
165 <                    lineNum,
166 <                    j);
167 <            
193 <            painCave.isFatal = 1;
194 <            simError();
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 <        
198 <        currPos = new fpos_t;
199 <        fgetpos(inFile_, currPos);
200 <        fgets(readBuffer, sizeof(readBuffer), inFile_);
201 <        lineNum++;
180 >        prevPos = currPos;
181        }
182        
183 <      delete currPos;
184 <      rewind(inFile_);
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 < #ifdef IS_MPI
195 <    }
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 <    MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
206 > #endif // is_mpi
207      
208 <    strcpy(checkPointMsg, "Successfully scanned DumpFile\n");
209 <    MPIcheckPoint();
210 <    
211 < #endif // is_mpi
212 <    
213 <    isScanned_ = true;
214 <  }
215 <  
216 <  void DumpReader::readFrame(int whichFrame) {
217 <    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
218 <    
219 <    if (storageLayout & DataStorage::dslPosition) {
220 <      needPos_ = true;
221 <    } else {
222 <      needPos_ = false;
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 <    if (storageLayout & DataStorage::dslVelocity) {
259 <      needVel_ = true;
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 <      needVel_ = false;
294 <    }
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 <    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
317 <      needQuaternion_ = true;
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 <      needQuaternion_ = false;
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 <    if (storageLayout & DataStorage::dslAngularMomentum) {
338 <      needAngMom_ = true;
339 <    } else {
340 <      needAngMom_ = false;    
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 <    readSet(whichFrame);
366 <  }
367 <  
368 <  void DumpReader::readSet(int whichFrame) {
369 <    int i;
370 <    int nTotObjs;                  // the number of atoms
371 <    char read_buffer[maxBufferSize];  //the line buffer for reading
372 <    char * eof_test;               // ptr to see when we reach the end of the file
373 <    
374 <    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();
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 <    nTotObjs = atoi(read_buffer);
378 <    
379 <    if (nTotObjs != info_->getNGlobalIntegrableObjects()) {
380 <      sprintf(painCave.errMsg,
381 <              "DumpReader error. %s nIntegrable, %d, "
382 <              "does not match the meta-data file's nIntegrable, %d.\n",
383 <              filename_.c_str(),
384 <              nTotObjs,
385 <              info_->getNGlobalIntegrableObjects());
386 <      
387 <      painCave.isFatal = 1;
388 <      simError();
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 <    //read the box mat from the comment line
392 <    
291 <    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();
298 <    }
299 <    
300 <    parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
301 <    
302 <    //parse dump lines
303 <    
304 <    i = 0;
305 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
306 <      
307 <      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
308 <           integrableObject = mol->nextIntegrableObject(ii)) {          
390 >
391 >    for(int i = 0; i < size; ++i) {
392 >      switch(type[i]) {
393          
394 <        eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
395 <        
396 <        if (eof_test == NULL) {
397 <          sprintf(painCave.errMsg,
398 <                  "DumpReader Error: error in reading file %s\n"
399 <                  "natoms  = %d; index = %d\n"
400 <                  "error reading the line from the file.\n",
401 <                  filename_.c_str(),
402 <                  nTotObjs,
319 <                  i);
320 <          
321 <          painCave.isFatal = 1;
322 <          simError();
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 <        
405 <        parseDumpLine(read_buffer, integrableObject);
406 <        i++;
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 <    // MPI Section of code..........
570 <    
571 < #else //IS_MPI
572 <    
573 <    // first thing first, suspend fatalities.
574 <    int masterNode = 0;
575 <    int nCurObj;
576 <    painCave.isEventLoop = 1;
577 <    
339 <    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
340 <    int haveError;
341 <    
342 <    MPI_Status istatus;
343 <    int nitems;
344 <    
345 <    nTotObjs = info_->getNGlobalIntegrableObjects();
346 <    haveError = 0;
347 <    
348 <    if (worldRank == masterNode) {
349 <      fsetpos(inFile_, framePos_[whichFrame]);
350 <      
351 <      eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
352 <      
353 <      if (eof_test == NULL) {
354 <        sprintf(painCave.errMsg, "DumpReader Error: Error reading 1st line of %s \n ",
355 <                filename_.c_str());
356 <        painCave.isFatal = 1;
357 <        simError();
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 <      
360 <      nitems = atoi(read_buffer);
361 <      
362 <      // Check to see that the number of integrable objects in the
363 <      // intial configuration file is the same as derived from the
364 <      // meta-data file.
365 <      
366 <      if (nTotObjs != nitems) {
367 <        sprintf(painCave.errMsg,
368 <                "DumpReader Error. %s nIntegrable, %d, "
369 <                "does not match the meta-data file's nIntegrable, %d.\n",
370 <                filename_.c_str(),
371 <                nTotObjs,
372 <                info_->getNGlobalIntegrableObjects());
579 >      case 'c' : {
580          
581 <        painCave.isFatal = 1;
582 <        simError();
581 >        RealType flucQPos;
582 >        flucQPos = tokenizer.nextTokenAsDouble();
583 >        sd->setFlucQPos(flucQPos);
584 >        break;
585        }
586 <      
587 <      //read the boxMat from the comment line
588 <      
589 <      eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
590 <      
591 <      if (eof_test == NULL) {
383 <        sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
384 <                filename_.c_str());
385 <        painCave.isFatal = 1;
386 <        simError();
586 >      case 'w' : {
587 >        
588 >        RealType flucQVel;
589 >        flucQVel = tokenizer.nextTokenAsDouble();
590 >        sd->setFlucQVel(flucQVel);
591 >        break;
592        }
593 <      
389 <      //Every single processor will parse the comment line by itself
390 <      //By using this way, we might lose some efficiency, but if we want to add
391 <      //more parameters into comment line, we only need to modify function
392 <      //parseCommentLine
393 <      
394 <      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
395 <      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
396 <      
397 <      for(i = 0; i < info_->getNGlobalMolecules(); i++) {
398 <        int which_node = info_->getMolToProc(i);
593 >      case 'g' : {
594          
595 <        if (which_node == masterNode) {
596 <          //molecules belong to master node
597 <          
598 <          mol = info_->getMoleculeByGlobalIndex(i);
404 <          
405 <          if (mol == NULL) {
406 <            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
407 <            painCave.isFatal = 1;
408 <            simError();
409 <          }
410 <          
411 <          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
412 <               integrableObject = mol->nextIntegrableObject(ii)){
413 <            
414 <            eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
415 <            
416 <            if (eof_test == NULL) {
417 <              sprintf(painCave.errMsg,
418 <                      "DumpReader Error: error in reading file %s\n"
419 <                      "natoms  = %d; index = %d\n"
420 <                      "error reading the line from the file.\n",
421 <                      filename_.c_str(),
422 <                      nTotObjs,
423 <                      i);
424 <              
425 <              painCave.isFatal = 1;
426 <              simError();
427 <            }
428 <            
429 <            parseDumpLine(read_buffer, integrableObject);
430 <          }
431 <        } else {
432 <          //molecule belongs to slave nodes
433 <          
434 <          MPI_Recv(&nCurObj, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
435 <                   MPI_COMM_WORLD, &istatus);
436 <          
437 <          for(int j = 0; j < nCurObj; j++) {
438 <            eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
439 <            
440 <            if (eof_test == NULL) {
441 <              sprintf(painCave.errMsg,
442 <                      "DumpReader Error: error in reading file %s\n"
443 <                      "natoms  = %d; index = %d\n"
444 <                      "error reading the line from the file.\n",
445 <                      filename_.c_str(),
446 <                      nTotObjs,
447 <                      i);
448 <              
449 <              painCave.isFatal = 1;
450 <              simError();
451 <            }
452 <            
453 <            MPI_Send(read_buffer, maxBufferSize, MPI_CHAR, which_node,
454 <                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
455 <          }
456 <        }
595 >        RealType flucQFrc;
596 >        flucQFrc = tokenizer.nextTokenAsDouble();
597 >        sd->setFlucQFrc(flucQFrc);
598 >        break;
599        }
600 <    } else {
459 <      //actions taken at slave nodes
460 <      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
461 <      
462 <      /**@todo*/
463 <      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
464 <      
465 <      for(i = 0; i < info_->getNGlobalMolecules(); i++) {
466 <        int which_node = info_->getMolToProc(i);
600 >      case 'e' : {
601          
602 <        if (which_node == worldRank) {
603 <          //molecule with global index i belongs to this processor
604 <          
605 <          mol = info_->getMoleculeByGlobalIndex(i);
606 <          if (mol == NULL) {
607 <            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
474 <            painCave.isFatal = 1;
475 <            simError();
476 <          }
477 <          
478 <          nCurObj = mol->getNIntegrableObjects();
479 <          
480 <          MPI_Send(&nCurObj, 1, MPI_INT, masterNode, TAKE_THIS_TAG_INT,
481 <                   MPI_COMM_WORLD);
482 <          
483 <          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
484 <               integrableObject = mol->nextIntegrableObject(ii)){
485 <            
486 <            MPI_Recv(read_buffer, maxBufferSize, MPI_CHAR, masterNode,
487 <                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
488 <            
489 <            parseDumpLine(read_buffer, integrableObject);
490 <          }
491 <          
492 <        }
493 <        
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 <      
610 <    }
611 <    
612 < #endif
613 <    
614 <  }
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 <  void DumpReader::parseDumpLine(char *line, StuntDouble *integrableObject) {
620 >  
621 >  void  DumpReader::readStuntDoubles(std::istream& inputStream) {
622      
623 <    Vector3d pos;  // position place holders
624 <    Vector3d vel;  // velocity placeholders
506 <    Quat4d q;    // the quaternions
507 <    Vector3d ji;   // angular velocity placeholders;
508 <    StringTokenizer tokenizer(line);
509 <    int nTokens;
623 >    inputStream.getline(buffer, bufferSize);
624 >    std::string line(buffer);
625      
626 <    nTokens = tokenizer.countTokens();
627 <    
628 <    if (nTokens < 14) {
629 <      sprintf(painCave.errMsg,
630 <              "DumpReader Error: Not enough Tokens.\n%s\n", line);
516 <      painCave.isFatal = 1;
517 <      simError();
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 <    std::string name = tokenizer.nextToken();
634 <    
522 <    if (name != integrableObject->getType()) {
632 >
633 >    while(inputStream.getline(buffer, bufferSize)) {
634 >      line = buffer;
635        
636 <      sprintf(painCave.errMsg,
637 <              "DumpReader Error: Atom type [%s] in %s does not match Atom Type [%s] in .md file.\n",
638 <              name.c_str(), filename_.c_str(), integrableObject->getType().c_str());
639 <      painCave.isFatal = 1;
640 <      simError();        
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 <    pos[0] = tokenizer.nextTokenAsDouble();
651 <    pos[1] = tokenizer.nextTokenAsDouble();
652 <    pos[2] = tokenizer.nextTokenAsDouble();
534 <    if (needPos_) {
535 <      integrableObject->setPos(pos);
650 >    if (line.find("<SiteData>") == std::string::npos) {
651 >      // site data isn't required for a simulation, so skip
652 >      return;
653      }
654 <    
655 <    vel[0] = tokenizer.nextTokenAsDouble();
656 <    vel[1] = tokenizer.nextTokenAsDouble();
540 <    vel[2] = tokenizer.nextTokenAsDouble();
541 <    if (needVel_) {
542 <      integrableObject->setVel(vel);
543 <    }
544 <    
545 <    if (integrableObject->isDirectional()) {
654 >
655 >    while(inputStream.getline(buffer, bufferSize)) {
656 >      line = buffer;
657        
658 <      q[0] = tokenizer.nextTokenAsDouble();
659 <      q[1] = tokenizer.nextTokenAsDouble();
549 <      q[2] = tokenizer.nextTokenAsDouble();
550 <      q[3] = tokenizer.nextTokenAsDouble();
551 <      
552 <      double qlen = q.length();
553 <      if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
554 <        
555 <        sprintf(painCave.errMsg,
556 <                "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0).\n");
557 <        painCave.isFatal = 1;
558 <        simError();
559 <        
560 <      }
561 <      
562 <      q.normalize();
563 <      if (needQuaternion_) {          
564 <        integrableObject->setQ(q);
658 >      if(line.find("</SiteData>") != std::string::npos) {
659 >        break;
660        }
661 <      
662 <      ji[0] = tokenizer.nextTokenAsDouble();
568 <      ji[1] = tokenizer.nextTokenAsDouble();
569 <      ji[2] = tokenizer.nextTokenAsDouble();
570 <      if (needAngMom_) {
571 <        integrableObject->setJ(ji);
572 <      }
661 >
662 >      parseSiteLine(line);
663      }
574    
575  }
664    
665 <  
666 <  void DumpReader::parseCommentLine(char* line, Snapshot* s) {
667 <    double currTime;
668 <    Mat3x3d hmat;
669 <    double chi;
670 <    double integralOfChiDt;
671 <    Mat3x3d eta;
672 <    
673 <    StringTokenizer tokenizer(line);
674 <    int nTokens;
675 <    
676 <    nTokens = tokenizer.countTokens();
677 <    
590 <    //comment line should at least contain 10 tokens: current time(1 token) and  h-matrix(9 tokens)
591 <    if (nTokens < 10) {
592 <      sprintf(painCave.errMsg,
593 <              "DumpReader Error: Not enough tokens in comment line: %s", line);
594 <      painCave.isFatal = 1;
595 <      simError();  
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 <    //read current time
681 <    currTime = tokenizer.nextTokenAsDouble();
600 <    s->setTime(currTime);
601 <    
602 <    //read h-matrix
603 <    hmat(0, 0) = tokenizer.nextTokenAsDouble();
604 <    hmat(0, 1) = tokenizer.nextTokenAsDouble();
605 <    hmat(0, 2) = tokenizer.nextTokenAsDouble();
606 <    hmat(1, 0) = tokenizer.nextTokenAsDouble();
607 <    hmat(1, 1) = tokenizer.nextTokenAsDouble();
608 <    hmat(1, 2) = tokenizer.nextTokenAsDouble();
609 <    hmat(2, 0) = tokenizer.nextTokenAsDouble();
610 <    hmat(2, 1) = tokenizer.nextTokenAsDouble();
611 <    hmat(2, 2) = tokenizer.nextTokenAsDouble();
612 <    s->setHmat(hmat);
613 <    
614 <    //read chi and integralOfChidt, they should apprear in pair
615 <    if (tokenizer.countTokens() >= 2) {
616 <      chi = tokenizer.nextTokenAsDouble();
617 <      integralOfChiDt = tokenizer.nextTokenAsDouble();            
679 >
680 >    while(inputStream.getline(buffer, bufferSize)) {
681 >      line = buffer;
682        
683 <      s->setChi(chi);
684 <      s->setIntegralOfChiDt(integralOfChiDt);
685 <    }
622 <    
623 <    //read eta (eta is 3x3 matrix)
624 <    if (tokenizer.countTokens() >= 9) {
625 <      eta(0, 0) = tokenizer.nextTokenAsDouble();
626 <      eta(0, 1) = tokenizer.nextTokenAsDouble();
627 <      eta(0, 2) = tokenizer.nextTokenAsDouble();
628 <      eta(1, 0) = tokenizer.nextTokenAsDouble();
629 <      eta(1, 1) = tokenizer.nextTokenAsDouble();
630 <      eta(1, 2) = tokenizer.nextTokenAsDouble();
631 <      eta(2, 0) = tokenizer.nextTokenAsDouble();
632 <      eta(2, 1) = tokenizer.nextTokenAsDouble();
633 <      eta(2, 2) = tokenizer.nextTokenAsDouble();      
683 >      if(line.find("</FrameData>") != std::string::npos) {
684 >        break;
685 >      }
686        
687 <      s->setEta(eta);
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 <    
638 <    
736 >
737    }
738 <  
739 < }//end namespace oopse
738 >
739 >  
740 > }//end namespace OpenMD

Comparing:
trunk/src/io/DumpReader.cpp (property svn:keywords), Revision 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
branches/development/src/io/DumpReader.cpp (property svn:keywords), Revision 1764 by gezelter, Tue Jul 3 18:32:27 2012 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines