src/io/DumpReader.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
#define _LARGEFILE_SOURCE64
#define _FILE_OFFSET_BITS 64

#include <sys/types.h>
#include <sys/stat.h>

#include <iostream>
#include <math.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "io/DumpReader.hpp"
#include "primitives/Molecule.hpp"
#include "utils/simError.h"
#include "utils/MemoryUtils.hpp"
#include "utils/StringTokenizer.hpp"

#ifdef IS_MPI

#include <mpi.h>
#define TAKE_THIS_TAG_CHAR 0
#define TAKE_THIS_TAG_INT 1

#endif // is_mpi


namespace oopse {
  
  DumpReader::DumpReader(SimInfo* info, const std::string& filename)
    : info_(info), filename_(filename), isScanned_(false), nframes_(0) {
    
#ifdef IS_MPI
    
      if (worldRank == 0) {
#endif
      
      inFile_ = new std::ifstream(filename_.c_str());
      
        if (inFile_->fail()) {
          sprintf(painCave.errMsg, "DumpReader: Cannot open file: %s\n", filename_.c_str());
          painCave.isFatal = 1;
          simError();
        }
      
#ifdef IS_MPI
      
      }
    
      strcpy(checkPointMsg, "Dump file opened for reading successfully.");
      MPIcheckPoint();
    
#endif
    
      return;
    }
  
  DumpReader::~DumpReader() {
    
#ifdef IS_MPI
    
    if (worldRank == 0) {
#endif
      
      delete inFile_;
      
#ifdef IS_MPI
      
    }
    
    strcpy(checkPointMsg, "Dump file closed successfully.");
    MPIcheckPoint();
    
#endif
    
    return;
  }
  
  int DumpReader::getNFrames(void) {
    
    if (!isScanned_)
      scanFile();
    
    return nframes_;
  }
  
  void DumpReader::scanFile(void) {
    int i, j;
    int lineNum = 0;
    char readBuffer[maxBufferSize];
    std::streampos  currPos;
    
#ifdef IS_MPI
    
    if (worldRank == 0) {
#endif // is_mpi
      
      inFile_->seekg (0, std::ios::beg);
      

      currPos = inFile_->tellg();
      inFile_->getline(readBuffer, sizeof(readBuffer));
      lineNum++;
      
      if (inFile_->eof()) {
        sprintf(painCave.errMsg,
                "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
                filename_.c_str(),
                lineNum);
        painCave.isFatal = 1;
        simError();
      }
      
      while (!inFile_->eof()) {
        framePos_.push_back(currPos);
        
        i = atoi(readBuffer);
        
        inFile_->getline(readBuffer, sizeof(readBuffer));
        lineNum++;
        
        if (inFile_->eof()) {
          sprintf(painCave.errMsg,
                  "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
                  filename_.c_str(),
                  lineNum);
          painCave.isFatal = 1;
          simError();
        }
        
        for(j = 0; j < i; j++) {
          inFile_->getline(readBuffer, sizeof(readBuffer));
          lineNum++;
          
          if (inFile_->eof()) {
            sprintf(painCave.errMsg,
                    "DumpReader Error: File \"%s\" ended unexpectedly at line %d,"
                    " with atom %d\n", filename_.c_str(),
                    lineNum,
                    j);
            
            painCave.isFatal = 1;
            simError();
          }
        }
        
        currPos = inFile_->tellg();
        inFile_->getline(readBuffer, sizeof(readBuffer));
        lineNum++;
      }

      inFile_->seekg (0, std::ios::beg);
      
      nframes_ = framePos_.size();
#ifdef IS_MPI
    }
    
    MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
    
    strcpy(checkPointMsg, "Successfully scanned DumpFile\n");
    MPIcheckPoint();
    
#endif // is_mpi
    
    isScanned_ = true;
  }
  
