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 601 by tim, Thu Sep 15 19:17:04 2005 UTC vs.
branches/development/src/io/DumpReader.cpp (file contents), Revision 1794 by gezelter, Thu Sep 6 19:44:06 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 >                                  ifstream::in | ifstream::binary);
80        
81 <        if (inFile_ == NULL) {
82 <          sprintf(painCave.errMsg, "DumpReader: Cannot open file: %s\n", filename_.c_str());
83 <          painCave.isFatal = 1;
84 <          simError();
85 <        }
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
89 > #ifdef IS_MPI
90        
91 <      }
91 >    }
92      
93 <      strcpy(checkPointMsg, "Dump file opened for reading successfully.");
94 <      MPIcheckPoint();
93 >    strcpy(checkPointMsg, "Dump file opened for reading successfully.");
94 >    errorCheckPoint();
95      
96 < #endif
96 > #endif
97      
98 <      return;
99 <    }
98 >    return;
99 >  }
100    
101 <  DumpReader::~DumpReader() {
101 >  DumpReader::~DumpReader() {
102      
103 < #ifdef IS_MPI
103 > #ifdef IS_MPI
104      
105 <    if (worldRank == 0) {
106 < #endif
105 >    if (worldRank == 0) {
106 > #endif
107        
108 <      int error;
108 <      error = fclose(inFile_);
108 >      delete inFile_;
109        
110 <      if (error) {
111 <        sprintf(painCave.errMsg, "DumpReader Error: Error closing %s\n", filename_.c_str());
112 <        painCave.isFatal = 1;            
113 <        simError();
114 <      }
110 > #ifdef IS_MPI
111        
112 <      MemoryUtils::deletePointers(framePos_);
117 <      
118 < #ifdef IS_MPI
119 <      
120 <    }
112 >    }
113      
114 <    strcpy(checkPointMsg, "Dump file closed successfully.");
115 <    MPIcheckPoint();
114 >    strcpy(checkPointMsg, "Dump file closed successfully.");
115 >    errorCheckPoint();
116      
117 < #endif
117 > #endif
118      
119 <    return;
120 <  }
119 >    return;
120 >  }
121    
122 <  int DumpReader::getNFrames(void) {
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 <    if (!isScanned_)
133 <      scanFile();
135 > #ifdef IS_MPI
136      
137 <    return nframes_;
138 <  }
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
137 >    if (worldRank == 0) {
138 > #endif // is_mpi
139        
140 <      rewind(inFile_);
141 <      
142 <      currPos = new fpos_t;
143 <      fgetpos(inFile_, currPos);
144 <      fgets(readBuffer, sizeof(readBuffer), inFile_);
145 <      lineNum++;
146 <      
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);
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 <        i = atoi(readBuffer);
149 <        
150 <        fgets(readBuffer, sizeof(readBuffer), inFile_);
151 <        lineNum++;
152 <        
153 <        if (feof(inFile_)) {
154 <          sprintf(painCave.errMsg,
155 <                  "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
156 <                  filename_.c_str(),
157 <                  lineNum);
158 <          painCave.isFatal = 1;
159 <          simError();
160 <        }
181 <        
182 <        for(j = 0; j < i; j++) {
183 <          fgets(readBuffer, sizeof(readBuffer), inFile_);
184 <          lineNum++;
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 <          if (feof(inFile_)) {
163 <            sprintf(painCave.errMsg,
164 <                    "DumpReader Error: File \"%s\" ended unexpectedly at line %d,"
165 <                    " with atom %d\n", filename_.c_str(),
166 <                    lineNum,
167 <                    j);
168 <            
193 <            painCave.isFatal = 1;
194 <            simError();
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 <        
198 <        currPos = new fpos_t;
199 <        fgetpos(inFile_, currPos);
200 <        fgets(readBuffer, sizeof(readBuffer), inFile_);
201 <        lineNum++;
181 >        prevPos = currPos;
182        }
183        
184 <      delete currPos;
185 <      rewind(inFile_);
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 < #ifdef IS_MPI
196 <    }
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 <    MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
207 > #endif // is_mpi
208      
209 <    strcpy(checkPointMsg, "Successfully scanned DumpFile\n");
210 <    MPIcheckPoint();
211 <    
212 < #endif // is_mpi
213 <    
214 <    isScanned_ = true;
215 <  }
216 <  
217 <  void DumpReader::readFrame(int whichFrame) {
218 <    if (!isScanned_)
219 <      scanFile();
220 <        
221 <    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
222 <    
223 <    if (storageLayout & DataStorage::dslPosition) {
224 <      needPos_ = true;
225 <    } else {
226 <      needPos_ = false;
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 <    if (storageLayout & DataStorage::dslVelocity) {
262 <      needVel_ = true;
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 <      needVel_ = false;
297 <    }
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 <    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
320 <      needQuaternion_ = true;
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 <      needQuaternion_ = false;
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 <    if (storageLayout & DataStorage::dslAngularMomentum) {
341 <      needAngMom_ = true;
342 <    } else {
343 <      needAngMom_ = false;    
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 <    readSet(whichFrame);
369 <  }
370 <  
371 <  void DumpReader::readSet(int whichFrame) {
372 <    int i;
373 <    int nTotObjs;                  // the number of atoms
374 <    char read_buffer[maxBufferSize];  //the line buffer for reading
375 <    char * eof_test;               // ptr to see when we reach the end of the file
376 <    
377 <    Molecule* mol;
261 <    StuntDouble* integrableObject;
262 <    SimInfo::MoleculeIterator mi;
263 <    Molecule::IntegrableObjectIterator ii;
264 <    
265 < #ifndef IS_MPI
266 <    
267 <    fsetpos(inFile_, framePos_[whichFrame]);
268 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
269 <    
270 <    if (eof_test == NULL) {
271 <      sprintf(painCave.errMsg,
272 <              "DumpReader error: error reading 1st line of \"%s\"\n",
273 <              filename_.c_str());
274 <      painCave.isFatal = 1;
275 <      simError();
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 <    nTotObjs = atoi(read_buffer);
381 <    
382 <    if (nTotObjs != info_->getNGlobalIntegrableObjects()) {
383 <      sprintf(painCave.errMsg,
384 <              "DumpReader error. %s nIntegrable, %d, "
385 <              "does not match the meta-data file's nIntegrable, %d.\n",
386 <              filename_.c_str(),
387 <              nTotObjs,
388 <              info_->getNGlobalIntegrableObjects());
389 <      
390 <      painCave.isFatal = 1;
391 <      simError();
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 <    //read the box mat from the comment line
395 <    
294 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
295 <    
296 <    if (eof_test == NULL) {
297 <      sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
298 <              filename_.c_str());
299 <      painCave.isFatal = 1;
300 <      simError();
301 <    }
302 <    
303 <    parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
304 <    
305 <    //parse dump lines
306 <    
307 <    i = 0;
308 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
309 <      
310 <      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
311 <           integrableObject = mol->nextIntegrableObject(ii)) {          
393 >
394 >    for(int i = 0; i < size; ++i) {
395 >      switch(type[i]) {
396          
397 <        eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
398 <        
399 <        if (eof_test == NULL) {
400 <          sprintf(painCave.errMsg,
401 <                  "DumpReader Error: error in reading file %s\n"
402 <                  "natoms  = %d; index = %d\n"
403 <                  "error reading the line from the file.\n",
404 <                  filename_.c_str(),
405 <                  nTotObjs,
322 <                  i);
323 <          
324 <          painCave.isFatal = 1;
325 <          simError();
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 <        
408 <        parseDumpLine(read_buffer, integrableObject);
409 <        i++;
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 <    // MPI Section of code..........
572 <    
573 < #else //IS_MPI
574 <    
575 <    // first thing first, suspend fatalities.
576 <    int masterNode = 0;
577 <    int nCurObj;
578 <    painCave.isEventLoop = 1;
579 <    
342 <    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
343 <    int haveError;
344 <    
345 <    MPI_Status istatus;
346 <    int nitems;
347 <    
348 <    nTotObjs = info_->getNGlobalIntegrableObjects();
349 <    haveError = 0;
350 <    
351 <    if (worldRank == masterNode) {
352 <      fsetpos(inFile_, framePos_[whichFrame]);
353 <      
354 <      eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
355 <      
356 <      if (eof_test == NULL) {
357 <        sprintf(painCave.errMsg, "DumpReader Error: Error reading 1st line of %s \n ",
358 <                filename_.c_str());
359 <        painCave.isFatal = 1;
360 <        simError();
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 <      
363 <      nitems = atoi(read_buffer);
364 <      
365 <      // Check to see that the number of integrable objects in the
366 <      // intial configuration file is the same as derived from the
367 <      // meta-data file.
368 <      
369 <      if (nTotObjs != nitems) {
370 <        sprintf(painCave.errMsg,
371 <                "DumpReader Error. %s nIntegrable, %d, "
372 <                "does not match the meta-data file's nIntegrable, %d.\n",
373 <                filename_.c_str(),
374 <                nTotObjs,
375 <                info_->getNGlobalIntegrableObjects());
581 >      case 'c' : {
582          
583 <        painCave.isFatal = 1;
584 <        simError();
583 >        RealType flucQPos;
584 >        flucQPos = tokenizer.nextTokenAsDouble();
585 >        sd->setFlucQPos(flucQPos);
586 >        break;
587        }
588 <      
589 <      //read the boxMat from the comment line
590 <      
591 <      eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
592 <      
593 <      if (eof_test == NULL) {
386 <        sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
387 <                filename_.c_str());
388 <        painCave.isFatal = 1;
389 <        simError();
588 >      case 'w' : {
589 >        
590 >        RealType flucQVel;
591 >        flucQVel = tokenizer.nextTokenAsDouble();
592 >        sd->setFlucQVel(flucQVel);
593 >        break;
594        }
595 <      
392 <      //Every single processor will parse the comment line by itself
393 <      //By using this way, we might lose some efficiency, but if we want to add
394 <      //more parameters into comment line, we only need to modify function
395 <      //parseCommentLine
396 <      
397 <      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
398 <      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
399 <      
400 <      for(i = 0; i < info_->getNGlobalMolecules(); i++) {
401 <        int which_node = info_->getMolToProc(i);
595 >      case 'g' : {
596          
597 <        if (which_node == masterNode) {
598 <          //molecules belong to master node
599 <          
600 <          mol = info_->getMoleculeByGlobalIndex(i);
407 <          
408 <          if (mol == NULL) {
409 <            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
410 <            painCave.isFatal = 1;
411 <            simError();
412 <          }
413 <          
414 <          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
415 <               integrableObject = mol->nextIntegrableObject(ii)){
416 <            
417 <            eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
418 <            
419 <            if (eof_test == NULL) {
420 <              sprintf(painCave.errMsg,
421 <                      "DumpReader Error: error in reading file %s\n"
422 <                      "natoms  = %d; index = %d\n"
423 <                      "error reading the line from the file.\n",
424 <                      filename_.c_str(),
425 <                      nTotObjs,
426 <                      i);
427 <              
428 <              painCave.isFatal = 1;
429 <              simError();
430 <            }
431 <            
432 <            parseDumpLine(read_buffer, integrableObject);
433 <          }
434 <        } else {
435 <          //molecule belongs to slave nodes
436 <          
437 <          MPI_Recv(&nCurObj, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
438 <                   MPI_COMM_WORLD, &istatus);
439 <          
440 <          for(int j = 0; j < nCurObj; j++) {
441 <            eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
442 <            
443 <            if (eof_test == NULL) {
444 <              sprintf(painCave.errMsg,
445 <                      "DumpReader Error: error in reading file %s\n"
446 <                      "natoms  = %d; index = %d\n"
447 <                      "error reading the line from the file.\n",
448 <                      filename_.c_str(),
449 <                      nTotObjs,
450 <                      i);
451 <              
452 <              painCave.isFatal = 1;
453 <              simError();
454 <            }
455 <            
456 <            MPI_Send(read_buffer, maxBufferSize, MPI_CHAR, which_node,
457 <                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
458 <          }
459 <        }
597 >        RealType flucQFrc;
598 >        flucQFrc = tokenizer.nextTokenAsDouble();
599 >        sd->setFlucQFrc(flucQFrc);
600 >        break;
601        }
602 <    } else {
462 <      //actions taken at slave nodes
463 <      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
464 <      
465 <      /**@todo*/
466 <      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
467 <      
468 <      for(i = 0; i < info_->getNGlobalMolecules(); i++) {
469 <        int which_node = info_->getMolToProc(i);
602 >      case 'e' : {
603          
604 <        if (which_node == worldRank) {
605 <          //molecule with global index i belongs to this processor
606 <          
607 <          mol = info_->getMoleculeByGlobalIndex(i);
608 <          if (mol == NULL) {
609 <            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
477 <            painCave.isFatal = 1;
478 <            simError();
479 <          }
480 <          
481 <          nCurObj = mol->getNIntegrableObjects();
482 <          
483 <          MPI_Send(&nCurObj, 1, MPI_INT, masterNode, TAKE_THIS_TAG_INT,
484 <                   MPI_COMM_WORLD);
485 <          
486 <          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
487 <               integrableObject = mol->nextIntegrableObject(ii)){
488 <            
489 <            MPI_Recv(read_buffer, maxBufferSize, MPI_CHAR, masterNode,
490 <                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
491 <            
492 <            parseDumpLine(read_buffer, integrableObject);
493 <          }
494 <          
495 <        }
496 <        
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 <      
612 <    }
613 <    
614 < #endif
615 <    
616 <  }
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 <  void DumpReader::parseDumpLine(char *line, StuntDouble *integrableObject) {
622 >  
623 >  void  DumpReader::readStuntDoubles(std::istream& inputStream) {
624      
625 <    Vector3d pos;  // position place holders
626 <    Vector3d vel;  // velocity placeholders
509 <    Quat4d q;    // the quaternions
510 <    Vector3d ji;   // angular velocity placeholders;
511 <    StringTokenizer tokenizer(line);
512 <    int nTokens;
625 >    inputStream.getline(buffer, bufferSize);
626 >    std::string line(buffer);
627      
628 <    nTokens = tokenizer.countTokens();
629 <    
630 <    if (nTokens < 14) {
631 <      sprintf(painCave.errMsg,
632 <              "DumpReader Error: Not enough Tokens.\n%s\n", line);
519 <      painCave.isFatal = 1;
520 <      simError();
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 <    std::string name = tokenizer.nextToken();
636 <    
525 <    if (name != integrableObject->getType()) {
634 >
635 >    while(inputStream.getline(buffer, bufferSize)) {
636 >      line = buffer;
637        
638 <      sprintf(painCave.errMsg,
639 <              "DumpReader Error: Atom type [%s] in %s does not match Atom Type [%s] in .md file.\n",
640 <              name.c_str(), filename_.c_str(), integrableObject->getType().c_str());
641 <      painCave.isFatal = 1;
642 <      simError();        
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 <    pos[0] = tokenizer.nextTokenAsDouble();
653 <    pos[1] = tokenizer.nextTokenAsDouble();
654 <    pos[2] = tokenizer.nextTokenAsDouble();
537 <    if (needPos_) {
538 <      integrableObject->setPos(pos);
652 >    if (line.find("<SiteData>") == std::string::npos) {
653 >      // site data isn't required for a simulation, so skip
654 >      return;
655      }
656 <    
657 <    vel[0] = tokenizer.nextTokenAsDouble();
658 <    vel[1] = tokenizer.nextTokenAsDouble();
543 <    vel[2] = tokenizer.nextTokenAsDouble();
544 <    if (needVel_) {
545 <      integrableObject->setVel(vel);
546 <    }
547 <    
548 <    if (integrableObject->isDirectional()) {
656 >
657 >    while(inputStream.getline(buffer, bufferSize)) {
658 >      line = buffer;
659        
660 <      q[0] = tokenizer.nextTokenAsDouble();
661 <      q[1] = tokenizer.nextTokenAsDouble();
552 <      q[2] = tokenizer.nextTokenAsDouble();
553 <      q[3] = tokenizer.nextTokenAsDouble();
554 <      
555 <      double qlen = q.length();
556 <      if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
557 <        
558 <        sprintf(painCave.errMsg,
559 <                "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0).\n");
560 <        painCave.isFatal = 1;
561 <        simError();
562 <        
563 <      }
564 <      
565 <      q.normalize();
566 <      if (needQuaternion_) {          
567 <        integrableObject->setQ(q);
660 >      if(line.find("</SiteData>") != std::string::npos) {
661 >        break;
662        }
663 <      
664 <      ji[0] = tokenizer.nextTokenAsDouble();
571 <      ji[1] = tokenizer.nextTokenAsDouble();
572 <      ji[2] = tokenizer.nextTokenAsDouble();
573 <      if (needAngMom_) {
574 <        integrableObject->setJ(ji);
575 <      }
663 >
664 >      parseSiteLine(line);
665      }
577    
578  }
666    
667 <  
668 <  void DumpReader::parseCommentLine(char* line, Snapshot* s) {
669 <    double currTime;
670 <    Mat3x3d hmat;
671 <    double chi;
672 <    double integralOfChiDt;
673 <    Mat3x3d eta;
674 <    
675 <    StringTokenizer tokenizer(line);
676 <    int nTokens;
677 <    
678 <    nTokens = tokenizer.countTokens();
679 <    
593 <    //comment line should at least contain 10 tokens: current time(1 token) and  h-matrix(9 tokens)
594 <    if (nTokens < 10) {
595 <      sprintf(painCave.errMsg,
596 <              "DumpReader Error: Not enough tokens in comment line: %s", line);
597 <      painCave.isFatal = 1;
598 <      simError();  
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 <    //read current time
683 <    currTime = tokenizer.nextTokenAsDouble();
603 <    s->setTime(currTime);
604 <    
605 <    //read h-matrix
606 <    hmat(0, 0) = tokenizer.nextTokenAsDouble();
607 <    hmat(0, 1) = tokenizer.nextTokenAsDouble();
608 <    hmat(0, 2) = tokenizer.nextTokenAsDouble();
609 <    hmat(1, 0) = tokenizer.nextTokenAsDouble();
610 <    hmat(1, 1) = tokenizer.nextTokenAsDouble();
611 <    hmat(1, 2) = tokenizer.nextTokenAsDouble();
612 <    hmat(2, 0) = tokenizer.nextTokenAsDouble();
613 <    hmat(2, 1) = tokenizer.nextTokenAsDouble();
614 <    hmat(2, 2) = tokenizer.nextTokenAsDouble();
615 <    s->setHmat(hmat);
616 <    
617 <    //read chi and integralOfChidt, they should apprear in pair
618 <    if (tokenizer.countTokens() >= 2) {
619 <      chi = tokenizer.nextTokenAsDouble();
620 <      integralOfChiDt = tokenizer.nextTokenAsDouble();            
681 >
682 >    while(inputStream.getline(buffer, bufferSize)) {
683 >      line = buffer;
684        
685 <      s->setChi(chi);
686 <      s->setIntegralOfChiDt(integralOfChiDt);
687 <    }
625 <    
626 <    //read eta (eta is 3x3 matrix)
627 <    if (tokenizer.countTokens() >= 9) {
628 <      eta(0, 0) = tokenizer.nextTokenAsDouble();
629 <      eta(0, 1) = tokenizer.nextTokenAsDouble();
630 <      eta(0, 2) = tokenizer.nextTokenAsDouble();
631 <      eta(1, 0) = tokenizer.nextTokenAsDouble();
632 <      eta(1, 1) = tokenizer.nextTokenAsDouble();
633 <      eta(1, 2) = tokenizer.nextTokenAsDouble();
634 <      eta(2, 0) = tokenizer.nextTokenAsDouble();
635 <      eta(2, 1) = tokenizer.nextTokenAsDouble();
636 <      eta(2, 2) = tokenizer.nextTokenAsDouble();      
685 >      if(line.find("</FrameData>") != std::string::npos) {
686 >        break;
687 >      }
688        
689 <      s->setEta(eta);
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 <    
641 <    
738 >
739    }
740 <  
741 < }//end namespace oopse
740 >
741 >  
742 > }//end namespace OpenMD

Comparing:
trunk/src/io/DumpReader.cpp (property svn:keywords), Revision 601 by tim, Thu Sep 15 19:17:04 2005 UTC vs.
branches/development/src/io/DumpReader.cpp (property svn:keywords), Revision 1794 by gezelter, Thu Sep 6 19:44:06 2012 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines