src/io/DumpReader.cpp

/*
 * Copyright (c) 2009 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. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. 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.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */
  
#define _LARGEFILE_SOURCE64 
#define _FILE_OFFSET_BITS 64 
 
#ifdef IS_MPI 
#include <mpi.h> 
#endif
 
#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" 
#include "brains/Thermo.hpp"
 
 
namespace OpenMD { 
   
  DumpReader::DumpReader(SimInfo* info, const std::string& filename) 
    : info_(info), filename_(filename), isScanned_(false), nframes_(0), needCOMprops_(false) { 
    
#ifdef IS_MPI 
    
    if (worldRank == 0) { 
#endif 
      
      inFile_ = new std::ifstream(filename_.c_str(),   
                                  ifstream::in | ifstream::binary); 
      
      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."); 
    errorCheckPoint(); 
    
#endif 
    
    return; 
  } 
  
  DumpReader::~DumpReader() { 
    
#ifdef IS_MPI 
    
    if (worldRank == 0) { 
#endif

      delete inFile_; 
      
#ifdef IS_MPI 
      
    } 
    
    strcpy(checkPointMsg, "Dump file closed successfully."); 
    errorCheckPoint(); 
    
#endif 
    
    return; 
  } 
  
  int DumpReader::getNFrames(void) { 
     
    if (!isScanned_) 
      scanFile(); 
     
    return nframes_; 
  } 
   
  void DumpReader::scanFile(void) { 

    std::streampos prevPos;
    std::streampos  currPos; 
    
#ifdef IS_MPI 
    
    if (worldRank == 0) { 
#endif // is_mpi 
      
      currPos = inFile_->tellg();
      prevPos = currPos;
      bool foundOpenSnapshotTag = false;
      bool foundClosedSnapshotTag = false;

      int lineNo = 0; 
      while(inFile_->getline(buffer, bufferSize)) {
        ++lineNo;
        
        std::string line = buffer;
        currPos = inFile_->tellg(); 
        if (line.find("<Snapshot>")!= std::string::npos) {
          if (foundOpenSnapshotTag) {
            sprintf(painCave.errMsg, 
                    "DumpReader:<Snapshot> is multiply nested at line %d in %s \n", lineNo, 
                    filename_.c_str()); 
            painCave.isFatal = 1; 
            simError();           
          }
          foundOpenSnapshotTag = true;
          foundClosedSnapshotTag = false;
          framePos_.push_back(prevPos);
          
        } else if (line.find("</Snapshot>") != std::string::npos){
          if (!foundOpenSnapshotTag) {
            sprintf(painCave.errMsg, 
                    "DumpReader:</Snapshot> appears before <Snapshot> at line %d in %s \n", lineNo, 
                    filename_.c_str()); 
            painCave.isFatal = 1; 
            simError(); 
          }
          
          if (foundClosedSnapshotTag) {
            sprintf(painCave.errMsg, 
                    "DumpReader:</Snapshot> appears multiply nested at line %d in %s \n", lineNo, 
                    filename_.c_str()); 
            painCave.isFatal = 1; 
            simError(); 
          }
          foundClosedSnapshotTag = true;
          foundOpenSnapshotTag = false;
        }
        prevPos = currPos;
      }
      
      // only found <Snapshot> for the last frame means the file is corrupted, we should discard
      // it and give a warning message
      if (foundOpenSnapshotTag) {
        sprintf(painCave.errMsg, 
                "DumpReader: last frame in %s is invalid\n", filename_.c_str()); 
        painCave.isFatal = 0; 
        simError();       
        framePos_.pop_back();
      }
      
      nframes_ = framePos_.size(); 
      
      if (nframes_ == 0) {
        sprintf(painCave.errMsg, 
                "DumpReader: %s does not contain a valid frame\n", filename_.c_str()); 
        painCave.isFatal = 1; 
        simError();      
      }
#ifdef IS_MPI 
    } 
     
    MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD); 
    
#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::dslDipole || 
        storageLayout & DataStorage::dslQuadrupole) { 
      needQuaternion_ = true; 
    } else { 
      needQuaternion_ = false; 
    } 
     
    if (storageLayout & DataStorage::dslAngularMomentum) { 
      needAngMom_ = true; 
    } else { 
      needAngMom_ = false;     
    } 
     
    readSet(whichFrame); 