  void DumpReader::readFrame(int whichFrame) {
    if (!isScanned_)
      scanFile();
        
    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
    
    if (storageLayout & DataStorage::dslPosition) {
      needPos_ = true;
    } else {
      needPos_ = false;
    }
    
    if (storageLayout & DataStorage::dslVelocity) {
      needVel_ = true;
    } else {
      needVel_ = false;
    }
    
    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
      needQuaternion_ = true;
    } else {
      needQuaternion_ = false;
    }
    
    if (storageLayout & DataStorage::dslAngularMomentum) {
      needAngMom_ = true;
    } else {
      needAngMom_ = false;    
    }
    
    readSet(whichFrame);
  }
  
  void DumpReader::readSet(int whichFrame) {
    int i;
    int nTotObjs;                  // the number of atoms
    char read_buffer[maxBufferSize];  //the line buffer for reading
    char * eof_test;               // ptr to see when we reach the end of the file
    
    Molecule* mol;
    StuntDouble* integrableObject;
    SimInfo::MoleculeIterator mi;
    Molecule::IntegrableObjectIterator ii;
    
#ifndef IS_MPI
    inFile_->clear(); 
    inFile_->seekg(framePos_[whichFrame]);
        
    if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
      sprintf(painCave.errMsg,
              "DumpReader error: error reading 1st line of \"%s\"\n",
              filename_.c_str());
      painCave.isFatal = 1;
      simError();
    }
    
    nTotObjs = atoi(read_buffer);
    
    if (nTotObjs != info_->getNGlobalIntegrableObjects()) {
      sprintf(painCave.errMsg,
              "DumpReader error. %s nIntegrable, %d, "
              "does not match the meta-data file's nIntegrable, %d.\n",
              filename_.c_str(),
              nTotObjs,
              info_->getNGlobalIntegrableObjects());
      
      painCave.isFatal = 1;
      simError();
    }
    
    //read the box mat from the comment line
    
    
    if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
      sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
              filename_.c_str());
      painCave.isFatal = 1;
      simError();
    }
    
    parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
    
    //parse dump lines
    
    i = 0;
    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
      
      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
           integrableObject = mol->nextIntegrableObject(ii)) {           
        
        
        if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
          sprintf(painCave.errMsg,
                  "DumpReader Error: error in reading file %s\n"
                  "natoms  = %d; index = %d\n"
                  "error reading the line from the file.\n",
                  filename_.c_str(),
                  nTotObjs,
                  i);
          
          painCave.isFatal = 1;
          simError();
        }
        
        parseDumpLine(read_buffer, integrableObject);
        i++;
      }
    }
    
    // MPI Section of code..........
    
