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, 24107 (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 
 
#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"
 
#ifdef IS_MPI 
#include <mpi.h> 
#endif
 
 
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()); 
      
      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) { 
    int lineNo = 0; 
    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;
      bool foundOpenSiteDataTag = false;
      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::dslElectroFrame) { 
      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* integrableObject = info_->getIOIndexToIntegrableObject(index);

    if (integrableObject == 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 (integrableObject->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_) { 
              integrableObject->setPos(pos); 
            }             
            break;
        }
        case 'v' : {
            Vector3d vel;
            vel[0] = tokenizer.nextTokenAsDouble(); 
            vel[1] = tokenizer.nextTokenAsDouble(); 
            vel[2] = tokenizer.nextTokenAsDouble(); 
            if (needVel_) { 
              integrableObject->setVel(vel); 
            } 
            break;
        }

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

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

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

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

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

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

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

           Vector3d eField;
           eField[0] = tokenizer.nextTokenAsDouble(); 
           eField[1] = tokenizer.nextTokenAsDouble(); 
           eField[2] = tokenizer.nextTokenAsDouble();           
           integrableObject->setElectricField(eField);          
           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* integrableObject = info_->getIOIndexToIntegrableObject(index);
    if (integrableObject == NULL) {
      return;
    }
    StuntDouble* sd = integrableObject;

    /**
     * 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();
      RigidBody* rb = static_cast<RigidBody*>(integrableObject);
      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;
      }
      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) {

    inputStream.getline(buffer, bufferSize);
    std::string line(buffer);
    
    if (line.find("<SiteData>") == std::string::npos) {
      // site data isn't required for a simulation, so skip
      return;
    }

    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:	1764
Committed:	Tue Jul 3 18:32:27 2012 UTC (12 years, 9 months ago) by gezelter
File size:	21617 byte(s)
Log Message:	Refactored Snapshot and Stats to use the Accumulator classes. Collected a number of methods into Thermo that belonged there.
#	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			* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	gezelter	1665	* [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			#include <sys/types.h>
47			#include <sys/stat.h>
48
49			#include <iostream>
50			#include <math.h>
51
52			#include <stdio.h>
53			#include <stdlib.h>
54			#include <string.h>
55
56			#include "io/DumpReader.hpp"
57			#include "primitives/Molecule.hpp"
58			#include "utils/simError.h"
59			#include "utils/MemoryUtils.hpp"
60			#include "utils/StringTokenizer.hpp"
61	gezelter	1764	#include "brains/Thermo.hpp"
62	chrisfen	721
63			#ifdef IS_MPI
64			#include <mpi.h>
65	gezelter	1764	#endif
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	chrisfen	721	inFile_ = new std::ifstream(filename_.c_str());
79	chrisfen	996
80			if (inFile_->fail()) {
81			sprintf(painCave.errMsg,
82			"DumpReader: Cannot open file: %s\n",
83			filename_.c_str());
84			painCave.isFatal = 1;
85			simError();
86			}
87
88	chrisfen	721	#ifdef IS_MPI
89	chrisfen	996
90			}
91
92			strcpy(checkPointMsg, "Dump file opened for reading successfully.");
93	gezelter	1241	errorCheckPoint();
94	chrisfen	996
95	chrisfen	721	#endif
96	chrisfen	996
97			return;
98			}
99
100	chrisfen	721	DumpReader::~DumpReader() {
101	chrisfen	996
102	chrisfen	721	#ifdef IS_MPI
103	chrisfen	996
104	chrisfen	721	if (worldRank == 0) {
105			#endif
106	chrisfen	996
107	chrisfen	721	delete inFile_;
108	chrisfen	996
109	chrisfen	721	#ifdef IS_MPI
110	chrisfen	996
111	chrisfen	721	}
112	chrisfen	996
113	chrisfen	721	strcpy(checkPointMsg, "Dump file closed successfully.");
114	gezelter	1241	errorCheckPoint();
115	chrisfen	996
116	chrisfen	721	#endif
117	chrisfen	996
118	chrisfen	721	return;
119			}
120	chrisfen	996
121	chrisfen	721	int DumpReader::getNFrames(void) {
122
123			if (!isScanned_)
124			scanFile();
125
126			return nframes_;
127			}
128
129			void DumpReader::scanFile(void) {
130	tim	1024	int lineNo = 0;
131			std::streampos prevPos;
132	chrisfen	721	std::streampos currPos;
133	tim	1024
134	chrisfen	721	#ifdef IS_MPI
135	tim	1024
136	chrisfen	721	if (worldRank == 0) {
137			#endif // is_mpi
138	tim	1024
139			currPos = inFile_->tellg();
140			prevPos = currPos;
141			bool foundOpenSnapshotTag = false;
142			bool foundClosedSnapshotTag = false;
143	gezelter	1752	bool foundOpenSiteDataTag = false;
144	tim	1024	while(inFile_->getline(buffer, bufferSize)) {
145			++lineNo;
146
147			std::string line = buffer;
148			currPos = inFile_->tellg();
149			if (line.find("<Snapshot>")!= std::string::npos) {
150			if (foundOpenSnapshotTag) {
151	chrisfen	721	sprintf(painCave.errMsg,
152	tim	1024	"DumpReader:<Snapshot> is multiply nested at line %d in %s \n", lineNo,
153			filename_.c_str());
154	chrisfen	721	painCave.isFatal = 1;
155	tim	1024	simError();
156			}
157			foundOpenSnapshotTag = true;
158			foundClosedSnapshotTag = false;
159			framePos_.push_back(prevPos);
160
161			} else if (line.find("</Snapshot>") != std::string::npos){
162			if (!foundOpenSnapshotTag) {
163			sprintf(painCave.errMsg,
164			"DumpReader:</Snapshot> appears before <Snapshot> at line %d in %s \n", lineNo,
165			filename_.c_str());
166			painCave.isFatal = 1;
167	chrisfen	721	simError();
168	tim	1024	}
169
170			if (foundClosedSnapshotTag) {
171			sprintf(painCave.errMsg,
172			"DumpReader:</Snapshot> appears multiply nested at line %d in %s \n", lineNo,
173			filename_.c_str());
174			painCave.isFatal = 1;
175			simError();
176			}
177			foundClosedSnapshotTag = true;
178			foundOpenSnapshotTag = false;
179			}
180			prevPos = currPos;
181			}
182
183			// only found <Snapshot> for the last frame means the file is corrupted, we should discard
184			// it and give a warning message
185			if (foundOpenSnapshotTag) {
186			sprintf(painCave.errMsg,
187			"DumpReader: last frame in %s is invalid\n", filename_.c_str());
188			painCave.isFatal = 0;
189			simError();
190			framePos_.pop_back();
191			}
192
193	chrisfen	721	nframes_ = framePos_.size();
194	tim	1024
195			if (nframes_ == 0) {
196			sprintf(painCave.errMsg,
197			"DumpReader: %s does not contain a valid frame\n", filename_.c_str());
198			painCave.isFatal = 1;
199			simError();
200			}
201	chrisfen	721	#ifdef IS_MPI
202			}
203
204			MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
205	tim	1024
206	chrisfen	721	#endif // is_mpi
207	tim	1024
208	chrisfen	721	isScanned_ = true;
209			}
210
211			void DumpReader::readFrame(int whichFrame) {
212			if (!isScanned_)
213			scanFile();
214
215			int storageLayout = info_->getSnapshotManager()->getStorageLayout();
216
217			if (storageLayout & DataStorage::dslPosition) {
218			needPos_ = true;
219			} else {
220			needPos_ = false;
221			}
222
223			if (storageLayout & DataStorage::dslVelocity) {
224			needVel_ = true;
225			} else {
226			needVel_ = false;
227			}
228
229			if (storageLayout & DataStorage::dslAmat \|\| storageLayout & DataStorage::dslElectroFrame) {
230			needQuaternion_ = true;
231			} else {
232			needQuaternion_ = false;
233			}
234
235			if (storageLayout & DataStorage::dslAngularMomentum) {
236			needAngMom_ = true;
237			} else {
238			needAngMom_ = false;
239			}
240
241			readSet(whichFrame);
242	gezelter	1104
243			if (needCOMprops_) {
244			Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
245	gezelter	1764	Thermo thermo(info_);
246	gezelter	1104	Vector3d com;
247	gezelter	1764
248			if (needPos_ && needVel_) {
249			Vector3d comvel;
250			Vector3d comw;
251			thermo.getComAll(com, comvel);
252			comw = thermo.getAngularMomentum();
253			} else {
254			com = thermo.getCom();
255			}
256	gezelter	1104	}
257	chrisfen	721	}
258
259	tim	1024	void DumpReader::readSet(int whichFrame) {
260			std::string line;
261
262	chrisfen	721	#ifndef IS_MPI
263			inFile_->clear();
264			inFile_->seekg(framePos_[whichFrame]);
265	tim	1024
266			std::istream& inputStream = *inFile_;
267
268			#else
269			int masterNode = 0;
270			std::stringstream sstream;
271			if (worldRank == masterNode) {
272			std::string sendBuffer;
273
274			inFile_->clear();
275			inFile_->seekg(framePos_[whichFrame]);
276
277			while (inFile_->getline(buffer, bufferSize)) {
278
279			line = buffer;
280			sendBuffer += line;
281			sendBuffer += '\n';
282			if (line.find("</Snapshot>") != std::string::npos) {
283			break;
284			}
285			}
286
287			int sendBufferSize = sendBuffer.size();
288			MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
289			MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
290
291			sstream.str(sendBuffer);
292			} else {
293			int sendBufferSize;
294			MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
295			char * recvBuffer = new char[sendBufferSize+1];
296	gezelter	1313	assert(recvBuffer);
297			recvBuffer[sendBufferSize] = '\0';
298	tim	1024	MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
299			sstream.str(recvBuffer);
300	gezelter	1313	delete [] recvBuffer;
301	tim	1024	}
302
303			std::istream& inputStream = sstream;
304			#endif
305
306			inputStream.getline(buffer, bufferSize);
307
308			line = buffer;
309			if (line.find("<Snapshot>") == std::string::npos) {
310	chrisfen	721	sprintf(painCave.errMsg,
311	tim	1024	"DumpReader Error: can not find <Snapshot>\n");
312	chrisfen	721	painCave.isFatal = 1;
313			simError();
314			}
315	tim	1024
316			//read frameData
317			readFrameProperties(inputStream);
318
319			//read StuntDoubles
320			readStuntDoubles(inputStream);
321
322			inputStream.getline(buffer, bufferSize);
323			line = buffer;
324	gezelter	1752
325			if (line.find("<SiteData>") != std::string::npos) {
326			//read SiteData
327			readSiteData(inputStream);
328			} else {
329			if (line.find("</Snapshot>") == std::string::npos) {
330			sprintf(painCave.errMsg,
331			"DumpReader Error: can not find </Snapshot>\n");
332			painCave.isFatal = 1;
333			simError();
334			}
335			}
336	chrisfen	721	}
337
338	tim	1024	void DumpReader::parseDumpLine(const std::string& line) {
339
340
341	chrisfen	721	StringTokenizer tokenizer(line);
342			int nTokens;
343
344			nTokens = tokenizer.countTokens();
345
346	gezelter	1450	if (nTokens < 2) {
347	chrisfen	721	sprintf(painCave.errMsg,
348	tim	1024	"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
349	chrisfen	721	painCave.isFatal = 1;
350			simError();
351			}
352	tim	1024
353			int index = tokenizer.nextTokenAsInt();
354
355			StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
356
357			if (integrableObject == NULL) {
358			return;
359			}
360			std::string type = tokenizer.nextToken();
361			int size = type.size();
362	gezelter	1037
363	gezelter	1450	size_t found;
364
365			if (needPos_) {
366			found = type.find("p");
367			if (found == std::string::npos) {
368			sprintf(painCave.errMsg,
369			"DumpReader Error: StuntDouble %d has no Position\n"
370			"\tField (\"p\") specified.\n%s\n", index,
371			line.c_str());
372			painCave.isFatal = 1;
373			simError();
374			}
375			}
376
377			if (integrableObject->isDirectional()) {
378			if (needQuaternion_) {
379			found = type.find("q");
380			if (found == std::string::npos) {
381			sprintf(painCave.errMsg,
382			"DumpReader Error: Directional StuntDouble %d has no\n"
383			"\tQuaternion Field (\"q\") specified.\n%s\n", index,
384			line.c_str());
385			painCave.isFatal = 1;
386			simError();
387			}
388			}
389			}
390
391	tim	1024	for(int i = 0; i < size; ++i) {
392			switch(type[i]) {
393
394			case 'p': {
395			Vector3d pos;
396			pos[0] = tokenizer.nextTokenAsDouble();
397			pos[1] = tokenizer.nextTokenAsDouble();
398			pos[2] = tokenizer.nextTokenAsDouble();
399			if (needPos_) {
400			integrableObject->setPos(pos);
401			}
402			break;
403			}
404			case 'v' : {
405			Vector3d vel;
406			vel[0] = tokenizer.nextTokenAsDouble();
407			vel[1] = tokenizer.nextTokenAsDouble();
408			vel[2] = tokenizer.nextTokenAsDouble();
409			if (needVel_) {
410			integrableObject->setVel(vel);
411			}
412			break;
413			}
414
415			case 'q' : {
416			Quat4d q;
417			if (integrableObject->isDirectional()) {
418
419			q[0] = tokenizer.nextTokenAsDouble();
420			q[1] = tokenizer.nextTokenAsDouble();
421			q[2] = tokenizer.nextTokenAsDouble();
422			q[3] = tokenizer.nextTokenAsDouble();
423
424			RealType qlen = q.length();
425	gezelter	1390	if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
426	tim	1024
427			sprintf(painCave.errMsg,
428			"DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
429			painCave.isFatal = 1;
430			simError();
431
432			}
433
434			q.normalize();
435			if (needQuaternion_) {
436			integrableObject->setQ(q);
437			}
438			}
439			break;
440			}
441			case 'j' : {
442			Vector3d ji;
443			if (integrableObject->isDirectional()) {
444			ji[0] = tokenizer.nextTokenAsDouble();
445			ji[1] = tokenizer.nextTokenAsDouble();
446			ji[2] = tokenizer.nextTokenAsDouble();
447			if (needAngMom_) {
448			integrableObject->setJ(ji);
449			}
450			}
451	gezelter	1037	break;
452	tim	1024	}
453			case 'f': {
454
455			Vector3d force;
456			force[0] = tokenizer.nextTokenAsDouble();
457			force[1] = tokenizer.nextTokenAsDouble();
458			force[2] = tokenizer.nextTokenAsDouble();
459			integrableObject->setFrc(force);
460			break;
461			}
462			case 't' : {
463
464			Vector3d torque;
465			torque[0] = tokenizer.nextTokenAsDouble();
466			torque[1] = tokenizer.nextTokenAsDouble();
467			torque[2] = tokenizer.nextTokenAsDouble();
468			integrableObject->setTrq(torque);
469			break;
470			}
471	gezelter	1629	case 'u' : {
472
473			RealType particlePot;
474			particlePot = tokenizer.nextTokenAsDouble();
475			integrableObject->setParticlePot(particlePot);
476			break;
477			}
478	gezelter	1714	case 'c' : {
479
480			RealType flucQPos;
481			flucQPos = tokenizer.nextTokenAsDouble();
482			integrableObject->setFlucQPos(flucQPos);
483			break;
484			}
485			case 'w' : {
486
487			RealType flucQVel;
488			flucQVel = tokenizer.nextTokenAsDouble();
489			integrableObject->setFlucQVel(flucQVel);
490			break;
491			}
492			case 'g' : {
493
494			RealType flucQFrc;
495			flucQFrc = tokenizer.nextTokenAsDouble();
496			integrableObject->setFlucQFrc(flucQFrc);
497			break;
498			}
499			case 'e' : {
500
501			Vector3d eField;
502			eField[0] = tokenizer.nextTokenAsDouble();
503			eField[1] = tokenizer.nextTokenAsDouble();
504			eField[2] = tokenizer.nextTokenAsDouble();
505			integrableObject->setElectricField(eField);
506			break;
507			}
508	tim	1024	default: {
509			sprintf(painCave.errMsg,
510			"DumpReader Error: %s is an unrecognized type\n", type.c_str());
511			painCave.isFatal = 1;
512			simError();
513			break;
514			}
515
516			}
517			}
518	chrisfen	721
519	tim	1024	}
520
521
522	gezelter	1752	void DumpReader::parseSiteLine(const std::string& line) {
523
524			StringTokenizer tokenizer(line);
525			int nTokens;
526
527			nTokens = tokenizer.countTokens();
528
529			if (nTokens < 2) {
530			sprintf(painCave.errMsg,
531			"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
532			painCave.isFatal = 1;
533			simError();
534			}
535
536			/**
537			* The first token is the global integrable object index.
538			*/
539
540			int index = tokenizer.nextTokenAsInt();
541			StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
542			if (integrableObject == NULL) {
543			return;
544			}
545			StuntDouble* sd = integrableObject;
546
547			/**
548			* Test to see if the next token is an integer or not. If not,
549			* we've got data on the integrable object itself. If there is an
550			* integer, we're parsing data for a site on a rigid body.
551			*/
552
553			std::string indexTest = tokenizer.peekNextToken();
554			std::istringstream i(indexTest);
555			int siteIndex;
556			if (i >> siteIndex) {
557			// chew up this token and parse as an int:
558			siteIndex = tokenizer.nextTokenAsInt();
559			RigidBody* rb = static_cast<RigidBody*>(integrableObject);
560			sd = rb->getAtoms()[siteIndex];
561			}
562
563			/**
564			* The next token contains information on what follows.
565			*/
566			std::string type = tokenizer.nextToken();
567			int size = type.size();
568
569			for(int i = 0; i < size; ++i) {
570			switch(type[i]) {
571
572			case 'u' : {
573
574			RealType particlePot;
575			particlePot = tokenizer.nextTokenAsDouble();
576			sd->setParticlePot(particlePot);
577			break;
578			}
579			case 'c' : {
580
581			RealType flucQPos;
582			flucQPos = tokenizer.nextTokenAsDouble();
583			sd->setFlucQPos(flucQPos);
584			break;
585			}
586			case 'w' : {
587
588			RealType flucQVel;
589			flucQVel = tokenizer.nextTokenAsDouble();
590			sd->setFlucQVel(flucQVel);
591			break;
592			}
593			case 'g' : {
594
595			RealType flucQFrc;
596			flucQFrc = tokenizer.nextTokenAsDouble();
597			sd->setFlucQFrc(flucQFrc);
598			break;
599			}
600			case 'e' : {
601
602			Vector3d eField;
603			eField[0] = tokenizer.nextTokenAsDouble();
604			eField[1] = tokenizer.nextTokenAsDouble();
605			eField[2] = tokenizer.nextTokenAsDouble();
606			sd->setElectricField(eField);
607			break;
608			}
609			default: {
610			sprintf(painCave.errMsg,
611			"DumpReader Error: %s is an unrecognized type\n", type.c_str());
612			painCave.isFatal = 1;
613			simError();
614			break;
615			}
616			}
617			}
618			}
619
620
621	tim	1024	void DumpReader::readStuntDoubles(std::istream& inputStream) {
622	gezelter	1752
623	tim	1024	inputStream.getline(buffer, bufferSize);
624			std::string line(buffer);
625
626			if (line.find("<StuntDoubles>") == std::string::npos) {
627	chrisfen	721	sprintf(painCave.errMsg,
628	tim	1024	"DumpReader Error: Missing <StuntDoubles>\n");
629	chrisfen	721	painCave.isFatal = 1;
630	tim	1024	simError();
631			}
632
633			while(inputStream.getline(buffer, bufferSize)) {
634			line = buffer;
635
636			if(line.find("</StuntDoubles>") != std::string::npos) {
637			break;
638			}
639
640			parseDumpLine(line);
641			}
642
643			}
644
645	gezelter	1752	void DumpReader::readSiteData(std::istream& inputStream) {
646
647			inputStream.getline(buffer, bufferSize);
648			std::string line(buffer);
649
650			if (line.find("<SiteData>") == std::string::npos) {
651			// site data isn't required for a simulation, so skip
652			return;
653			}
654
655			while(inputStream.getline(buffer, bufferSize)) {
656			line = buffer;
657
658			if(line.find("</SiteData>") != std::string::npos) {
659			break;
660			}
661
662			parseSiteLine(line);
663			}
664
665			}
666
667	tim	1024	void DumpReader::readFrameProperties(std::istream& inputStream) {
668
669			Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
670			inputStream.getline(buffer, bufferSize);
671			std::string line(buffer);
672
673			if (line.find("<FrameData>") == std::string::npos) {
674			sprintf(painCave.errMsg,
675			"DumpReader Error: Missing <FrameData>\n");
676			painCave.isFatal = 1;
677			simError();
678			}
679
680			while(inputStream.getline(buffer, bufferSize)) {
681			line = buffer;
682
683			if(line.find("</FrameData>") != std::string::npos) {
684			break;
685			}
686
687			StringTokenizer tokenizer(line, " ;\t\n\r{}:,");
688			if (!tokenizer.hasMoreTokens()) {
689	chrisfen	721	sprintf(painCave.errMsg,
690	tim	1024	"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
691	chrisfen	721	painCave.isFatal = 1;
692	tim	1024	simError();
693			}
694
695			std::string propertyName = tokenizer.nextToken();
696			if (propertyName == "Time") {
697			RealType currTime = tokenizer.nextTokenAsDouble();
698			s->setTime(currTime);
699			} else if (propertyName == "Hmat"){
700			Mat3x3d hmat;
701			hmat(0, 0) = tokenizer.nextTokenAsDouble();
702			hmat(0, 1) = tokenizer.nextTokenAsDouble();
703			hmat(0, 2) = tokenizer.nextTokenAsDouble();
704			hmat(1, 0) = tokenizer.nextTokenAsDouble();
705			hmat(1, 1) = tokenizer.nextTokenAsDouble();
706			hmat(1, 2) = tokenizer.nextTokenAsDouble();
707			hmat(2, 0) = tokenizer.nextTokenAsDouble();
708			hmat(2, 1) = tokenizer.nextTokenAsDouble();
709			hmat(2, 2) = tokenizer.nextTokenAsDouble();
710			s->setHmat(hmat);
711			} else if (propertyName == "Thermostat") {
712	gezelter	1764	pair<RealType, RealType> thermostat;
713			thermostat.first = tokenizer.nextTokenAsDouble();
714			thermostat.second = tokenizer.nextTokenAsDouble();
715			s->setThermostat(thermostat);
716	tim	1024	} else if (propertyName == "Barostat") {
717			Mat3x3d eta;
718			eta(0, 0) = tokenizer.nextTokenAsDouble();
719			eta(0, 1) = tokenizer.nextTokenAsDouble();
720			eta(0, 2) = tokenizer.nextTokenAsDouble();
721			eta(1, 0) = tokenizer.nextTokenAsDouble();
722			eta(1, 1) = tokenizer.nextTokenAsDouble();
723			eta(1, 2) = tokenizer.nextTokenAsDouble();
724			eta(2, 0) = tokenizer.nextTokenAsDouble();
725			eta(2, 1) = tokenizer.nextTokenAsDouble();
726			eta(2, 2) = tokenizer.nextTokenAsDouble();
727	gezelter	1764	s->setBarostat(eta);
728	tim	1024	} else {
729			sprintf(painCave.errMsg,
730			"DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str());
731			painCave.isFatal = 0;
732			simError();
733			}
734
735			}
736
737			}
738
739	chrisfen	721
740	gezelter	1390	}//end namespace OpenMD