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 132 by tim, Thu Oct 21 16:22:01 2004 UTC vs.
branches/development/src/io/DumpReader.cpp (file contents), Revision 1714 by gezelter, Sat May 19 18:12:46 2012 UTC

# Line 1 | Line 1
1 < #define _LARGEFILE_SOURCE64
2 < #define _FILE_OFFSET_BITS 64
3 <
4 < #include <sys/types.h>
5 < #include <sys/stat.h>
6 <
7 < #include <iostream>
8 < #include <math.h>
9 <
10 < #include <stdio.h>
11 < #include <stdlib.h>
12 < #include <string.h>
13 <
14 <
15 < #include "io/ReadWrite.hpp"
16 < #include "utils/simError.h"
17 <
18 < #ifdef IS_MPI
19 < #include <mpi.h>
20 < #include "brains/mpiSimulation.hpp"
21 < #define TAKE_THIS_TAG_CHAR 0
22 < #define TAKE_THIS_TAG_INT 1
23 < #endif // is_mpi
24 <
25 <
26 < DumpReader :: DumpReader(const char *in_name ){
27 <
28 <  isScanned = false;
29 <
30 < #ifdef IS_MPI
31 <  if (worldRank == 0) {
32 < #endif
33 <
34 <  inFile = fopen(in_name, "r");
35 <  if(inFile == NULL){
36 <    sprintf(painCave.errMsg,
37 <            "Cannot open file: %s\n", in_name);
38 <    painCave.isFatal = 1;
39 <    simError();
40 <  }
1 > /*
2 > * Copyright (c) 2009 The University of Notre Dame. All Rights Reserved.
3 > *
4 > * The University of Notre Dame grants you ("Licensee") a
5 > * non-exclusive, royalty free, license to use, modify and
6 > * redistribute this software in source and binary code form, provided
7 > * that the following conditions are met:
8 > *
9 > * 1. Redistributions of source code must retain the above copyright
10 > *    notice, this list of conditions and the following disclaimer.
11 > *
12 > * 2. Redistributions in binary form must reproduce the above copyright
13 > *    notice, this list of conditions and the following disclaimer in the
14 > *    documentation and/or other materials provided with the
15 > *    distribution.
16 > *
17 > * This software is provided "AS IS," without a warranty of any
18 > * kind. All express or implied conditions, representations and
19 > * warranties, including any implied warranty of merchantability,
20 > * fitness for a particular purpose or non-infringement, are hereby
21 > * excluded.  The University of Notre Dame and its licensors shall not
22 > * be liable for any damages suffered by licensee as a result of
23 > * using, modifying or distributing the software or its
24 > * derivatives. In no event will the University of Notre Dame or its
25 > * licensors be liable for any lost revenue, profit or data, or for
26 > * direct, indirect, special, consequential, incidental or punitive
27 > * damages, however caused and regardless of the theory of liability,
28 > * arising out of the use of or inability to use software, even if the
29 > * University of Notre Dame has been advised of the possibility of
30 > * such damages.
31 > *
32 > * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 > * research, please cite the appropriate papers when you publish your
34 > * work.  Good starting points are:
35 > *                                                                      
36 > * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 > * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 > */
42    
43 <  inFileName = in_name;
44 < #ifdef IS_MPI
45 <  }
46 <  strcpy( checkPointMsg, "Dump file opened for reading successfully." );
47 <  MPIcheckPoint();
48 < #endif
49 <  return;  
50 < }
51 <
52 < DumpReader :: ~DumpReader( ){
53 < #ifdef IS_MPI
54 <  if (worldRank == 0) {
55 < #endif
56 <  vector<fpos_t*>::iterator i;
57 <
58 <  int error;
59 <  error = fclose( inFile );
60 <  if( error ){
61 <    sprintf( painCave.errMsg,
62 <             "Error closing %s\n", inFileName.c_str());
63 <    simError();
64 <  }
65 <
66 <  for(i = framePos.begin(); i != framePos.end(); ++i)
67 <    delete *i;
68 <  framePos.clear();
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 >
62 > #ifdef IS_MPI
63 >
64 > #include <mpi.h>
65 > #define TAKE_THIS_TAG_CHAR 0
66 > #define TAKE_THIS_TAG_INT 1
67 >
68 > #endif // is_mpi
69 >
70 >
71 > namespace OpenMD {
72 >  
73 >  DumpReader::DumpReader(SimInfo* info, const std::string& filename)
74 >    : info_(info), filename_(filename), isScanned_(false), nframes_(0), needCOMprops_(false) {
75 >    
76 > #ifdef IS_MPI
77 >    
78 >    if (worldRank == 0) {
79 > #endif
80 >      
81 >      inFile_ = new std::ifstream(filename_.c_str());
82 >      
83 >      if (inFile_->fail()) {
84 >        sprintf(painCave.errMsg,
85 >                "DumpReader: Cannot open file: %s\n",
86 >                filename_.c_str());
87 >        painCave.isFatal = 1;
88 >        simError();
89 >      }
90 >      
91 > #ifdef IS_MPI
92 >      
93 >    }
94 >    
95 >    strcpy(checkPointMsg, "Dump file opened for reading successfully.");
96 >    errorCheckPoint();
97 >    
98 > #endif
99 >    
100 >    return;
101 >  }
102    
103 < #ifdef IS_MPI
70 <  }
71 <  strcpy( checkPointMsg, "Dump file closed successfully." );
72 <  MPIcheckPoint();
73 < #endif
74 <
75 <  return;
76 < }
77 <
78 < int DumpReader::getNframes( void ){
79 <
80 <  if( !isScanned )
81 <    scanFile();
82 <  return framePos.size();
83 < }
84 <
85 < void DumpReader::scanFile( void ){
86 <
87 <  int i, j;
88 <  int lineNum = 0;
89 <  char readBuffer[2000];
90 <  fpos_t *currPos;
91 <
92 < #ifdef IS_MPI
93 <  if( worldRank == 0 ){
94 < #endif // is_mpi
103 >  DumpReader::~DumpReader() {
104      
105 <    rewind( inFile );
105 > #ifdef IS_MPI
106      
107 <    currPos = new fpos_t;
108 <    fgetpos( inFile, currPos );
100 <    fgets( readBuffer, sizeof( readBuffer ), inFile );
101 <    lineNum++;
102 <    if( feof( inFile ) ){
103 <      sprintf( painCave.errMsg,
104 <               "File \"%s\" ended unexpectedly at line %d\n",
105 <               inFileName.c_str(),
106 <               lineNum );
107 <      painCave.isFatal = 1;
108 <      simError();
109 <    }
110 <
111 <    while( !feof( inFile ) ){
107 >    if (worldRank == 0) {
108 > #endif
109        
110 <      framePos.push_back(currPos);
114 <
115 <      i = atoi(readBuffer);
110 >      delete inFile_;
111        
112 <      fgets( readBuffer, sizeof( readBuffer ), inFile );
118 <      lineNum++;    
119 <      if( feof( inFile ) ){
120 <        sprintf( painCave.errMsg,
121 <                 "File \"%s\" ended unexpectedly at line %d\n",
122 <                 inFileName.c_str(),
123 <                 lineNum );
124 <        painCave.isFatal = 1;
125 <        simError();
126 <      }
127 <            
128 <      for(j=0; j<i; j++){
129 <        
130 <        fgets( readBuffer, sizeof( readBuffer ), inFile );
131 <        lineNum++;    
132 <        if( feof( inFile ) ){
133 <          sprintf( painCave.errMsg,
134 <                   "File \"%s\" ended unexpectedly at line %d,"
135 <                   " with atom %d\n",
136 <                   inFileName.c_str(),
137 <                   lineNum,
138 <                   j );
139 <          painCave.isFatal = 1;
140 <          simError();
141 <        }
142 <        
143 <      }
112 > #ifdef IS_MPI
113        
114 <      currPos = new fpos_t;
146 <      fgetpos( inFile, currPos );
147 <      fgets( readBuffer, sizeof( readBuffer ), inFile );
148 <      lineNum++;
149 <    }
114 >    }
115      
116 <    delete currPos;
117 <    rewind( inFile );
116 >    strcpy(checkPointMsg, "Dump file closed successfully.");
117 >    errorCheckPoint();
118 >    
119 > #endif
120 >    
121 >    return;
122 >  }
123 >  
124 >  int DumpReader::getNFrames(void) {
125      
126 <    isScanned = true;
127 <
128 < #ifdef IS_MPI
129 <  }
130 <  strcpy( checkPointMsg, "Successfully scanned DumpFile\n" );
131 <  MPIcheckPoint();
132 < #endif // is_mpi
133 < }
134 <
135 < void DumpReader :: readFrame( SimInfo* the_simnfo, int whichFrame){
136 <
137 <  simnfo = the_simnfo;
138 <
139 <  this->readSet( whichFrame );
140 < }
126 >    if (!isScanned_)
127 >      scanFile();
128 >    
129 >    return nframes_;
130 >  }
131 >  
132 >  void DumpReader::scanFile(void) {
133 >    int lineNo = 0;
134 >    std::streampos prevPos;
135 >    std::streampos  currPos;
136 >    
137 > #ifdef IS_MPI
138 >    
139 >    if (worldRank == 0) {
140 > #endif // is_mpi
141 >      
142 >      currPos = inFile_->tellg();
143 >      prevPos = currPos;
144 >      bool foundOpenSnapshotTag = false;
145 >      bool foundClosedSnapshotTag = false;
146 >      while(inFile_->getline(buffer, bufferSize)) {
147 >        ++lineNo;
148 >        
149 >        std::string line = buffer;
150 >        currPos = inFile_->tellg();
151 >        if (line.find("<Snapshot>")!= std::string::npos) {
152 >          if (foundOpenSnapshotTag) {
153 >            sprintf(painCave.errMsg,
154 >                    "DumpReader:<Snapshot> is multiply nested at line %d in %s \n", lineNo,
155 >                    filename_.c_str());
156 >            painCave.isFatal = 1;
157 >            simError();          
158 >          }
159 >          foundOpenSnapshotTag = true;
160 >          foundClosedSnapshotTag = false;
161 >          framePos_.push_back(prevPos);
162 >          
163 >        } else if (line.find("</Snapshot>") != std::string::npos){
164 >          if (!foundOpenSnapshotTag) {
165 >            sprintf(painCave.errMsg,
166 >                    "DumpReader:</Snapshot> appears before <Snapshot> at line %d in %s \n", lineNo,
167 >                    filename_.c_str());
168 >            painCave.isFatal = 1;
169 >            simError();
170 >          }
171 >          
172 >          if (foundClosedSnapshotTag) {
173 >            sprintf(painCave.errMsg,
174 >                    "DumpReader:</Snapshot> appears multiply nested at line %d in %s \n", lineNo,
175 >                    filename_.c_str());
176 >            painCave.isFatal = 1;
177 >            simError();
178 >          }
179 >          foundClosedSnapshotTag = true;
180 >          foundOpenSnapshotTag = false;
181 >        }
182 >        prevPos = currPos;
183 >      }
184 >      
185 >      // only found <Snapshot> for the last frame means the file is corrupted, we should discard
186 >      // it and give a warning message
187 >      if (foundOpenSnapshotTag) {
188 >        sprintf(painCave.errMsg,
189 >                "DumpReader: last frame in %s is invalid\n", filename_.c_str());
190 >        painCave.isFatal = 0;
191 >        simError();      
192 >        framePos_.pop_back();
193 >      }
194 >      
195 >      nframes_ = framePos_.size();
196 >      
197 >      if (nframes_ == 0) {
198 >        sprintf(painCave.errMsg,
199 >                "DumpReader: %s does not contain a valid frame\n", filename_.c_str());
200 >        painCave.isFatal = 1;
201 >        simError();      
202 >      }
203 > #ifdef IS_MPI
204 >    }
205 >    
206 >    MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
207 >    
208 > #endif // is_mpi
209 >    
210 >    isScanned_ = true;
211 >  }
212 >  
213 >  void DumpReader::readFrame(int whichFrame) {
214 >    if (!isScanned_)
215 >      scanFile();
216 >        
217 >    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
218 >    
219 >    if (storageLayout & DataStorage::dslPosition) {
220 >      needPos_ = true;
221 >    } else {
222 >      needPos_ = false;
223 >    }
224 >    
225 >    if (storageLayout & DataStorage::dslVelocity) {
226 >      needVel_ = true;
227 >    } else {
228 >      needVel_ = false;
229 >    }
230 >    
231 >    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
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 +    if (needCOMprops_) {
246 +      Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
247 +      Vector3d com;
248 +      Vector3d comvel;
249 +      Vector3d comw;
250 +      if (needPos_ && needVel_){
251 +        info_->getComAll(com, comvel);
252 +        comw = info_->getAngularMomentum();
253 +      }else{
254 +        com = info_->getCom();
255 +        comvel = 0.0;
256 +        comw   = 0.0;
257 +      }
258 +      s->setCOMprops(com, comvel, comw);      
259 +    }
260  
261 +  }
262 +  
263 +  void DumpReader::readSet(int whichFrame) {    
264 +    std::string line;
265  
266 < void DumpReader :: readSet( int whichFrame ){
266 > #ifndef IS_MPI
267 >    inFile_->clear();  
268 >    inFile_->seekg(framePos_[whichFrame]);
269  
270 <  int i;
175 <  unsigned int j;
270 >    std::istream& inputStream = *inFile_;    
271  
272 < #ifdef IS_MPI
273 <  int done, which_node, which_atom; // loop counter
274 < #endif //is_mpi
272 > #else
273 >    int masterNode = 0;
274 >    std::stringstream sstream;
275 >    if (worldRank == masterNode) {
276 >      std::string sendBuffer;
277  
278 <  const int BUFFERSIZE = 2000; // size of the read buffer
279 <  int nTotObjs; // the number of atoms
280 <  char read_buffer[BUFFERSIZE]; //the line buffer for reading
278 >      inFile_->clear();  
279 >      inFile_->seekg(framePos_[whichFrame]);
280 >      
281 >      while (inFile_->getline(buffer, bufferSize)) {
282  
283 <  char *eof_test; // ptr to see when we reach the end of the file
284 <  char *parseErr;
283 >        line = buffer;
284 >        sendBuffer += line;
285 >        sendBuffer += '\n';
286 >        if (line.find("</Snapshot>") != std::string::npos) {
287 >          break;
288 >        }        
289 >      }
290  
291 <  vector<StuntDouble*> integrableObjects;
291 >      int sendBufferSize = sendBuffer.size();
292 >      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);    
293 >      MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);    
294 >      
295 >      sstream.str(sendBuffer);
296 >    } else {
297 >      int sendBufferSize;
298 >      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);    
299 >      char * recvBuffer = new char[sendBufferSize+1];
300 >      assert(recvBuffer);
301 >      recvBuffer[sendBufferSize] = '\0';
302 >      MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);    
303 >      sstream.str(recvBuffer);
304 >      delete [] recvBuffer;
305 >    }      
306  
307 +    std::istream& inputStream = sstream;  
308 + #endif
309  
310 < #ifndef IS_MPI
310 >    inputStream.getline(buffer, bufferSize);
311  
312 <  fsetpos(inFile, framePos[whichFrame]);
313 <  eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
314 <  if( eof_test == NULL ){
315 <    sprintf( painCave.errMsg,
316 <             "DumpReader error: error reading 1st line of \"%s\"\n",
317 <             inFileName.c_str() );
318 <    painCave.isFatal = 1;
319 <    simError();
320 <  }
312 >    line = buffer;
313 >    if (line.find("<Snapshot>") == std::string::npos) {
314 >      sprintf(painCave.errMsg,
315 >              "DumpReader Error: can not find <Snapshot>\n");
316 >      painCave.isFatal = 1;
317 >      simError();
318 >    }
319 >    
320 >    //read frameData
321 >    readFrameProperties(inputStream);
322  
323 <  nTotObjs = atoi( read_buffer );
323 >    //read StuntDoubles
324 >    readStuntDoubles(inputStream);    
325  
326 <  if( nTotObjs != simnfo->getTotIntegrableObjects() ){
327 <    sprintf( painCave.errMsg,
328 <             "DumpReader error. %s nIntegrable, %d, "
329 <             "does not match the meta-data file's nIntegrable, %d.\n",
330 <             inFileName.c_str(), nTotObjs, simnfo->getTotIntegrableObjects());
331 <    painCave.isFatal = 1;
332 <    simError();
333 <  }
326 >    inputStream.getline(buffer, bufferSize);
327 >    line = buffer;
328 >    if (line.find("</Snapshot>") == std::string::npos) {
329 >      sprintf(painCave.errMsg,
330 >              "DumpReader Error: can not find </Snapshot>\n");
331 >      painCave.isFatal = 1;
332 >      simError();
333 >    }        
334 >  
335 >  }
336 >  
337 >  void DumpReader::parseDumpLine(const std::string& line) {
338  
339 <  //read the box mat from the comment line
339 >      
340 >    StringTokenizer tokenizer(line);
341 >    int nTokens;
342 >    
343 >    nTokens = tokenizer.countTokens();
344 >    
345 >    if (nTokens < 2) {  
346 >      sprintf(painCave.errMsg,
347 >              "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
348 >      painCave.isFatal = 1;
349 >      simError();
350 >    }
351  
352 <  eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
353 <  if(eof_test == NULL){
354 <    sprintf( painCave.errMsg,
219 <             "error in reading commment in %s\n", inFileName.c_str());
220 <    painCave.isFatal = 1;
221 <    simError();
222 <  }
352 >    int index = tokenizer.nextTokenAsInt();
353 >
354 >    StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
355  
356 <  parseErr = parseCommentLine( read_buffer, simnfo);
357 <  if( parseErr != NULL ){
358 <    strcpy( painCave.errMsg, parseErr );
359 <    painCave.isFatal = 1;
360 <    simError();
229 <  }
356 >    if (integrableObject == NULL) {
357 >      return;
358 >    }
359 >    std::string type = tokenizer.nextToken();
360 >    int size = type.size();
361  
362 <  //parse dump lines
362 >    size_t found;
363 >    
364 >    if (needPos_) {
365 >      found = type.find("p");      
366 >      if (found == std::string::npos) {
367 >        sprintf(painCave.errMsg,
368 >                "DumpReader Error: StuntDouble %d has no Position\n"
369 >                "\tField (\"p\") specified.\n%s\n", index,
370 >                line.c_str());  
371 >        painCave.isFatal = 1;
372 >        simError();
373 >      }
374 >    }
375 >    
376 >    if (integrableObject->isDirectional()) {
377 >      if (needQuaternion_) {
378 >        found = type.find("q");      
379 >        if (found == std::string::npos) {
380 >          sprintf(painCave.errMsg,
381 >                  "DumpReader Error: Directional StuntDouble %d has no\n"
382 >                  "\tQuaternion Field (\"q\") specified.\n%s\n", index,
383 >                  line.c_str());  
384 >          painCave.isFatal = 1;
385 >          simError();
386 >        }
387 >      }      
388 >    }
389  
390 <  for( i=0; i < simnfo->n_mol; i++){
390 >    for(int i = 0; i < size; ++i) {
391 >      switch(type[i]) {
392 >        
393 >        case 'p': {
394 >            Vector3d pos;
395 >            pos[0] = tokenizer.nextTokenAsDouble();
396 >            pos[1] = tokenizer.nextTokenAsDouble();
397 >            pos[2] = tokenizer.nextTokenAsDouble();
398 >            if (needPos_) {
399 >              integrableObject->setPos(pos);
400 >            }            
401 >            break;
402 >        }
403 >        case 'v' : {
404 >            Vector3d vel;
405 >            vel[0] = tokenizer.nextTokenAsDouble();
406 >            vel[1] = tokenizer.nextTokenAsDouble();
407 >            vel[2] = tokenizer.nextTokenAsDouble();
408 >            if (needVel_) {
409 >              integrableObject->setVel(vel);
410 >            }
411 >            break;
412 >        }
413  
414 <    integrableObjects = (simnfo->molecules[i]).getIntegrableObjects();
414 >        case 'q' : {
415 >           Quat4d q;
416 >           if (integrableObject->isDirectional()) {
417 >              
418 >             q[0] = tokenizer.nextTokenAsDouble();
419 >             q[1] = tokenizer.nextTokenAsDouble();
420 >             q[2] = tokenizer.nextTokenAsDouble();
421 >             q[3] = tokenizer.nextTokenAsDouble();
422 >              
423 >             RealType qlen = q.length();
424 >             if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
425 >                
426 >               sprintf(painCave.errMsg,
427 >                       "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
428 >               painCave.isFatal = 1;
429 >               simError();
430 >                
431 >             }  
432 >              
433 >             q.normalize();
434 >             if (needQuaternion_) {            
435 >               integrableObject->setQ(q);
436 >             }              
437 >           }            
438 >           break;
439 >        }  
440 >        case 'j' : {
441 >          Vector3d ji;
442 >          if (integrableObject->isDirectional()) {
443 >             ji[0] = tokenizer.nextTokenAsDouble();
444 >             ji[1] = tokenizer.nextTokenAsDouble();
445 >             ji[2] = tokenizer.nextTokenAsDouble();
446 >             if (needAngMom_) {
447 >               integrableObject->setJ(ji);
448 >             }
449 >          }
450 >          break;
451 >        }  
452 >        case 'f': {
453 >
454 >          Vector3d force;
455 >          force[0] = tokenizer.nextTokenAsDouble();
456 >          force[1] = tokenizer.nextTokenAsDouble();
457 >          force[2] = tokenizer.nextTokenAsDouble();          
458 >          integrableObject->setFrc(force);
459 >          break;
460 >        }
461 >        case 't' : {
462 >
463 >           Vector3d torque;
464 >           torque[0] = tokenizer.nextTokenAsDouble();
465 >           torque[1] = tokenizer.nextTokenAsDouble();
466 >           torque[2] = tokenizer.nextTokenAsDouble();          
467 >           integrableObject->setTrq(torque);          
468 >           break;
469 >        }
470 >        case 'u' : {
471 >
472 >           RealType particlePot;
473 >           particlePot = tokenizer.nextTokenAsDouble();
474 >           integrableObject->setParticlePot(particlePot);          
475 >           break;
476 >        }
477 >        case 'c' : {
478  
479 <    for(j = 0; j < integrableObjects.size(); j++){
479 >           RealType flucQPos;
480 >           flucQPos = tokenizer.nextTokenAsDouble();
481 >           integrableObject->setFlucQPos(flucQPos);          
482 >           break;
483 >        }
484 >        case 'w' : {
485  
486 <      eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
487 <      if(eof_test == NULL){
488 <        sprintf(painCave.errMsg,
489 <              "error in reading file %s\n"
490 <              "natoms  = %d; index = %d\n"
491 <              "error reading the line from the file.\n",
245 <              inFileName.c_str(), nTotObjs, i );
246 <        painCave.isFatal = 1;
247 <        simError();
248 <      }
249 <      
250 <      parseErr = parseDumpLine( read_buffer, integrableObjects[j]);
251 <      if( parseErr != NULL ){
252 <        strcpy( painCave.errMsg, parseErr );
253 <        painCave.isFatal = 1;
254 <        simError();
255 <      }
256 <    }
257 <  }
486 >           RealType flucQVel;
487 >           flucQVel = tokenizer.nextTokenAsDouble();
488 >           integrableObject->setFlucQVel(flucQVel);          
489 >           break;
490 >        }
491 >        case 'g' : {
492  
493 <  // MPI Section of code..........
494 < #else //IS_MPI
493 >           RealType flucQFrc;
494 >           flucQFrc = tokenizer.nextTokenAsDouble();
495 >           integrableObject->setFlucQFrc(flucQFrc);          
496 >           break;
497 >        }
498 >        case 'e' : {
499  
500 <  // first thing first, suspend fatalities.
501 <  painCave.isEventLoop = 1;
500 >           Vector3d eField;
501 >           eField[0] = tokenizer.nextTokenAsDouble();
502 >           eField[1] = tokenizer.nextTokenAsDouble();
503 >           eField[2] = tokenizer.nextTokenAsDouble();          
504 >           integrableObject->setElectricField(eField);          
505 >           break;
506 >        }
507 >        default: {
508 >               sprintf(painCave.errMsg,
509 >                       "DumpReader Error: %s is an unrecognized type\n", type.c_str());
510 >               painCave.isFatal = 1;
511 >               simError();
512 >          break;  
513 >        }
514  
515 <  int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
266 <  int haveError;
267 <
268 <  MPI_Status istatus;
269 <  int *MolToProcMap = mpiSim->getMolToProcMap();
270 <  int localIndex;
271 <  int nCurObj;
272 <  int nitems;
273 <
274 <  nTotObjs = simnfo->getTotIntegrableObjects();
275 <  haveError = 0;
276 <  if (worldRank == 0) {
277 <     fsetpos(inFile,  framePos[whichFrame]);
278 <
279 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
280 <    if( eof_test == NULL ){
281 <      sprintf( painCave.errMsg,
282 <               "Error reading 1st line of %s \n ",inFileName.c_str());
283 <      haveError = 1;
284 <      simError();
515 >      }
516      }
517 +    
518 +  }
519 +  
520  
521 <    nitems = atoi( read_buffer );
521 >  void  DumpReader::readStuntDoubles(std::istream& inputStream) {
522  
523 <    // Check to see that the number of integrable objects in the
524 <    // intial configuration file is the same as derived from the
525 <    // meta-data file.
526 <
527 <    if( nTotObjs != nitems){
528 <      sprintf( painCave.errMsg,
529 <               "DumpReader Error. %s nIntegrable, %d, "
530 <               "does not match the meta-data file's nIntegrable, %d.\n",
297 <               inFileName.c_str(), nTotObjs, simnfo->getTotIntegrableObjects());
298 <      haveError= 1;
299 <      simError();
523 >    inputStream.getline(buffer, bufferSize);
524 >    std::string line(buffer);
525 >    
526 >    if (line.find("<StuntDoubles>") == std::string::npos) {
527 >      sprintf(painCave.errMsg,
528 >              "DumpReader Error: Missing <StuntDoubles>\n");
529 >      painCave.isFatal = 1;
530 >      simError();
531      }
532  
533 <    //read the boxMat from the comment line
534 <
304 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
305 <    if(eof_test == NULL){
306 <      sprintf( painCave.errMsg,
307 <               "error in reading commment in %s\n", inFileName.c_str());
308 <      haveError = 1;
309 <      simError();
310 <    }
311 <
312 <    //Every single processor will parse the comment line by itself
313 <    //By using this way, we might lose some efficiency, but if we want to add
314 <    //more parameters into comment line, we only need to modify function
315 <    //parseCommentLine
316 <
317 <    MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, 0, MPI_COMM_WORLD);
318 <
319 <    parseErr = parseCommentLine( read_buffer, simnfo);
320 <
321 <    if( parseErr != NULL ){
322 <      strcpy( painCave.errMsg, parseErr );
323 <      haveError = 1;
324 <      simError();
325 <    }
326 <
327 <    for (i=0 ; i < mpiSim->getNMolGlobal(); i++) {
328 <      which_node = MolToProcMap[i];
329 <      if(which_node == 0){
330 <       //molecules belong to master node
331 <
332 <      localIndex = mpiSim->getGlobalToLocalMol(i);
333 <
334 <      if(localIndex == -1) {
335 <        strcpy(painCave.errMsg, "Molecule not found on node 0!");
336 <        haveError = 1;
337 <        simError();
338 <      }
339 <
340 <       integrableObjects = (simnfo->molecules[localIndex]).getIntegrableObjects();
341 <       for(j=0; j < integrableObjects.size(); j++){
342 <        
343 <          eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
344 <          if(eof_test == NULL){
345 <                sprintf(painCave.errMsg,
346 <                    "error in reading file %s\n"
347 <                    "natoms  = %d; index = %d\n"
348 <                    "error reading the line from the file.\n",
349 <                    inFileName.c_str(), nTotObjs, i );
350 <                haveError= 1;
351 <                simError();
352 <          }
353 <          
354 <          if(haveError) nodeZeroError();
355 <
356 <          parseDumpLine(read_buffer, integrableObjects[j]);
357 <          
358 <       }
359 <
360 <
361 <      }
362 <      else{
363 <      //molecule belongs to slave nodes
364 <
365 <        MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
366 <               TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
533 >    while(inputStream.getline(buffer, bufferSize)) {
534 >      line = buffer;
535        
536 <       for(j=0; j < nCurObj; j++){
537 <        
370 <          eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
371 <          if(eof_test == NULL){
372 <                sprintf(painCave.errMsg,
373 <                    "error in reading file %s\n"
374 <                    "natoms  = %d; index = %d\n"
375 <                    "error reading the line from the file.\n",
376 <                    inFileName.c_str(), nTotObjs, i );
377 <                haveError= 1;
378 <                simError();
379 <          }
380 <          
381 <          if(haveError) nodeZeroError();
382 <
383 <            MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
384 <                      TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
385 <          
386 <       }
387 <
536 >      if(line.find("</StuntDoubles>") != std::string::npos) {
537 >        break;
538        }
389      
390    }
391    
392  }
393  else{
394  //actions taken at slave nodes
395    MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, 0, MPI_COMM_WORLD);
539  
540 <    parseErr = parseCommentLine( read_buffer, simnfo);
398 <
399 <    if( parseErr != NULL ){
400 <      strcpy( painCave.errMsg, parseErr );
401 <      haveError = 1;
402 <      simError();
540 >      parseDumpLine(line);
541      }
542    
543 <    for (i=0 ; i < mpiSim->getNMolGlobal(); i++) {
406 <      which_node = MolToProcMap[i];
407 <      
408 <      if(which_node == worldRank){
409 <      //molecule with global index i belongs to this processor
410 <      
411 <        localIndex = mpiSim->getGlobalToLocalMol(i);
543 >  }
544  
545 <        if(localIndex == -1) {
414 <          sprintf(painCave.errMsg, "Molecule not found on node %d\n", worldRank);
415 <          haveError = 1;
416 <          simError();
417 <        }
545 >  void DumpReader::readFrameProperties(std::istream& inputStream) {
546  
547 <        integrableObjects = (simnfo->molecules[localIndex]).getIntegrableObjects();        
547 >    Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
548 >    inputStream.getline(buffer, bufferSize);
549 >    std::string line(buffer);
550  
551 <        nCurObj = integrableObjects.size();
552 <        
553 <        MPI_Send(&nCurObj, 1, MPI_INT, 0,
554 <                        TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
551 >    if (line.find("<FrameData>") == std::string::npos) {
552 >      sprintf(painCave.errMsg,
553 >              "DumpReader Error: Missing <FrameData>\n");
554 >      painCave.isFatal = 1;
555 >      simError();
556 >    }
557  
558 <        for(j = 0; j < integrableObjects.size(); j++){
559 <
560 <          MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, 0,
561 <                              TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
562 <
431 <          parseErr = parseDumpLine(read_buffer, integrableObjects[j]);
432 <
433 <          if( parseErr != NULL ){
434 <                strcpy( painCave.errMsg, parseErr );
435 <                simError();
436 <          }
437 <
438 <        }
439 <          
558 >    while(inputStream.getline(buffer, bufferSize)) {
559 >      line = buffer;
560 >      
561 >      if(line.find("</FrameData>") != std::string::npos) {
562 >        break;
563        }
564        
565 <    }
566 <
567 <  }
568 <
569 < #endif
570 < }
571 <
449 < char* DumpReader::parseDumpLine(char* readLine, StuntDouble* sd){
450 <
451 <  char *foo; // the pointer to the current string token
452 <
453 <  double pos[3]; // position place holders
454 <  double vel[3]; // velocity placeholders
455 <  double q[4]; // the quaternions
456 <  double ji[3]; // angular velocity placeholders;
457 <  double qSqr, qLength; // needed to normalize the quaternion vector.
458 <
459 <
460 <  // set the string tokenizer
461 <
462 <  foo = strtok(readLine, " ,;\t");
463 <
464 <  // check the atom name to the current atom
465 <
466 <  if( strcmp( foo, sd->getType() ) ){
467 <    sprintf( painCave.errMsg,
468 <             "DumpReader error.  Does not"
469 <             " match the meta-data atom %s.\n",
470 <             sd->getType() );
471 <    return strdup( painCave.errMsg );
472 <  }
473 <
474 <  // get the positions
475 <
476 <  foo = strtok(NULL, " ,;\t");
477 <  if(foo == NULL){
478 <    sprintf( painCave.errMsg,
479 <             "error in reading postition x from %s\n",
480 <             inFileName.c_str());
481 <    return strdup( painCave.errMsg );
482 <  }
483 <  pos[0] = atof( foo );
484 <
485 <  foo = strtok(NULL, " ,;\t");
486 <  if(foo == NULL){
487 <    sprintf( painCave.errMsg,
488 <             "error in reading postition y from %s\n",
489 <             inFileName.c_str());
490 <    return strdup( painCave.errMsg );
491 <  }
492 <  pos[1] = atof( foo );
565 >      StringTokenizer tokenizer(line, " ;\t\n\r{}:,");
566 >      if (!tokenizer.hasMoreTokens()) {
567 >        sprintf(painCave.errMsg,
568 >                "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
569 >        painCave.isFatal = 1;
570 >        simError();      
571 >      }
572  
573 <  foo = strtok(NULL, " ,;\t");
574 <  if(foo == NULL){
575 <    sprintf( painCave.errMsg,
576 <             "error in reading postition z from %s\n",
577 <             inFileName.c_str());
578 <    return strdup( painCave.errMsg );
579 <  }
580 <  pos[2] = atof( foo );
581 <
582 <
583 <  // get the velocities
584 <
585 <  foo = strtok(NULL, " ,;\t");
586 <  if(foo == NULL){
587 <    sprintf( painCave.errMsg,
588 <             "error in reading velocity x from %s\n",
589 <             inFileName.c_str() );
590 <    return strdup( painCave.errMsg );
591 <  }
592 <  vel[0] = atof( foo );
593 <
594 <  foo = strtok(NULL, " ,;\t");
595 <  if(foo == NULL){
596 <    sprintf( painCave.errMsg,
597 <             "error in reading velocity x from %s\n",
598 <             inFileName.c_str() );
599 <    return strdup( painCave.errMsg );
600 <  }
601 <  vel[1] = atof( foo );
602 <
603 <  foo = strtok(NULL, " ,;\t");
604 <  if(foo == NULL){
605 <    sprintf( painCave.errMsg,
606 <             "error in reading velocity x from %s\n",
607 <             inFileName.c_str() );
608 <    return strdup( painCave.errMsg );
609 <  }
610 <  vel[2] = atof( foo );
532 <
533 <
534 <  // add the positions and velocities to the atom
535 <
536 <  sd->setPos( pos );
537 <  sd->setVel( vel );
538 <
539 <  if (!sd->isDirectional())
540 <    return NULL;
541 <
542 <  // get the quaternions
543 <
544 <  if( sd->isDirectional() ){
545 <
546 <    foo = strtok(NULL, " ,;\t");
547 <    if(foo == NULL){
548 <      sprintf( painCave.errMsg,
549 <                     "error in reading velocity x from %s\n",
550 <                      inFileName.c_str() );
551 <      return strdup( painCave.errMsg );
552 <    }
553 <    q[0] = atof( foo );
554 <
555 <    foo = strtok(NULL, " ,;\t");
556 <    if(foo == NULL){
557 <      sprintf( painCave.errMsg,
558 <                     "error in reading velocity x from %s\n",
559 <                      inFileName.c_str() );
560 <      return strdup( painCave.errMsg );
561 <    }
562 <    q[1] = atof( foo );
563 <
564 <    foo = strtok(NULL, " ,;\t");
565 <    if(foo == NULL){
566 <      sprintf( painCave.errMsg,
567 <                     "error in reading velocity x from %s\n",
568 <                      inFileName.c_str() );
569 <      return strdup( painCave.errMsg );
570 <    }
571 <    q[2] = atof( foo );
572 <
573 <    foo = strtok(NULL, " ,;\t");
574 <    if(foo == NULL){
575 <      sprintf( painCave.errMsg,
576 <                     "error in reading velocity x from %s\n",
577 <                      inFileName.c_str() );
578 <      return strdup( painCave.errMsg );
579 <    }
580 <    q[3] = atof( foo );
581 <
582 <    // get the angular velocities
583 <
584 <    foo = strtok(NULL, " ,;\t");
585 <    if(foo == NULL){
586 <      sprintf( painCave.errMsg,
587 <                     "error in reading velocity x from %s\n",
588 <                      inFileName.c_str() );
589 <      return strdup( painCave.errMsg );
590 <    }
591 <    ji[0] = atof( foo );
592 <
593 <    foo = strtok(NULL, " ,;\t");
594 <    if(foo == NULL){
595 <      sprintf( painCave.errMsg,
596 <                     "error in reading velocity x from %s\n",
597 <                      inFileName.c_str() );
598 <      return strdup( painCave.errMsg );
599 <    }
600 <    ji[1] = atof(foo );
601 <
602 <    foo = strtok(NULL, " ,;\t");
603 <    if(foo == NULL){
604 <      sprintf( painCave.errMsg,
605 <                     "error in reading velocity x from %s\n",
606 <                      inFileName.c_str() );
607 <      return strdup( painCave.errMsg );
608 <    }
609 <    ji[2] = atof( foo );
610 <
611 <
612 <    // check that the quaternion vector is normalized
613 <
614 <    qSqr = (q[0] * q[0]) + (q[1] * q[1]) + (q[2] * q[2]) + (q[3] * q[3]);
615 <
616 <    if (fabs(qSqr) < 1e-6) {
617 <      sprintf(painCave.errMsg,
618 <          "initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0).\n");
619 <       return strdup(painCave.errMsg);
620 <    }
621 <
622 <    qLength = sqrt( qSqr );
623 <    q[0] = q[0] / qLength;
624 <    q[1] = q[1] / qLength;
625 <    q[2] = q[2] / qLength;
626 <    q[3] = q[3] / qLength;
627 <
628 <    // add quaternion and angular velocities
629 <
630 <    sd->setQ( q );
631 <    sd->setJ( ji );
632 <  }
633 <
634 <
635 <
636 <  return NULL;
637 < }
638 <
639 <
640 < char* DumpReader::parseCommentLine(char* readLine, SimInfo* entry_plug){
641 <
642 <  double currTime;
643 <  double boxMat[9];
644 <  double theBoxMat3[3][3];
645 <  double chi;
646 <  double integralOfChidt;
647 <  double eta[9];
648 <
649 <  char *foo; // the pointer to the current string token
650 <
651 <  // set the string tokenizer
652 <
653 <  foo = strtok(readLine, " ,;\t");
654 <  // set the timeToken.
655 <
656 <  if(foo == NULL){
657 <    sprintf( painCave.errMsg,
658 <             "error in reading Time from %s\n",
659 <             inFileName.c_str() );
660 <    return strdup( painCave.errMsg );
661 <  }
662 <
663 <  currTime = atof( foo );
664 <  entry_plug->setTime( currTime );
665 <
666 <  //get H-Matrix
667 <
668 <  for(int i = 0 ; i < 9; i++){
669 <    foo = strtok(NULL, " ,;\t");
670 <    if(foo == NULL){
671 <      sprintf( painCave.errMsg,
672 <               "error in reading H[%d] from %s\n", i, inFileName.c_str() );
673 <      return strdup( painCave.errMsg );
674 <    }
675 <    boxMat[i] = atof( foo );
676 <  }
677 <
678 <  for(int i=0;i<3;i++)
679 <    for(int j=0;j<3;j++) theBoxMat3[i][j] = boxMat[3*j+i];
680 <
681 <  //set H-Matrix
682 <  entry_plug->setBoxM( theBoxMat3 );
683 <
684 <  //get chi and integralOfChidt, they should appear by pair
685 <
686 <  if( entry_plug->useInitXSstate ){
687 <    foo = strtok(NULL, " ,;\t\n");
688 <    if(foo != NULL){
689 <      chi = atof(foo);
690 <      
691 <      foo = strtok(NULL, " ,;\t\n");
692 <      if(foo == NULL){
693 <        sprintf( painCave.errMsg,
694 <                 "chi and integralOfChidt should appear by pair in %s\n", inFileName.c_str() );
695 <        return strdup( painCave.errMsg );
573 >      std::string propertyName = tokenizer.nextToken();
574 >      if (propertyName == "Time") {
575 >        RealType currTime = tokenizer.nextTokenAsDouble();
576 >        s->setTime(currTime);
577 >      } else if (propertyName == "Hmat"){
578 >        Mat3x3d hmat;
579 >        hmat(0, 0) = tokenizer.nextTokenAsDouble();
580 >        hmat(0, 1) = tokenizer.nextTokenAsDouble();
581 >        hmat(0, 2) = tokenizer.nextTokenAsDouble();
582 >        hmat(1, 0) = tokenizer.nextTokenAsDouble();
583 >        hmat(1, 1) = tokenizer.nextTokenAsDouble();
584 >        hmat(1, 2) = tokenizer.nextTokenAsDouble();
585 >        hmat(2, 0) = tokenizer.nextTokenAsDouble();
586 >        hmat(2, 1) = tokenizer.nextTokenAsDouble();
587 >        hmat(2, 2) = tokenizer.nextTokenAsDouble();
588 >        s->setHmat(hmat);      
589 >      } else if (propertyName == "Thermostat") {
590 >        RealType chi = tokenizer.nextTokenAsDouble();
591 >        RealType integralOfChiDt = tokenizer.nextTokenAsDouble();
592 >        s->setChi(chi);
593 >        s->setIntegralOfChiDt(integralOfChiDt);        
594 >     } else if (propertyName == "Barostat") {
595 >        Mat3x3d eta;
596 >        eta(0, 0) = tokenizer.nextTokenAsDouble();
597 >        eta(0, 1) = tokenizer.nextTokenAsDouble();
598 >        eta(0, 2) = tokenizer.nextTokenAsDouble();
599 >        eta(1, 0) = tokenizer.nextTokenAsDouble();
600 >        eta(1, 1) = tokenizer.nextTokenAsDouble();
601 >        eta(1, 2) = tokenizer.nextTokenAsDouble();
602 >        eta(2, 0) = tokenizer.nextTokenAsDouble();
603 >        eta(2, 1) = tokenizer.nextTokenAsDouble();
604 >        eta(2, 2) = tokenizer.nextTokenAsDouble();
605 >        s->setEta(eta);
606 >      } else {
607 >        sprintf(painCave.errMsg,
608 >                "DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str());
609 >        painCave.isFatal = 0;
610 >        simError();        
611        }
697      integralOfChidt = atof( foo );
612        
699      //push chi and integralOfChidt into SimInfo::properties which can be
700      //retrieved by integrator later
701      DoubleGenericData* chiValue = new DoubleGenericData();
702      chiValue->setID(CHIVALUE_ID);
703      chiValue->setData(chi);
704      entry_plug->addProperty(chiValue);
705      
706      DoubleGenericData* integralOfChidtValue = new DoubleGenericData();
707      integralOfChidtValue->setID(INTEGRALOFCHIDT_ID);
708      integralOfChidtValue->setData(integralOfChidt);
709      entry_plug->addProperty(integralOfChidtValue);
710      
613      }
712    else
713      return NULL;
714    
715    //get eta
716    foo = strtok(NULL, " ,;\t\n");
717    if(foo != NULL ){
718  
719      for(int i = 0 ; i < 9; i++){
720        
721        if(foo == NULL){
722          sprintf( painCave.errMsg,
723                   "error in reading eta[%d] from %s\n", i, inFileName.c_str() );
724          return strdup( painCave.errMsg );
725        }
726        eta[i] = atof( foo );
727        foo = strtok(NULL, " ,;\t\n");
728      }
729    }
730    else
731      return NULL;
732    
733    //push eta into SimInfo::properties which can be
734    //retrieved by integrator later
735    //entry_plug->setBoxM( theBoxMat3 );
736    DoubleVectorGenericData* etaValue = new DoubleVectorGenericData(ETAVALUE_ID);
737    etaValue->insert(etaValue->end(), eta, eta+9);
738    entry_plug->addProperty(etaValue);
739  }
614  
741  return NULL;
742 }
743
744 #ifdef IS_MPI
745 void DumpReader::nodeZeroError( void ){
746  int j, myStatus;
747
748  myStatus = 0;
749  for (j = 0; j < mpiSim->getNProcessors(); j++) {
750    MPI_Send( &myStatus, 1, MPI_INT, j,
751              TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
615    }
616  
617 <
618 <  MPI_Finalize();
756 <  exit (0);
757 <
758 < }
759 <
760 < void DumpReader::anonymousNodeDie( void ){
761 <
762 <  MPI_Finalize();
763 <  exit (0);
764 < }
765 < #endif
617 >  
618 > }//end namespace OpenMD

Comparing:
trunk/src/io/DumpReader.cpp (property svn:keywords), Revision 132 by tim, Thu Oct 21 16:22:01 2004 UTC vs.
branches/development/src/io/DumpReader.cpp (property svn:keywords), Revision 1714 by gezelter, Sat May 19 18:12:46 2012 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines