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