    if (needCOMprops_) {
      Thermo thermo(info_);
      Vector3d com;

      if (needPos_ && needVel_) {
        Vector3d comvel;
        Vector3d comw;
        thermo.getComAll(com, comvel);
        comw = thermo.getAngularMomentum();
      } else {
        com = thermo.getCom();
      }                    
    }
  } 
   
  void DumpReader::readSet(int whichFrame) {     
    std::string line;

#ifndef IS_MPI 
    inFile_->clear();  
    inFile_->seekg(framePos_[whichFrame]); 

    std::istream& inputStream = *inFile_;     

#else 
    int masterNode = 0;
    std::stringstream sstream;
    if (worldRank == masterNode) {
      std::string sendBuffer;

      inFile_->clear();  
      inFile_->seekg(framePos_[whichFrame]); 
      
      while (inFile_->getline(buffer, bufferSize)) {

        line = buffer;
        sendBuffer += line;
        sendBuffer += '\n';
        if (line.find("</Snapshot>") != std::string::npos) {
          break;
        }        
      }

      int sendBufferSize = sendBuffer.size();
      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);     
      MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, 
                MPI_CHAR, masterNode, MPI_COMM_WORLD);     
      
      sstream.str(sendBuffer);
    } else {
      int sendBufferSize;
      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
      char * recvBuffer = new char[sendBufferSize+1];
      assert(recvBuffer);
      recvBuffer[sendBufferSize] = '\0';
      MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
      sstream.str(recvBuffer);
      delete [] recvBuffer;
    }      

    std::istream& inputStream = sstream;  