#else //IS_MPI
    
    // first thing first, suspend fatalities.
    int masterNode = 0;
    int nCurObj;
    painCave.isEventLoop = 1;
    
    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
    int haveError;
    
    MPI_Status istatus;
    int nitems;
    
    nTotObjs = info_->getNGlobalIntegrableObjects();
    haveError = 0;
    
    if (worldRank == masterNode) {
      inFile_->clear();            
      inFile_->seekg(framePos_[whichFrame]);
      
      if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
        sprintf(painCave.errMsg, "DumpReader Error: Error reading 1st line of %s \n ",
                filename_.c_str());
        painCave.isFatal = 1;
        simError();
      }
      
      nitems = atoi(read_buffer);
      
      // Check to see that the number of integrable objects in the
      // intial configuration file is the same as derived from the
      // meta-data file.
      
      if (nTotObjs != nitems) {
        sprintf(painCave.errMsg,
                "DumpReader Error. %s nIntegrable, %d, "
                "does not match the meta-data file's nIntegrable, %d.\n",
                filename_.c_str(),
                nTotObjs,
                info_->getNGlobalIntegrableObjects());
        
        painCave.isFatal = 1;
        simError();
      }
      
      //read the boxMat from the comment line
      
      
      if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
        sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
                filename_.c_str());
        painCave.isFatal = 1;
        simError();
      }
      
      //Every single processor will parse the comment line by itself
      //By using this way, we might lose some efficiency, but if we want to add
      //more parameters into comment line, we only need to modify function
      //parseCommentLine
      
      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
      
      for(i = 0; i < info_->getNGlobalMolecules(); i++) {
        int which_node = info_->getMolToProc(i);
        
        if (which_node == masterNode) {
          //molecules belong to master node
          
          mol = info_->getMoleculeByGlobalIndex(i);
          
          if (mol == NULL) {
            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
            painCave.isFatal = 1;
            simError();
          }
          
          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
               integrableObject = mol->nextIntegrableObject(ii)){
            
            
            if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
              sprintf(painCave.errMsg,
                      "DumpReader Error: error in reading file %s\n"
                      "natoms  = %d; index = %d\n"
                      "error reading the line from the file.\n",
                      filename_.c_str(),
                      nTotObjs,
                      i);
              
              painCave.isFatal = 1;
              simError();
            }
            
            parseDumpLine(read_buffer, integrableObject);
          }
        } else {
          //molecule belongs to slave nodes
          
          MPI_Recv(&nCurObj, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
                   MPI_COMM_WORLD, &istatus);
          
          for(int j = 0; j < nCurObj; j++) {
            
            
            if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
              sprintf(painCave.errMsg,
                      "DumpReader Error: error in reading file %s\n"
                      "natoms  = %d; index = %d\n"
                      "error reading the line from the file.\n",
                      filename_.c_str(),
                      nTotObjs,
                      i);
              
              painCave.isFatal = 1;
              simError();
            }
            
            MPI_Send(read_buffer, maxBufferSize, MPI_CHAR, which_node,
                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
          }
        }
      }
    } else {
      //actions taken at slave nodes
      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
      
      /**@todo*/
      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
      
      for(i = 0; i < info_->getNGlobalMolecules(); i++) {
        int which_node = info_->getMolToProc(i);
        
        if (which_node == worldRank) {
          //molecule with global index i belongs to this processor
          
          mol = info_->getMoleculeByGlobalIndex(i);
          if (mol == NULL) {
            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
            painCave.isFatal = 1;
            simError();
          }
          
          nCurObj = mol->getNIntegrableObjects();
          
          MPI_Send(&nCurObj, 1, MPI_INT, masterNode, TAKE_THIS_TAG_INT,
                   MPI_COMM_WORLD);
          
          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
               integrableObject = mol->nextIntegrableObject(ii)){
            
            MPI_Recv(read_buffer, maxBufferSize, MPI_CHAR, masterNode,
                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
            
            parseDumpLine(read_buffer, integrableObject);
          }
          
        }
        
      }
      
    }
    
#endif
    
  }
  
  void DumpReader::parseDumpLine(char *line, StuntDouble *integrableObject) {
    
    Vector3d pos;  // position place holders
    Vector3d vel;  // velocity placeholders
    Quat4d q;    // the quaternions
    Vector3d ji;   // angular velocity placeholders;
    StringTokenizer tokenizer(line);
    int nTokens;
    
    nTokens = tokenizer.countTokens();
    
    if (nTokens < 14) {
      sprintf(painCave.errMsg,
              "DumpReader Error: Not enough Tokens.\n%s\n", line);
      painCave.isFatal = 1;
      simError();
    }
    
    std::string name = tokenizer.nextToken();
    
    if (name != integrableObject->getType()) {
      
      sprintf(painCave.errMsg,
              "DumpReader Error: Atom type [%s] in %s does not match Atom Type [%s] in .md file.\n",
              name.c_str(), filename_.c_str(), integrableObject->getType().c_str());
      painCave.isFatal = 1;
      simError();        
    }
    
    pos[0] = tokenizer.nextTokenAsDouble();
    pos[1] = tokenizer.nextTokenAsDouble();
    pos[2] = tokenizer.nextTokenAsDouble();
    if (needPos_) {
      integrableObject->setPos(pos);
    }
    
    vel[0] = tokenizer.nextTokenAsDouble();
    vel[1] = tokenizer.nextTokenAsDouble();
    vel[2] = tokenizer.nextTokenAsDouble();
    if (needVel_) {
      integrableObject->setVel(vel);
    }
    
    if (integrableObject->isDirectional()) {
      
      q[0] = tokenizer.nextTokenAsDouble();
      q[1] = tokenizer.nextTokenAsDouble();
      q[2] = tokenizer.nextTokenAsDouble();
      q[3] = tokenizer.nextTokenAsDouble();
      
      double qlen = q.length();
      if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
        
        sprintf(painCave.errMsg,
                "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0\n");
        painCave.isFatal = 1;
        simError();
        
      } 
      
      q.normalize();
      if (needQuaternion_) {           
        integrableObject->setQ(q);
      }
      
      ji[0] = tokenizer.nextTokenAsDouble();
      ji[1] = tokenizer.nextTokenAsDouble();
      ji[2] = tokenizer.nextTokenAsDouble();
      if (needAngMom_) {
        integrableObject->setJ(ji);
      }
    }
    
  }
  
  
  void DumpReader::parseCommentLine(char* line, Snapshot* s) {
    double currTime;
    Mat3x3d hmat;
    double chi;
    double integralOfChiDt;
    Mat3x3d eta;
    
    StringTokenizer tokenizer(line);
    int nTokens;
    
    nTokens = tokenizer.countTokens();
    
    //comment line should at least contain 10 tokens: current time(1 token) and  h-matrix(9 tokens)
    if (nTokens < 10) {
      sprintf(painCave.errMsg,
              "DumpReader Error: Not enough tokens in comment line: %s", line);
      painCave.isFatal = 1;
      simError();   
    }
    
    //read current time
    currTime = tokenizer.nextTokenAsDouble();
    s->setTime(currTime);
    
    //read h-matrix
    hmat(0, 0) = tokenizer.nextTokenAsDouble();
    hmat(0, 1) = tokenizer.nextTokenAsDouble();
    hmat(0, 2) = tokenizer.nextTokenAsDouble();
    hmat(1, 0) = tokenizer.nextTokenAsDouble();
    hmat(1, 1) = tokenizer.nextTokenAsDouble();
    hmat(1, 2) = tokenizer.nextTokenAsDouble();
    hmat(2, 0) = tokenizer.nextTokenAsDouble();
    hmat(2, 1) = tokenizer.nextTokenAsDouble();
    hmat(2, 2) = tokenizer.nextTokenAsDouble();
    s->setHmat(hmat);
    
    //read chi and integralOfChidt, they should apprear in pair
    if (tokenizer.countTokens() >= 2) {
      chi = tokenizer.nextTokenAsDouble();
      integralOfChiDt = tokenizer.nextTokenAsDouble();            
      
      s->setChi(chi);
      s->setIntegralOfChiDt(integralOfChiDt);
    }
    
    //read eta (eta is 3x3 matrix)
    if (tokenizer.countTokens() >= 9) {
      eta(0, 0) = tokenizer.nextTokenAsDouble();
      eta(0, 1) = tokenizer.nextTokenAsDouble();
      eta(0, 2) = tokenizer.nextTokenAsDouble();
      eta(1, 0) = tokenizer.nextTokenAsDouble();
      eta(1, 1) = tokenizer.nextTokenAsDouble();
      eta(1, 2) = tokenizer.nextTokenAsDouble();
      eta(2, 0) = tokenizer.nextTokenAsDouble();
      eta(2, 1) = tokenizer.nextTokenAsDouble();
      eta(2, 2) = tokenizer.nextTokenAsDouble();      
      
      s->setEta(eta);
    }
    
    
  }
  
}//end namespace oopse
Revision:	2334
Committed:	Wed Sep 28 16:32:30 2005 UTC (18 years, 9 months ago) by gezelter
File size:	19397 byte(s)
Log Message:	tracking down quaternion bug
#	User	Rev	Content
1	tim	2318	/*
2			* Copyright (c) 2005 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
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40			*/
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 {
71
72			DumpReader::DumpReader(SimInfo* info, const std::string& filename)
73			: info_(info), filename_(filename), isScanned_(false), nframes_(0) {
74
75			#ifdef IS_MPI
76
77			if (worldRank == 0) {
78			#endif
79
80			inFile_ = new std::ifstream(filename_.c_str());
81
82			if (inFile_->fail()) {
83			sprintf(painCave.errMsg, "DumpReader: Cannot open file: %s\n", filename_.c_str());
84			painCave.isFatal = 1;
85			simError();
86			}
87
88			#ifdef IS_MPI
89
90			}
91
92			strcpy(checkPointMsg, "Dump file opened for reading successfully.");
93			MPIcheckPoint();
94
95			#endif
96
97			return;
98			}
99
100			DumpReader::~DumpReader() {
101
102			#ifdef IS_MPI
103
104			if (worldRank == 0) {
105			#endif
106
107			delete inFile_;
108
109			#ifdef IS_MPI
110
111			}
112
113			strcpy(checkPointMsg, "Dump file closed successfully.");
114			MPIcheckPoint();
115
116			#endif
117
118			return;
119			}
120
121			int DumpReader::getNFrames(void) {
122
123			if (!isScanned_)
124			scanFile();
125
126			return nframes_;
127			}
128
129			void DumpReader::scanFile(void) {
130			int i, j;
131			int lineNum = 0;
132			char readBuffer[maxBufferSize];
133			std::streampos currPos;
134
135			#ifdef IS_MPI
136
137			if (worldRank == 0) {
138			#endif // is_mpi
139
140			inFile_->seekg (0, std::ios::beg);
141
142
143			currPos = inFile_->tellg();
144			inFile_->getline(readBuffer, sizeof(readBuffer));
145			lineNum++;
146
147			if (inFile_->eof()) {
148			sprintf(painCave.errMsg,
149			"DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
150			filename_.c_str(),
151			lineNum);
152			painCave.isFatal = 1;
153			simError();
154			}
155
156			while (!inFile_->eof()) {
157			framePos_.push_back(currPos);
158
159			i = atoi(readBuffer);
160
161			inFile_->getline(readBuffer, sizeof(readBuffer));
162			lineNum++;
163
164			if (inFile_->eof()) {
165			sprintf(painCave.errMsg,
166			"DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
167			filename_.c_str(),
168			lineNum);
169			painCave.isFatal = 1;
170			simError();
171			}
172
173			for(j = 0; j < i; j++) {
174			inFile_->getline(readBuffer, sizeof(readBuffer));
175			lineNum++;
176
177			if (inFile_->eof()) {
178			sprintf(painCave.errMsg,
179			"DumpReader Error: File \"%s\" ended unexpectedly at line %d,"
180			" with atom %d\n", filename_.c_str(),
181			lineNum,
182			j);
183
184			painCave.isFatal = 1;
185			simError();
186			}
187			}
188
189			currPos = inFile_->tellg();
190			inFile_->getline(readBuffer, sizeof(readBuffer));
191			lineNum++;
192			}
193
194			inFile_->seekg (0, std::ios::beg);
195
196			nframes_ = framePos_.size();
197			#ifdef IS_MPI
198			}
199
200			MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
201
202			strcpy(checkPointMsg, "Successfully scanned DumpFile\n");
203			MPIcheckPoint();
204
205			#endif // is_mpi
206
207			isScanned_ = true;
208			}
209
210			void DumpReader::readFrame(int whichFrame) {
211			if (!isScanned_)
212			scanFile();
213
214			int storageLayout = info_->getSnapshotManager()->getStorageLayout();
215
216			if (storageLayout & DataStorage::dslPosition) {
217			needPos_ = true;
218			} else {
219			needPos_ = false;
220			}
221
222			if (storageLayout & DataStorage::dslVelocity) {
223			needVel_ = true;
224			} else {
225			needVel_ = false;
226			}
227
228			if (storageLayout & DataStorage::dslAmat \|\| storageLayout & DataStorage::dslElectroFrame) {
229			needQuaternion_ = true;
230			} else {
231			needQuaternion_ = false;
232			}
233
234			if (storageLayout & DataStorage::dslAngularMomentum) {
235			needAngMom_ = true;
236			} else {
237			needAngMom_ = false;
238			}
239
240			readSet(whichFrame);
241			}
242
243			void DumpReader::readSet(int whichFrame) {
244			int i;
245			int nTotObjs; // the number of atoms
246			char read_buffer[maxBufferSize]; //the line buffer for reading
247			char * eof_test; // ptr to see when we reach the end of the file
248
249			Molecule* mol;
250			StuntDouble* integrableObject;
251			SimInfo::MoleculeIterator mi;
252			Molecule::IntegrableObjectIterator ii;
253
254			#ifndef IS_MPI
255			inFile_->clear();
256			inFile_->seekg(framePos_[whichFrame]);
257
258			if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
259			sprintf(painCave.errMsg,
260			"DumpReader error: error reading 1st line of \"%s\"\n",
261			filename_.c_str());
262			painCave.isFatal = 1;
263			simError();
264			}
265
266			nTotObjs = atoi(read_buffer);
267
268			if (nTotObjs != info_->getNGlobalIntegrableObjects()) {
269			sprintf(painCave.errMsg,
270			"DumpReader error. %s nIntegrable, %d, "
271			"does not match the meta-data file's nIntegrable, %d.\n",
272			filename_.c_str(),
273			nTotObjs,
274			info_->getNGlobalIntegrableObjects());
275
276			painCave.isFatal = 1;
277			simError();
278			}
279
280			//read the box mat from the comment line
281
282
283			if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
284			sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
285			filename_.c_str());
286			painCave.isFatal = 1;
287			simError();
288			}
289
290			parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
291
292			//parse dump lines
293
294			i = 0;
295			for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
296
297			for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
298			integrableObject = mol->nextIntegrableObject(ii)) {
299
300
301
302			if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
303			sprintf(painCave.errMsg,
304			"DumpReader Error: error in reading file %s\n"
305			"natoms = %d; index = %d\n"
306			"error reading the line from the file.\n",
307			filename_.c_str(),
308			nTotObjs,
309			i);
310
311			painCave.isFatal = 1;
312			simError();
313			}
314
315			parseDumpLine(read_buffer, integrableObject);
316			i++;
317			}
318			}
319
320			// MPI Section of code..........
321
322			#else //IS_MPI
323
324			// first thing first, suspend fatalities.
325			int masterNode = 0;
326			int nCurObj;
327			painCave.isEventLoop = 1;
328
329			int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
330			int haveError;
331
332			MPI_Status istatus;
333			int nitems;
334
335			nTotObjs = info_->getNGlobalIntegrableObjects();
336			haveError = 0;
337
338			if (worldRank == masterNode) {
339			inFile_->clear();
340			inFile_->seekg(framePos_[whichFrame]);
341
342	tim	2327	if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
343	tim	2318	sprintf(painCave.errMsg, "DumpReader Error: Error reading 1st line of %s \n ",
344			filename_.c_str());
345			painCave.isFatal = 1;
346			simError();
347			}
348
349			nitems = atoi(read_buffer);
350
351			// Check to see that the number of integrable objects in the
352			// intial configuration file is the same as derived from the
353			// meta-data file.
354
355			if (nTotObjs != nitems) {
356			sprintf(painCave.errMsg,
357			"DumpReader Error. %s nIntegrable, %d, "
358			"does not match the meta-data file's nIntegrable, %d.\n",
359			filename_.c_str(),
360			nTotObjs,
361			info_->getNGlobalIntegrableObjects());
362
363			painCave.isFatal = 1;
364			simError();
365			}
366
367			//read the boxMat from the comment line
368
369
370
371	tim	2327	if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
372	tim	2318	sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
373			filename_.c_str());
374			painCave.isFatal = 1;
375			simError();
376			}
377
378			//Every single processor will parse the comment line by itself
379			//By using this way, we might lose some efficiency, but if we want to add
380			//more parameters into comment line, we only need to modify function
381			//parseCommentLine
382
383			MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
384			parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
385
386			for(i = 0; i < info_->getNGlobalMolecules(); i++) {
387			int which_node = info_->getMolToProc(i);
388
389			if (which_node == masterNode) {
390			//molecules belong to master node
391
392			mol = info_->getMoleculeByGlobalIndex(i);
393
394			if (mol == NULL) {
395			sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
396			painCave.isFatal = 1;
397			simError();
398			}
399
400			for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
401			integrableObject = mol->nextIntegrableObject(ii)){
402
403
404
405	tim	2327	if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
406	tim	2318	sprintf(painCave.errMsg,
407			"DumpReader Error: error in reading file %s\n"
408			"natoms = %d; index = %d\n"
409			"error reading the line from the file.\n",
410			filename_.c_str(),
411			nTotObjs,
412			i);
413
414			painCave.isFatal = 1;
415			simError();
416			}
417
418			parseDumpLine(read_buffer, integrableObject);
419			}
420			} else {
421			//molecule belongs to slave nodes
422
423			MPI_Recv(&nCurObj, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
424			MPI_COMM_WORLD, &istatus);
425
426			for(int j = 0; j < nCurObj; j++) {
427
428
429	tim	2327	if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
430	tim	2318	sprintf(painCave.errMsg,
431			"DumpReader Error: error in reading file %s\n"
432			"natoms = %d; index = %d\n"
433			"error reading the line from the file.\n",
434			filename_.c_str(),
435			nTotObjs,
436			i);
437
438			painCave.isFatal = 1;
439			simError();
440			}
441
442			MPI_Send(read_buffer, maxBufferSize, MPI_CHAR, which_node,
443			TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
444			}
445			}
446			}
447			} else {
448			//actions taken at slave nodes
449			MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
450
451			/*@todo/
452			parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
453
454			for(i = 0; i < info_->getNGlobalMolecules(); i++) {
455			int which_node = info_->getMolToProc(i);
456
457			if (which_node == worldRank) {
458			//molecule with global index i belongs to this processor
459
460			mol = info_->getMoleculeByGlobalIndex(i);
461			if (mol == NULL) {
462			sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
463			painCave.isFatal = 1;
464			simError();
465			}
466
467			nCurObj = mol->getNIntegrableObjects();
468
469			MPI_Send(&nCurObj, 1, MPI_INT, masterNode, TAKE_THIS_TAG_INT,
470			MPI_COMM_WORLD);
471
472			for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
473			integrableObject = mol->nextIntegrableObject(ii)){
474
475			MPI_Recv(read_buffer, maxBufferSize, MPI_CHAR, masterNode,
476			TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
477
478			parseDumpLine(read_buffer, integrableObject);
479			}
480
481			}
482
483			}
484
485			}
486
487			#endif
488
489			}
490
491			void DumpReader::parseDumpLine(char line, StuntDouble integrableObject) {
492
493			Vector3d pos; // position place holders
494			Vector3d vel; // velocity placeholders
495			Quat4d q; // the quaternions
496			Vector3d ji; // angular velocity placeholders;
497			StringTokenizer tokenizer(line);
498			int nTokens;
499
500			nTokens = tokenizer.countTokens();
501
502			if (nTokens < 14) {
503			sprintf(painCave.errMsg,
504			"DumpReader Error: Not enough Tokens.\n%s\n", line);
505			painCave.isFatal = 1;
506			simError();
507			}
508
509			std::string name = tokenizer.nextToken();
510
511			if (name != integrableObject->getType()) {
512
513			sprintf(painCave.errMsg,
514			"DumpReader Error: Atom type [%s] in %s does not match Atom Type [%s] in .md file.\n",
515			name.c_str(), filename_.c_str(), integrableObject->getType().c_str());
516			painCave.isFatal = 1;
517			simError();
518			}
519
520			pos[0] = tokenizer.nextTokenAsDouble();
521			pos[1] = tokenizer.nextTokenAsDouble();
522			pos[2] = tokenizer.nextTokenAsDouble();
523			if (needPos_) {
524			integrableObject->setPos(pos);
525			}
526
527			vel[0] = tokenizer.nextTokenAsDouble();
528			vel[1] = tokenizer.nextTokenAsDouble();
529			vel[2] = tokenizer.nextTokenAsDouble();
530			if (needVel_) {
531			integrableObject->setVel(vel);
532			}
533
534			if (integrableObject->isDirectional()) {
535
536			q[0] = tokenizer.nextTokenAsDouble();
537			q[1] = tokenizer.nextTokenAsDouble();
538			q[2] = tokenizer.nextTokenAsDouble();
539			q[3] = tokenizer.nextTokenAsDouble();
540
541			double qlen = q.length();
542			if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
543
544			sprintf(painCave.errMsg,
545	gezelter	2334	"DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0\n");
546	tim	2318	painCave.isFatal = 1;
547			simError();
548
549			}
550
551			q.normalize();
552			if (needQuaternion_) {
553			integrableObject->setQ(q);
554			}
555
556			ji[0] = tokenizer.nextTokenAsDouble();
557			ji[1] = tokenizer.nextTokenAsDouble();
558			ji[2] = tokenizer.nextTokenAsDouble();
559			if (needAngMom_) {
560			integrableObject->setJ(ji);
561			}
562			}
563
564			}
565
566
567			void DumpReader::parseCommentLine(char* line, Snapshot* s) {
568			double currTime;
569			Mat3x3d hmat;
570			double chi;
571			double integralOfChiDt;
572			Mat3x3d eta;
573
574			StringTokenizer tokenizer(line);
575			int nTokens;
576
577			nTokens = tokenizer.countTokens();
578
579			//comment line should at least contain 10 tokens: current time(1 token) and h-matrix(9 tokens)
580			if (nTokens < 10) {
581			sprintf(painCave.errMsg,
582			"DumpReader Error: Not enough tokens in comment line: %s", line);
583			painCave.isFatal = 1;
584			simError();
585			}
586
587			//read current time
588			currTime = tokenizer.nextTokenAsDouble();
589			s->setTime(currTime);
590
591			//read h-matrix
592			hmat(0, 0) = tokenizer.nextTokenAsDouble();
593			hmat(0, 1) = tokenizer.nextTokenAsDouble();
594			hmat(0, 2) = tokenizer.nextTokenAsDouble();
595			hmat(1, 0) = tokenizer.nextTokenAsDouble();
596			hmat(1, 1) = tokenizer.nextTokenAsDouble();
597			hmat(1, 2) = tokenizer.nextTokenAsDouble();
598			hmat(2, 0) = tokenizer.nextTokenAsDouble();
599			hmat(2, 1) = tokenizer.nextTokenAsDouble();
600			hmat(2, 2) = tokenizer.nextTokenAsDouble();
601			s->setHmat(hmat);
602
603			//read chi and integralOfChidt, they should apprear in pair
604			if (tokenizer.countTokens() >= 2) {
605			chi = tokenizer.nextTokenAsDouble();
606			integralOfChiDt = tokenizer.nextTokenAsDouble();
607
608			s->setChi(chi);
609			s->setIntegralOfChiDt(integralOfChiDt);
610			}
611
612			//read eta (eta is 3x3 matrix)
613			if (tokenizer.countTokens() >= 9) {
614			eta(0, 0) = tokenizer.nextTokenAsDouble();
615			eta(0, 1) = tokenizer.nextTokenAsDouble();
616			eta(0, 2) = tokenizer.nextTokenAsDouble();
617			eta(1, 0) = tokenizer.nextTokenAsDouble();
618			eta(1, 1) = tokenizer.nextTokenAsDouble();
619			eta(1, 2) = tokenizer.nextTokenAsDouble();
620			eta(2, 0) = tokenizer.nextTokenAsDouble();
621			eta(2, 1) = tokenizer.nextTokenAsDouble();
622			eta(2, 2) = tokenizer.nextTokenAsDouble();
623
624			s->setEta(eta);
625			}
626
627
628			}
629
630			}//end namespace oopse