#endif

    inputStream.getline(buffer, bufferSize);

    line = buffer;
    if (line.find("<Snapshot>") == std::string::npos) {
      sprintf(painCave.errMsg, 
              "DumpReader Error: can not find <Snapshot>\n"); 
      painCave.isFatal = 1; 
      simError(); 
    } 
    
    //read frameData
    readFrameProperties(inputStream);

    //read StuntDoubles
    readStuntDoubles(inputStream);     

    inputStream.getline(buffer, bufferSize);
    line = buffer;

    if (line.find("<SiteData>") != std::string::npos) {
      //read SiteData
      readSiteData(inputStream);         
    } else {
      if (line.find("</Snapshot>") == std::string::npos) {
        sprintf(painCave.errMsg, 
                "DumpReader Error: can not find </Snapshot>\n"); 
        painCave.isFatal = 1; 
        simError(); 
      }        
    }
  } 
   
  void DumpReader::parseDumpLine(const std::string& line) { 

       
    StringTokenizer tokenizer(line); 
    int nTokens; 
     
    nTokens = tokenizer.countTokens(); 
     
    if (nTokens < 2) {  
      sprintf(painCave.errMsg, 
              "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str()); 
      painCave.isFatal = 1; 
      simError(); 
    } 

    int index = tokenizer.nextTokenAsInt();
 
    StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);

    if (sd == NULL) {
      return;
    }
    std::string type = tokenizer.nextToken(); 
    int size = type.size();

    size_t found;
    
    if (needPos_) {
      found = type.find("p");      
      if (found == std::string::npos) {
        sprintf(painCave.errMsg, 
                "DumpReader Error: StuntDouble %d has no Position\n"
                "\tField (\"p\") specified.\n%s\n", index, 
                line.c_str());  
        painCave.isFatal = 1; 
        simError(); 
      }
    }
    
    if (sd->isDirectional()) {
      if (needQuaternion_) {
        found = type.find("q");      
        if (found == std::string::npos) {
          sprintf(painCave.errMsg, 
                  "DumpReader Error: Directional StuntDouble %d has no\n"
                  "\tQuaternion Field (\"q\") specified.\n%s\n", index, 
                  line.c_str());  
          painCave.isFatal = 1; 
          simError(); 
        }
      }      
    }

    for(int i = 0; i < size; ++i) {
      switch(type[i]) {
        
        case 'p': {
            Vector3d pos;
            pos[0] = tokenizer.nextTokenAsDouble(); 
            pos[1] = tokenizer.nextTokenAsDouble(); 
            pos[2] = tokenizer.nextTokenAsDouble(); 
            if (needPos_) { 
              sd->setPos(pos); 
            }             
            break;
        }
        case 'v' : {
            Vector3d vel;
            vel[0] = tokenizer.nextTokenAsDouble(); 
            vel[1] = tokenizer.nextTokenAsDouble(); 
            vel[2] = tokenizer.nextTokenAsDouble(); 
            if (needVel_) { 
              sd->setVel(vel); 
            } 
            break;
        }

        case 'q' : {
           Quat4d q;
           if (sd->isDirectional()) { 
              
             q[0] = tokenizer.nextTokenAsDouble(); 
             q[1] = tokenizer.nextTokenAsDouble(); 
             q[2] = tokenizer.nextTokenAsDouble(); 
             q[3] = tokenizer.nextTokenAsDouble(); 
              
             RealType qlen = q.length(); 
             if (qlen < OpenMD::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_) {            
               sd->setQ(q); 
             }               
           }            
           break;
        }  
        case 'j' : {
          Vector3d ji;
          if (sd->isDirectional()) {
             ji[0] = tokenizer.nextTokenAsDouble(); 
             ji[1] = tokenizer.nextTokenAsDouble(); 
             ji[2] = tokenizer.nextTokenAsDouble(); 
             if (needAngMom_) { 
               sd->setJ(ji); 
             } 
          }
          break;
        }  
        case 'f': {

          Vector3d force;
          force[0] = tokenizer.nextTokenAsDouble(); 
          force[1] = tokenizer.nextTokenAsDouble(); 
          force[2] = tokenizer.nextTokenAsDouble();           
          sd->setFrc(force); 
          break;
        }
        case 't' : {

           Vector3d torque;
           torque[0] = tokenizer.nextTokenAsDouble(); 
           torque[1] = tokenizer.nextTokenAsDouble(); 
           torque[2] = tokenizer.nextTokenAsDouble();           
           sd->setTrq(torque);          
           break;
        }
        case 'u' : {

           RealType particlePot;
           particlePot = tokenizer.nextTokenAsDouble(); 
           sd->setParticlePot(particlePot);          
           break;
        }
        case 'c' : {

           RealType flucQPos;
           flucQPos = tokenizer.nextTokenAsDouble(); 
           sd->setFlucQPos(flucQPos);          
           break;
        }
        case 'w' : {

           RealType flucQVel;
           flucQVel = tokenizer.nextTokenAsDouble(); 
           sd->setFlucQVel(flucQVel);          
           break;
        }
        case 'g' : {

           RealType flucQFrc;
           flucQFrc = tokenizer.nextTokenAsDouble(); 
           sd->setFlucQFrc(flucQFrc);          
           break;
        }
        case 'e' : {

           Vector3d eField;
           eField[0] = tokenizer.nextTokenAsDouble(); 
           eField[1] = tokenizer.nextTokenAsDouble(); 
           eField[2] = tokenizer.nextTokenAsDouble();           
           sd->setElectricField(eField);          
           break;
        }
        case 's' : {

           RealType sPot;
           sPot = tokenizer.nextTokenAsDouble(); 
           sd->setSitePotential(sPot);          
           break;
        }
        default: {
               sprintf(painCave.errMsg, 
                       "DumpReader Error: %s is an unrecognized type\n", type.c_str()); 
               painCave.isFatal = 1; 
               simError(); 
          break;   
        }

      }
    }
  } 
   

  void DumpReader::parseSiteLine(const std::string& line) { 

    StringTokenizer tokenizer(line); 
    int nTokens; 
         
    nTokens = tokenizer.countTokens(); 
     
    if (nTokens < 2) {  
      sprintf(painCave.errMsg, 
              "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str()); 
      painCave.isFatal = 1; 
      simError(); 
    } 

    /**
     * The first token is the global integrable object index.
     */

    int index = tokenizer.nextTokenAsInt();
    StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
    if (sd == NULL) {
      return;
    }

    /**
     * Test to see if the next token is an integer or not.  If not,
     * we've got data on the integrable object itself.  If there is an
     * integer, we're parsing data for a site on a rigid body.
     */

    std::string indexTest = tokenizer.peekNextToken();
    std::istringstream i(indexTest);
    int siteIndex;
    if (i >> siteIndex) {
      // chew up this token and parse as an int:
      siteIndex = tokenizer.nextTokenAsInt();

      if (sd->isRigidBody()) {
        RigidBody* rb = static_cast<RigidBody*>(sd);
        sd = rb->getAtoms()[siteIndex];
      }
    }

    /**
     * The next token contains information on what follows.
     */
    std::string type = tokenizer.nextToken(); 
    int size = type.size();
    
    for(int i = 0; i < size; ++i) {
      switch(type[i]) {
        
      case 'u' : {
        
        RealType particlePot;
        particlePot = tokenizer.nextTokenAsDouble(); 
        sd->setParticlePot(particlePot);
        break;
      }
      case 'c' : {
        
        RealType flucQPos;
        flucQPos = tokenizer.nextTokenAsDouble(); 
        sd->setFlucQPos(flucQPos);
        break;
      }
      case 'w' : {
        
        RealType flucQVel;
        flucQVel = tokenizer.nextTokenAsDouble(); 
        sd->setFlucQVel(flucQVel);
        break;
      }
      case 'g' : {
        
        RealType flucQFrc;
        flucQFrc = tokenizer.nextTokenAsDouble(); 
        sd->setFlucQFrc(flucQFrc);
        break;
      }
      case 'e' : {
        
        Vector3d eField;
        eField[0] = tokenizer.nextTokenAsDouble(); 
        eField[1] = tokenizer.nextTokenAsDouble(); 
        eField[2] = tokenizer.nextTokenAsDouble();  
        sd->setElectricField(eField);          
        break;
      }
      case 's' : {
        
        RealType sPot;
        sPot = tokenizer.nextTokenAsDouble(); 
        sd->setSitePotential(sPot);          
        break;
      }
      default: {
        sprintf(painCave.errMsg, 
                "DumpReader Error: %s is an unrecognized type\n", type.c_str()); 
        painCave.isFatal = 1; 
        simError(); 
        break;   
      }
      }
    }    
  } 
  
  
    void  DumpReader::readStuntDoubles(std::istream& inputStream) {
    
    inputStream.getline(buffer, bufferSize);
    std::string line(buffer);
    
    if (line.find("<StuntDoubles>") == std::string::npos) {
      sprintf(painCave.errMsg, 
              "DumpReader Error: Missing <StuntDoubles>\n"); 
      painCave.isFatal = 1; 
      simError(); 
    }

    while(inputStream.getline(buffer, bufferSize)) {
      line = buffer;
      
      if(line.find("</StuntDoubles>") != std::string::npos) {
        break;
      }

      parseDumpLine(line);
    }
  
  }

  void  DumpReader::readSiteData(std::istream& inputStream) {

    std::string line(buffer);

    // We already found the starting <SiteData> tag or we wouldn't be
    // here, so just start parsing until we get to the ending
    // </SiteData> tag:
    
    while(inputStream.getline(buffer, bufferSize)) {
      line = buffer;
      
      if(line.find("</SiteData>") != std::string::npos) {
        break;
      }

      parseSiteLine(line);
    }
  
  }

  void DumpReader::readFrameProperties(std::istream& inputStream) {

    Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
    inputStream.getline(buffer, bufferSize);
    std::string line(buffer);

    if (line.find("<FrameData>") == std::string::npos) {
      sprintf(painCave.errMsg, 
              "DumpReader Error: Missing <FrameData>\n"); 
      painCave.isFatal = 1; 
      simError(); 
    }

    while(inputStream.getline(buffer, bufferSize)) {
      line = buffer;
      
      if(line.find("</FrameData>") != std::string::npos) {
        break;
      }
      
      StringTokenizer tokenizer(line, " ;\t\n\r{}:,");
      if (!tokenizer.hasMoreTokens()) {
        sprintf(painCave.errMsg, 
                "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str()); 
        painCave.isFatal = 1; 
        simError();      
      }

      std::string propertyName = tokenizer.nextToken();
      if (propertyName == "Time") {
        RealType currTime = tokenizer.nextTokenAsDouble(); 
        s->setTime(currTime); 
      } else if (propertyName == "Hmat"){
        Mat3x3d hmat;
        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);      
      } else if (propertyName == "Thermostat") {
        pair<RealType, RealType> thermostat;
        thermostat.first = tokenizer.nextTokenAsDouble();
        thermostat.second = tokenizer.nextTokenAsDouble();
        s->setThermostat(thermostat); 
     } else if (propertyName == "Barostat") {
        Mat3x3d eta;
        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->setBarostat(eta); 
      } else {
        sprintf(painCave.errMsg, 
                "DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str()); 
        painCave.isFatal = 0; 
        simError();         
      }
      
    }

  }

   
}//end namespace OpenMD 
Revision:	2069
Committed:	Thu Mar 5 16:30:23 2015 UTC (10 years, 1 month ago) by gezelter
File size:	21693 byte(s)
Log Message:	More g++ warning fixes
#	User	Rev	Content
1	gezelter	1390	/*
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	gezelter	1879	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	gezelter	1782	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	gezelter	1390	*/
42	chrisfen	721
43			#define _LARGEFILE_SOURCE64
44			#define _FILE_OFFSET_BITS 64
45
46	gezelter	1938	#ifdef IS_MPI
47			#include <mpi.h>
48			#endif
49
50	chrisfen	721	#include <sys/types.h>
51			#include <sys/stat.h>
52
53			#include <iostream>
54			#include <math.h>
55
56			#include <stdio.h>
57			#include <stdlib.h>
58			#include <string.h>
59
60			#include "io/DumpReader.hpp"
61			#include "primitives/Molecule.hpp"
62			#include "utils/simError.h"
63			#include "utils/MemoryUtils.hpp"
64			#include "utils/StringTokenizer.hpp"
65	gezelter	1782	#include "brains/Thermo.hpp"
66	chrisfen	721
67
68	gezelter	1390	namespace OpenMD {
69	chrisfen	721
70			DumpReader::DumpReader(SimInfo* info, const std::string& filename)
71	gezelter	1104	: info_(info), filename_(filename), isScanned_(false), nframes_(0), needCOMprops_(false) {
72	chrisfen	996
73	chrisfen	721	#ifdef IS_MPI
74	chrisfen	996
75			if (worldRank == 0) {
76	chrisfen	721	#endif
77	chrisfen	996
78	gezelter	1790	inFile_ = new std::ifstream(filename_.c_str(),
79			ifstream::in \| ifstream::binary);
80	chrisfen	996
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	chrisfen	721	#ifdef IS_MPI
90	chrisfen	996
91			}
92
93			strcpy(checkPointMsg, "Dump file opened for reading successfully.");
94	gezelter	1241	errorCheckPoint();
95	chrisfen	996
96	chrisfen	721	#endif
97	chrisfen	996
98			return;
99			}
100
101	chrisfen	721	DumpReader::~DumpReader() {
102	chrisfen	996
103	chrisfen	721	#ifdef IS_MPI
104	chrisfen	996
105	chrisfen	721	if (worldRank == 0) {
106	gezelter	2000	#endif
107
108	gezelter	2001	delete inFile_;
109	chrisfen	996
110	chrisfen	721	#ifdef IS_MPI
111	chrisfen	996
112	chrisfen	721	}
113	chrisfen	996
114	chrisfen	721	strcpy(checkPointMsg, "Dump file closed successfully.");
115	gezelter	1241	errorCheckPoint();
116	chrisfen	996
117	chrisfen	721	#endif
118	chrisfen	996
119	chrisfen	721	return;
120			}
121	chrisfen	996
122	chrisfen	721	int DumpReader::getNFrames(void) {
123
124			if (!isScanned_)
125			scanFile();
126
127			return nframes_;
128			}
129
130			void DumpReader::scanFile(void) {
131	gezelter	1879
132	tim	1024	std::streampos prevPos;
133	chrisfen	721	std::streampos currPos;
134	tim	1024
135	chrisfen	721	#ifdef IS_MPI
136	tim	1024
137	chrisfen	721	if (worldRank == 0) {
138			#endif // is_mpi
139	tim	1024
140			currPos = inFile_->tellg();
141			prevPos = currPos;
142			bool foundOpenSnapshotTag = false;
143			bool foundClosedSnapshotTag = false;
144	gezelter	1793
145	gezelter	1879	int lineNo = 0;
146	tim	1024	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	chrisfen	721	sprintf(painCave.errMsg,
154	tim	1024	"DumpReader:<Snapshot> is multiply nested at line %d in %s \n", lineNo,
155			filename_.c_str());
156	chrisfen	721	painCave.isFatal = 1;
157	tim	1024	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	chrisfen	721	simError();
170	tim	1024	}
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	chrisfen	721	nframes_ = framePos_.size();
196	tim	1024
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	chrisfen	721	#ifdef IS_MPI
204			}
205
206	gezelter	1969	MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
207	tim	1024
208	chrisfen	721	#endif // is_mpi
209	tim	1024
210	chrisfen	721	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	gezelter	1879	if (storageLayout & DataStorage::dslAmat \|\|
232			storageLayout & DataStorage::dslDipole \|\|
233			storageLayout & DataStorage::dslQuadrupole) {
234	chrisfen	721	needQuaternion_ = true;
235			} else {
236			needQuaternion_ = false;
237			}
238
239			if (storageLayout & DataStorage::dslAngularMomentum) {
240			needAngMom_ = true;
241			} else {
242			needAngMom_ = false;
243			}
244
245			readSet(whichFrame);
246	gezelter	1104
247			if (needCOMprops_) {
248	gezelter	1782	Thermo thermo(info_);
249	gezelter	1104	Vector3d com;
250	gezelter	1782
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	gezelter	1104	}
260	chrisfen	721	}
261
262	tim	1024	void DumpReader::readSet(int whichFrame) {
263			std::string line;
264
265	chrisfen	721	#ifndef IS_MPI
266			inFile_->clear();
267			inFile_->seekg(framePos_[whichFrame]);
268	tim	1024
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	gezelter	1969	MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
292			MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize,
293			MPI_CHAR, masterNode, MPI_COMM_WORLD);
294	tim	1024
295			sstream.str(sendBuffer);
296			} else {
297			int sendBufferSize;
298	gezelter	1969	MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
299	tim	1024	char * recvBuffer = new char[sendBufferSize+1];
300	gezelter	1313	assert(recvBuffer);
301			recvBuffer[sendBufferSize] = '\0';
302	gezelter	1969	MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
303	tim	1024	sstream.str(recvBuffer);
304	gezelter	1313	delete [] recvBuffer;
305	tim	1024	}
306
307			std::istream& inputStream = sstream;
308			#endif
309
310			inputStream.getline(buffer, bufferSize);
311
312			line = buffer;
313			if (line.find("<Snapshot>") == std::string::npos) {
314	chrisfen	721	sprintf(painCave.errMsg,
315	tim	1024	"DumpReader Error: can not find <Snapshot>\n");
316	chrisfen	721	painCave.isFatal = 1;
317			simError();
318			}
319	tim	1024
320			//read frameData
321			readFrameProperties(inputStream);
322
323			//read StuntDoubles
324			readStuntDoubles(inputStream);
325
326			inputStream.getline(buffer, bufferSize);
327			line = buffer;
328	gezelter	1782
329			if (line.find("<SiteData>") != std::string::npos) {
330			//read SiteData
331			readSiteData(inputStream);
332			} else {
333			if (line.find("</Snapshot>") == std::string::npos) {
334			sprintf(painCave.errMsg,
335			"DumpReader Error: can not find </Snapshot>\n");
336			painCave.isFatal = 1;
337			simError();
338			}
339			}
340	chrisfen	721	}
341
342	tim	1024	void DumpReader::parseDumpLine(const std::string& line) {
343
344
345	chrisfen	721	StringTokenizer tokenizer(line);
346			int nTokens;
347
348			nTokens = tokenizer.countTokens();
349
350	gezelter	1450	if (nTokens < 2) {
351	chrisfen	721	sprintf(painCave.errMsg,
352	tim	1024	"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
353	chrisfen	721	painCave.isFatal = 1;
354			simError();
355			}
356	tim	1024
357			int index = tokenizer.nextTokenAsInt();
358
359	gezelter	1782	StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
360	tim	1024
361	gezelter	1782	if (sd == NULL) {
362	tim	1024	return;
363			}
364			std::string type = tokenizer.nextToken();
365			int size = type.size();
366	gezelter	1037
367	gezelter	1450	size_t found;
368
369			if (needPos_) {
370			found = type.find("p");
371			if (found == std::string::npos) {
372			sprintf(painCave.errMsg,
373			"DumpReader Error: StuntDouble %d has no Position\n"
374			"\tField (\"p\") specified.\n%s\n", index,
375			line.c_str());
376			painCave.isFatal = 1;
377			simError();
378			}
379			}
380
381	gezelter	1782	if (sd->isDirectional()) {
382	gezelter	1450	if (needQuaternion_) {
383			found = type.find("q");
384			if (found == std::string::npos) {
385			sprintf(painCave.errMsg,
386			"DumpReader Error: Directional StuntDouble %d has no\n"
387			"\tQuaternion Field (\"q\") specified.\n%s\n", index,
388			line.c_str());
389			painCave.isFatal = 1;
390			simError();
391			}
392			}
393			}
394
395	tim	1024	for(int i = 0; i < size; ++i) {
396			switch(type[i]) {
397
398			case 'p': {
399			Vector3d pos;
400			pos[0] = tokenizer.nextTokenAsDouble();
401			pos[1] = tokenizer.nextTokenAsDouble();
402			pos[2] = tokenizer.nextTokenAsDouble();
403			if (needPos_) {
404	gezelter	1782	sd->setPos(pos);
405	tim	1024	}
406			break;
407			}
408			case 'v' : {
409			Vector3d vel;
410			vel[0] = tokenizer.nextTokenAsDouble();
411			vel[1] = tokenizer.nextTokenAsDouble();
412			vel[2] = tokenizer.nextTokenAsDouble();
413			if (needVel_) {
414	gezelter	1782	sd->setVel(vel);
415	tim	1024	}
416			break;
417			}
418
419			case 'q' : {
420			Quat4d q;
421	gezelter	1782	if (sd->isDirectional()) {
422	tim	1024
423			q[0] = tokenizer.nextTokenAsDouble();
424			q[1] = tokenizer.nextTokenAsDouble();
425			q[2] = tokenizer.nextTokenAsDouble();
426			q[3] = tokenizer.nextTokenAsDouble();
427
428			RealType qlen = q.length();
429	gezelter	1390	if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
430	tim	1024
431			sprintf(painCave.errMsg,
432			"DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
433			painCave.isFatal = 1;
434			simError();
435
436			}
437
438			q.normalize();
439			if (needQuaternion_) {
440	gezelter	1782	sd->setQ(q);
441	tim	1024	}
442			}
443			break;
444			}
445			case 'j' : {
446			Vector3d ji;
447	gezelter	1782	if (sd->isDirectional()) {
448	tim	1024	ji[0] = tokenizer.nextTokenAsDouble();
449			ji[1] = tokenizer.nextTokenAsDouble();
450			ji[2] = tokenizer.nextTokenAsDouble();
451			if (needAngMom_) {
452	gezelter	1782	sd->setJ(ji);
453	tim	1024	}
454			}
455	gezelter	1037	break;
456	tim	1024	}
457			case 'f': {
458
459			Vector3d force;
460			force[0] = tokenizer.nextTokenAsDouble();
461			force[1] = tokenizer.nextTokenAsDouble();
462			force[2] = tokenizer.nextTokenAsDouble();
463	gezelter	1782	sd->setFrc(force);
464	tim	1024	break;
465			}
466			case 't' : {
467
468			Vector3d torque;
469			torque[0] = tokenizer.nextTokenAsDouble();
470			torque[1] = tokenizer.nextTokenAsDouble();
471			torque[2] = tokenizer.nextTokenAsDouble();
472	gezelter	1782	sd->setTrq(torque);
473	tim	1024	break;
474			}
475	gezelter	1610	case 'u' : {
476
477			RealType particlePot;
478			particlePot = tokenizer.nextTokenAsDouble();
479	gezelter	1782	sd->setParticlePot(particlePot);
480	gezelter	1610	break;
481			}
482	gezelter	1782	case 'c' : {
483
484			RealType flucQPos;
485			flucQPos = tokenizer.nextTokenAsDouble();
486			sd->setFlucQPos(flucQPos);
487			break;
488			}
489			case 'w' : {
490
491			RealType flucQVel;
492			flucQVel = tokenizer.nextTokenAsDouble();
493			sd->setFlucQVel(flucQVel);
494			break;
495			}
496			case 'g' : {
497
498			RealType flucQFrc;
499			flucQFrc = tokenizer.nextTokenAsDouble();
500			sd->setFlucQFrc(flucQFrc);
501			break;
502			}
503			case 'e' : {
504
505			Vector3d eField;
506			eField[0] = tokenizer.nextTokenAsDouble();
507			eField[1] = tokenizer.nextTokenAsDouble();
508			eField[2] = tokenizer.nextTokenAsDouble();
509			sd->setElectricField(eField);
510			break;
511			}
512	gezelter	1993	case 's' : {
513
514			RealType sPot;
515			sPot = tokenizer.nextTokenAsDouble();
516			sd->setSitePotential(sPot);
517			break;
518			}
519	tim	1024	default: {
520			sprintf(painCave.errMsg,
521			"DumpReader Error: %s is an unrecognized type\n", type.c_str());
522			painCave.isFatal = 1;
523			simError();
524			break;
525			}
526
527			}
528			}
529			}
530
531
532	gezelter	1782	void DumpReader::parseSiteLine(const std::string& line) {
533
534			StringTokenizer tokenizer(line);
535			int nTokens;
536	gezelter	1879
537	gezelter	1782	nTokens = tokenizer.countTokens();
538
539			if (nTokens < 2) {
540			sprintf(painCave.errMsg,
541			"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
542			painCave.isFatal = 1;
543			simError();
544			}
545
546			/**
547			* The first token is the global integrable object index.
548			*/
549
550			int index = tokenizer.nextTokenAsInt();
551			StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
552			if (sd == NULL) {
553			return;
554			}
555
556			/**
557			* Test to see if the next token is an integer or not. If not,
558			* we've got data on the integrable object itself. If there is an
559			* integer, we're parsing data for a site on a rigid body.
560			*/
561
562			std::string indexTest = tokenizer.peekNextToken();
563			std::istringstream i(indexTest);
564			int siteIndex;
565			if (i >> siteIndex) {
566			// chew up this token and parse as an int:
567			siteIndex = tokenizer.nextTokenAsInt();
568	gezelter	1897
569			if (sd->isRigidBody()) {
570			RigidBody* rb = static_cast<RigidBody*>(sd);
571			sd = rb->getAtoms()[siteIndex];
572			}
573	gezelter	1782	}
574
575			/**
576			* The next token contains information on what follows.
577			*/
578			std::string type = tokenizer.nextToken();
579			int size = type.size();
580
581			for(int i = 0; i < size; ++i) {
582			switch(type[i]) {
583
584			case 'u' : {
585
586			RealType particlePot;
587			particlePot = tokenizer.nextTokenAsDouble();
588			sd->setParticlePot(particlePot);
589			break;
590			}
591			case 'c' : {
592
593			RealType flucQPos;
594			flucQPos = tokenizer.nextTokenAsDouble();
595			sd->setFlucQPos(flucQPos);
596			break;
597			}
598			case 'w' : {
599
600			RealType flucQVel;
601			flucQVel = tokenizer.nextTokenAsDouble();
602			sd->setFlucQVel(flucQVel);
603			break;
604			}
605			case 'g' : {
606
607			RealType flucQFrc;
608			flucQFrc = tokenizer.nextTokenAsDouble();
609			sd->setFlucQFrc(flucQFrc);
610			break;
611			}
612			case 'e' : {
613
614			Vector3d eField;
615			eField[0] = tokenizer.nextTokenAsDouble();
616			eField[1] = tokenizer.nextTokenAsDouble();
617			eField[2] = tokenizer.nextTokenAsDouble();
618			sd->setElectricField(eField);
619			break;
620			}
621	gezelter	1993	case 's' : {
622
623			RealType sPot;
624			sPot = tokenizer.nextTokenAsDouble();
625			sd->setSitePotential(sPot);
626			break;
627			}
628	gezelter	1782	default: {
629			sprintf(painCave.errMsg,
630			"DumpReader Error: %s is an unrecognized type\n", type.c_str());
631			painCave.isFatal = 1;
632			simError();
633			break;
634			}
635			}
636			}
637			}
638
639
640	gezelter	1993	void DumpReader::readStuntDoubles(std::istream& inputStream) {
641	gezelter	1782
642	tim	1024	inputStream.getline(buffer, bufferSize);
643			std::string line(buffer);
644
645			if (line.find("<StuntDoubles>") == std::string::npos) {
646	chrisfen	721	sprintf(painCave.errMsg,
647	tim	1024	"DumpReader Error: Missing <StuntDoubles>\n");
648	chrisfen	721	painCave.isFatal = 1;
649	tim	1024	simError();
650			}
651
652			while(inputStream.getline(buffer, bufferSize)) {
653			line = buffer;
654
655			if(line.find("</StuntDoubles>") != std::string::npos) {
656			break;
657			}
658
659			parseDumpLine(line);
660			}
661
662			}
663
664	gezelter	1782	void DumpReader::readSiteData(std::istream& inputStream) {
665
666			std::string line(buffer);
667	gezelter	1879
668			// We already found the starting <SiteData> tag or we wouldn't be
669			// here, so just start parsing until we get to the ending
670			// </SiteData> tag:
671	gezelter	1782
672			while(inputStream.getline(buffer, bufferSize)) {
673			line = buffer;
674
675			if(line.find("</SiteData>") != std::string::npos) {
676			break;
677			}
678
679			parseSiteLine(line);
680			}
681
682			}
683
684	tim	1024	void DumpReader::readFrameProperties(std::istream& inputStream) {
685
686			Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
687			inputStream.getline(buffer, bufferSize);
688			std::string line(buffer);
689
690			if (line.find("<FrameData>") == std::string::npos) {
691			sprintf(painCave.errMsg,
692			"DumpReader Error: Missing <FrameData>\n");
693			painCave.isFatal = 1;
694			simError();
695			}
696
697			while(inputStream.getline(buffer, bufferSize)) {
698			line = buffer;
699
700			if(line.find("</FrameData>") != std::string::npos) {
701			break;
702			}
703
704			StringTokenizer tokenizer(line, " ;\t\n\r{}:,");
705			if (!tokenizer.hasMoreTokens()) {
706	chrisfen	721	sprintf(painCave.errMsg,
707	tim	1024	"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
708	chrisfen	721	painCave.isFatal = 1;
709	tim	1024	simError();
710			}
711
712			std::string propertyName = tokenizer.nextToken();
713			if (propertyName == "Time") {
714			RealType currTime = tokenizer.nextTokenAsDouble();
715			s->setTime(currTime);
716			} else if (propertyName == "Hmat"){
717			Mat3x3d hmat;
718			hmat(0, 0) = tokenizer.nextTokenAsDouble();
719			hmat(0, 1) = tokenizer.nextTokenAsDouble();
720			hmat(0, 2) = tokenizer.nextTokenAsDouble();
721			hmat(1, 0) = tokenizer.nextTokenAsDouble();
722			hmat(1, 1) = tokenizer.nextTokenAsDouble();
723			hmat(1, 2) = tokenizer.nextTokenAsDouble();
724			hmat(2, 0) = tokenizer.nextTokenAsDouble();
725			hmat(2, 1) = tokenizer.nextTokenAsDouble();
726			hmat(2, 2) = tokenizer.nextTokenAsDouble();
727			s->setHmat(hmat);
728			} else if (propertyName == "Thermostat") {
729	gezelter	1782	pair<RealType, RealType> thermostat;
730			thermostat.first = tokenizer.nextTokenAsDouble();
731			thermostat.second = tokenizer.nextTokenAsDouble();
732			s->setThermostat(thermostat);
733	tim	1024	} else if (propertyName == "Barostat") {
734			Mat3x3d eta;
735			eta(0, 0) = tokenizer.nextTokenAsDouble();
736			eta(0, 1) = tokenizer.nextTokenAsDouble();
737			eta(0, 2) = tokenizer.nextTokenAsDouble();
738			eta(1, 0) = tokenizer.nextTokenAsDouble();
739			eta(1, 1) = tokenizer.nextTokenAsDouble();
740			eta(1, 2) = tokenizer.nextTokenAsDouble();
741			eta(2, 0) = tokenizer.nextTokenAsDouble();
742			eta(2, 1) = tokenizer.nextTokenAsDouble();
743			eta(2, 2) = tokenizer.nextTokenAsDouble();
744	gezelter	1782	s->setBarostat(eta);
745	tim	1024	} else {
746			sprintf(painCave.errMsg,
747			"DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str());
748			painCave.isFatal = 0;
749			simError();
750			}
751
752			}
753
754			}
755
756	chrisfen	721
757	gezelter	1390	}//end namespace OpenMD