OOPSE/libmdtools/DumpReader.cpp

#define _FILE_OFFSET_BITS 64

#include <iostream>
#include <cmath>
#include <vector>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>

#include "ReadWrite.hpp"
#include "simError.h"

#ifdef IS_MPI
#include <mpi.h>
#include "mpiSimulation.hpp"
#define TAKE_THIS_TAG_CHAR 0
#define TAKE_THIS_TAG_INT 1
#endif // is_mpi

namespace dumpRead{

#ifdef IS_MPI
  void nodeZeroError( void );
  void anonymousNodeDie( void );
#endif // is_mpi
  
  class FilePos{
    
  public:
    FilePos(){ myPos = NULL; }
    FilePos( fpos_t* thePos ) { myPos = thePos; }
    ~FilePos(){ if( myPos != NULL ) delete myPos; }

    FilePos &operator=(fpos_t *thePos){ myPos = thePos; return *this; }

    void setPos( fpos_t *thePos ){ myPos = thePos; }
    fpos_t *getPos( void ){ return myPos; }
   
  private:
    
    fpos_t *myPos

  };

  bool operator<(FilePos a, FilePos b){
    return (a.getPos())->__pos < (b.getPos())->__pos; }

  bool operator==(FilePos a, FilePos b){
    return (a.getPos())->__pos == (b.getPos())->__pos; }

  vector<FilePos> frameStart;
}

using namespace dumpRead;


DumpReader :: DumpReader( char *in_name ){

  isScanned = false;

#ifdef IS_MPI
  if (worldRank == 0) {
#endif

  inFile = fopen(in_name, "r");
  if(inFile == NULL){
    sprintf(painCave.errMsg,
            "Cannot open file: %s\n", in_name);
    painCave.isFatal = 1;
    simError();
  }
  
  strcpy( inName, in_name);
#ifdef IS_MPI
  }
  strcpy( checkPointMsg, "Dump file opened for reading successfully." );
  MPIcheckPoint();
#endif
  return;  
}

DumpReader :: ~DumpReader( ){
#ifdef IS_MPI
  if (worldRank == 0) {
#endif
  int error;
  error = fclose( inFile );
  if( error ){
    sprintf( painCave.errMsg,
             "Error closing %s\n", inName );
    simError();
  }
#ifdef IS_MPI
  }
  strcpy( checkPointMsg, "Dump file closed successfully." );
  MPIcheckPoint();
#endif

  return;
}

int DumpReader::getNframes( void ){

  if( !isScanned ) scanFile();
  return nFrames;
}

void DumpReader::scanFile( void ){

  int vectorSize;
  int i, j, k;
  int lineNum = 0;
  char readBuffer[2000];
  char* foo;
  fpos_t *currPos;

#ifdef IS_MPI
  if( worldRank == 0 ){
#endif // is_mpi
    
    rewind( inFile );
    
    currPos = new fpos_t;
    fgetpos( inFile, currPos );
    fgets( readBuffer, sizeof( readBuffer ), in_file );
    lineNum++;
    if( feof( in_file ) ){
      sprintf( painCave.errMsg,
               "File \"%s\" ended unexpectedly at line %d\n", 
               inName,
               lineNum );
      painCave.isFatal = 1;
      simError();
    }
    
    while( !feof( in_file ) ){
      
      frameStart.push_back( FilePos( currPos ) );
      i = atoi(readBuffer);
      
      fgets( readBuffer, sizeof( readBuffer ), in_file );
      lineNum++;    
      if( feof( in_file ) ){
        sprintf( painCave.errMsg,
                 "File \"%s\" ended unexpectedly at line %d\n",
                 inName,
                 lineNum );
        painCave.isFatal = 1;
        simError();
      }
      
      if(outTime){
        foo = strtok( readBuffer, " ,;\t" );
        time = atof( foo );
      }
      
      for(j=0; j<i; j++){
        
        fgets( readBuffer, sizeof( readBuffer ), in_file );
        lineNum++;    
        if( feof( in_file ) ){
          sprintf( painCave.errMsg,
                   "File \"%s\" ended unexpectedly at line %d,"
                   " with atom %d\n", 
                   inName,
                   lineNum, 
                   j );
          painCave.isFatal = 1;
          simError();
        }
        
      }
      
      currPos = new fpos_t;
      fgetpos( inFile, currPos );
      fgets( readBuffer, sizeof( readBuffer ), in_file );
      lineNum++;
    }
    
    delete currPos;
    vectorSize = frameStart.size();
    rewind( inFile );
    
#ifdef IS_MPI
  }
  strcpy( checkPointMsg, "Successfully scanned DumpFile\n" );
  MPIcheckPoint();
#endif // is_mpi

  nFrames = vectorSize;
  isScanned = true;
}

void DumpReader :: readFrame( SimInfo* the_simnfo, int whichFrame){

  simnfo = the_simnfo;

  this->readSet( whichFrame );
}


void DumpReader :: readSet( int whichFrame ){

  int i, j, done, which_node, which_atom; // loop counter

  const int BUFFERSIZE = 2000; // size of the read buffer
  int n_atoms; // the number of atoms
  char read_buffer[BUFFERSIZE]; //the line buffer for reading 
#ifdef IS_MPI
  char send_buffer[BUFFERSIZE];
#endif

  char *eof_test; // ptr to see when we reach the end of the file 
  char *parseErr;
  int procIndex;
  double boxMat[9];
  double theBoxMat3[3][3];
  
  fpos_t *framePos;
  
#ifndef IS_MPI
  
  framePos = startFrame[whichFrame].getPos();

  
  eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
  if( eof_test == NULL ){
    sprintf( painCave.errMsg,
             "DumpReader error: error reading 1st line of \"%s\"\n",
             inName );
    painCave.isFatal = 1;
    simError();
  }
  
  n_atoms = atoi( read_buffer );

  Atom **atoms = simnfo->atoms;
  DirectionalAtom* dAtom;

  if( n_atoms != simnfo->n_atoms ){
    sprintf( painCave.errMsg,
             "DumpReader error. %s n_atoms, %d, "
             "does not match the BASS file's n_atoms, %d.\n",
             inName, n_atoms, simnfo->n_atoms );
    painCave.isFatal = 1;
    simError();
  }
  
  //read the box mat from the comment line 
  
  eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
  if(eof_test == NULL){
    sprintf( painCave.errMsg,
             "error in reading commment in %s\n", inName);
    painCave.isFatal = 1;
    simError();
  }

  parseErr = parseCommentLine( read_buffer, time, boxMat );
  if( parseErr != NULL ){
    strcpy( painCave.errMsg, parseErr );
    painCave.isFatal = 1;
    simError();
  }

  simnfo->setTime( time );
  
  for(i=0;i<3;i++)
    for(j=0;j<3;j++) theBoxMat3[i][j] = boxMat[3*j+i];
  

  simnfo->setBoxM( theBoxMat3 );
  

  for( i=0; i < n_atoms; i++){
    
    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
    if(eof_test == NULL){
      sprintf(painCave.errMsg,
              "error in reading file %s\n"
              "natoms  = %d; index = %d\n"
              "error reading the line from the file.\n",
              inName, n_atoms, i );
      painCave.isFatal = 1;
      simError();
    }

    
    parseErr = parseDumpLine( read_buffer, i );
    if( parseErr != NULL ){
      strcpy( painCave.errMsg, parseErr );
      painCave.isFatal = 1;
      simError();
    }
  }


  // MPI Section of code..........
#else //IS_MPI

  // first thing first, suspend fatalities.
  painCave.isEventLoop = 1;

  int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
  int haveError;
  
  MPI_Status istatus;
  int *AtomToProcMap = mpiSim->getAtomToProcMap();

  
  haveError = 0;
  if (worldRank == 0) {

    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
    if( eof_test == NULL ){
      sprintf( painCave.errMsg,
               "Error reading 1st line of %d \n ",inName);
      haveError = 1;
      simError();
    }
    
    n_atoms = atoi( read_buffer );
    
    Atom **atoms = simnfo->atoms;
    DirectionalAtom* dAtom;
 
    // Check to see that the number of atoms in the intial configuration file is the
    // same as declared in simBass.
    
    if( n_atoms != mpiSim->getTotAtoms() ){
      sprintf( painCave.errMsg,
               "Initialize from File error. %s n_atoms, %d, "
               "does not match the BASS file's n_atoms, %d.\n",
               inName, n_atoms, simnfo->n_atoms );
      haveError= 1;
      simError();
    }
    
    //read the time and boxMat from the comment line 
    
    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
    if(eof_test == NULL){
      sprintf( painCave.errMsg,
               "error in reading commment in %s\n", inName);
      haveError = 1;
      simError();
    }
    
    parseErr = parseCommentLine( read_buffer, time, boxMat );
    if( parseErr != NULL ){
      strcpy( painCave.errMsg, parseErr );
      haveError = 1;
      simError();
    }

    MPI_Bcast(time, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD );
    
    MPI_Bcast(boxMat, 9, MPI_DOUBLE, 0, MPI_COMM_WORLD );
    
    if(haveError) nodeZeroError();
    
    for (i=0 ; i < mpiSim->getTotAtoms(); i++) {
      
      eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
      if(eof_test == NULL){
        sprintf(painCave.errMsg,
                "error in reading file %s\n"
                "natoms  = %d; index = %d\n"
                "error reading the line from the file.\n",
                inName, n_atoms, i );
        haveError= 1;
        simError();
      }
    
      if(haveError) nodeZeroError();

      // Get the Node number which wants this atom:
      which_node = AtomToProcMap[i];     
      if (which_node == 0) {
        parseErr = parseDumpLine( read_buffer, i );
        if( parseErr != NULL ){
          strcpy( painCave.errMsg, parseErr );
          haveError = 1;
          simError();
        }    
        if(haveError) nodeZeroError();
      } 
      
      else {
       
        myStatus = 1;
        MPI_Send(&myStatus, 1, MPI_INT, which_node, 
                 TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
        MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node, 
                 TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
        MPI_Send(&i, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
                 MPI_COMM_WORLD);
        MPI_Recv(&myStatus, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT, 
                 MPI_COMM_WORLD, &istatus);
        
        if(!myStatus) nodeZeroError();
      }
    }
    myStatus = -1;
    for (j = 0; j < mpiSim->getNumberProcessors(); j++) {      
      MPI_Send( &myStatus, 1, MPI_INT, j, 
                TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
    }
    
  } else {
    
    MPI_Bcast(time, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(boxMat, 9, MPI_DOUBLE, 0, MPI_COMM_WORLD);

    done = 0;
    while (!done) {

      MPI_Recv(&myStatus, 1, MPI_INT, 0, 
               TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
      
      if(!myStatus) anonymousNodeDie();
      
      if(myStatus < 0) break;

      MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, 0, 
               TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
      MPI_Recv(&which_atom, 1, MPI_INT, 0, 
               TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
      
      myStatus = 1;
      parseErr = parseDumpLine( read_buffer, which_atom );
      if( parseErr != NULL ){
        strcpy( painCave.errMsg, parseErr );
        myStatus = 0;;
        simError();
      }
      
      MPI_Send( &myStatus, 1, MPI_INT, 0, 
                TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
      
    }
  }
     
 // last  thing last, enable  fatalities.
  painCave.isEventLoop = 0;

  simnfo->setTime( time );

  for(i=0;i<3;i++)
    for(j=0;j<3;j++) theBoxMat3[i][j] = boxMat[3*j+i];

  simnfo->setBoxM( theBoxMat3 );

   
#endif
}

char* DumpReader::parseDumpLine(char* readLine, int globalIndex){

  char *foo; // the pointer to the current string token 
  
  double rx, ry, rz; // position place holders
  double vx, vy, vz; // velocity placeholders
  double q[4]; // the quaternions
  double jx, jy, jz; // angular velocity placeholders;
  double qSqr, qLength; // needed to normalize the quaternion vector.
  
  Atom **atoms = simnfo->atoms;
  DirectionalAtom* dAtom;
  
  int j, n_atoms, atomIndex;

#ifdef IS_MPI
  n_atoms = mpiSim->getTotAtoms();
  atomIndex=-1;        
  for (j=0; j < mpiSim->getMyNlocal(); j++) {
    if (atoms[j]->getGlobalIndex() == globalIndex) atomIndex = j;
  }
  if (atomIndex == -1) {
    sprintf( painCave.errMsg,
             "Initialize from file error. Atom at index %d "
             "in file %s does not exist on processor %d .\n",
             globalIndex, inName, mpiSim->getMyNode() );
    return strdup( painCave.errMsg );
  }  
#else
  n_atoms = simnfo->n_atoms;
  atomIndex = globalIndex;
#endif // is_mpi

  // set the string tokenizer
  
  foo = strtok(readLine, " ,;\t");
  
  // check the atom name to the current atom
  
  if( strcmp( foo, atoms[atomIndex]->getType() ) ){
    sprintf( painCave.errMsg,
             "Initialize from file error. Atom %s at index %d "
             "in file %s does not"
             " match the BASS atom %s.\n",
             foo, atomIndex, inName, atoms[atomIndex]->getType() );
    return strdup( painCave.errMsg );
  }
    
  // get the positions

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading postition x from %s\n"
             "natoms  = %d, index = %d\n",
             inName, n_atoms, atomIndex );
    return strdup( painCave.errMsg );
  }
  rx = atof( foo );
   
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading postition y from %s\n"
             "natoms  = %d, index = %d\n",
             inName, n_atoms, atomIndex );
    return strdup( painCave.errMsg );
  }
  ry = atof( foo );
     
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading postition z from %s\n"
             "natoms  = %d, index = %d\n",
             inName, n_atoms, atomIndex );
    return strdup( painCave.errMsg );
  }
  rz = atof( foo );    


  // get the velocities

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading velocity x from %s\n"
             "natoms  = %d, index = %d\n",
             inName, n_atoms, atomIndex );
    return strdup( painCave.errMsg );
  }
  vx = atof( foo );
    
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading velocity y from %s\n"
             "natoms  = %d, index = %d\n",
             inName, n_atoms, atomIndex );
    return strdup( painCave.errMsg );
  }
  vy = atof( foo );
    
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading velocity z from %s\n"
             "natoms  = %d, index = %d\n",
             inName, n_atoms, atomIndex );
    return strdup( painCave.errMsg );
  }
  vz = atof( foo );
    
    
  // get the quaternions
    
  if( atoms[atomIndex]->isDirectional() ){
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf(painCave.errMsg,
              "error in reading quaternion 0 from %s\n"
              "natoms  = %d, index = %d\n",
              inName, n_atoms, atomIndex );
      return strdup( painCave.errMsg );
    }
    q[0] = atof( foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
               "error in reading quaternion 1 from %s\n"
               "natoms  = %d, index = %d\n",
               inName, n_atoms, atomIndex );
      return strdup( painCave.errMsg );
    }
    q[1] = atof( foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
               "error in reading quaternion 2 from %s\n"
               "natoms  = %d, index = %d\n",
               inName, n_atoms, atomIndex );
      return strdup( painCave.errMsg );
    }
    q[2] = atof( foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
               "error in reading quaternion 3 from %s\n"
               "natoms  = %d, index = %d\n",
               inName, n_atoms, atomIndex );
      return strdup( painCave.errMsg );
    }
    q[3] = atof( foo );
      
    // get the angular velocities
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
               "error in reading angular momentum jx from %s\n"
               "natoms  = %d, index = %d\n",
               inName, n_atoms, atomIndex );
      return strdup( painCave.errMsg );
    }
    jx = atof( foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
               "error in reading angular momentum jy from %s\n"
               "natoms  = %d, index = %d\n",
               inName, n_atoms, atomIndex );
      return strdup( painCave.errMsg );
    }
    jy = atof(foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
               "error in reading angular momentum jz from %s\n"
               "natoms  = %d, index = %d\n",
               inName, n_atoms, atomIndex );
      return strdup( painCave.errMsg );
    }
    jz = atof( foo );
      
    dAtom = ( DirectionalAtom* )atoms[atomIndex];

    // check that the quaternion vector is normalized

    qSqr = (q[0] * q[0]) + (q[1] * q[1]) + (q[2] * q[2]) + (q[3] * q[3]);
      
    qLength = sqrt( qSqr );
    q[0] = q[0] / qLength;
    q[1] = q[1] / qLength;
    q[2] = q[2] / qLength;
    q[3] = q[3] / qLength;
      
    dAtom->setQ( q );
      
    // add the angular velocities

    dAtom->setJx( jx );
    dAtom->setJy( jy );
    dAtom->setJz( jz );
  }
    
  // add the positions and velocities to the atom
    
  atoms[atomIndex]->setX( rx );
  atoms[atomIndex]->setY( ry );
  atoms[atomIndex]->setZ( rz );
    
  atoms[atomIndex]->set_vx( vx );
  atoms[atomIndex]->set_vy( vy );
  atoms[atomIndex]->set_vz( vz );

  return NULL;
}


char* DumpReader::parseCommentLine(char* readLine, double time, 
                                   double boxMat[9]){

  char *foo; // the pointer to the current string token 
  int j;

  // set the string tokenizer
  
  foo = strtok(readLine, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading time from %s\n",
             inName );
    return strdup( painCave.errMsg );
  }
  time = atof( foo );

  // get the Hx vector

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading Hx[0] from %s\n",
             inName );
    return strdup( painCave.errMsg );
  }
  boxMat[0] = atof( foo );
   
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading Hx[1] from %s\n",
             inName );
    return strdup( painCave.errMsg );
  }
  boxMat[1] = atof( foo );
     
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading Hx[2] from %s\n",
             inName );
    return strdup( painCave.errMsg );
  }
  boxMat[2] = atof( foo );    

  // get the Hy vector

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading Hy[0] from %s\n",
             inName );
    return strdup( painCave.errMsg );
  }
  boxMat[3] = atof( foo );
   
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading Hy[1] from %s\n",
             inName );
    return strdup( painCave.errMsg );
  }
  boxMat[4] = atof( foo );
     
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading Hy[2] from %s\n",
             inName );
    return strdup( painCave.errMsg );
  }
  boxMat[5] = atof( foo );    

  // get the Hz vector

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading Hz[0] from %s\n",
             inName );
    return strdup( painCave.errMsg );
  }
  boxMat[6] = atof( foo );
   
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading Hz[1] from %s\n",
             inName );
    return strdup( painCave.errMsg );
  }
  boxMat[7] = atof( foo );
     
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading Hz[2] from %s\n",
             inName );
    return strdup( painCave.errMsg );
  }
  boxMat[8] = atof( foo );    

  return NULL;
}


#ifdef IS_MPI

// a couple of functions to let us escape the read loop

void initFile::nodeZeroError( void ){
  int j, myStatus;
  
  myStatus = 0;
  for (j = 0; j < mpiSim->getNumberProcessors(); j++) {      
    MPI_Send( &myStatus, 1, MPI_INT, j, 
              TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
  }  
  

  MPI_Finalize();
  exit (0);
  
}

void initFile::anonymousNodeDie( void ){

  MPI_Finalize();
  exit (0);
}

#endif //is_mpi
Revision:	656
Committed:	Tue Jul 29 16:32:37 2003 UTC (22 years, 3 months ago) by mmeineke
File size:	19190 byte(s)
Log Message:	working on the props code
#	Content
1	#define _FILE_OFFSET_BITS 64
2
3	#include <iostream>
4	#include <cmath>
5	#include <vector>
6
7	#include <stdio.h>
8	#include <stdlib.h>
9	#include <string.h>
10	#include <unistd.h>
11	#include <sys/types.h>
12	#include <sys/stat.h>
13
14	#include "ReadWrite.hpp"
15	#include "simError.h"
16
17	#ifdef IS_MPI
18	#include <mpi.h>
19	#include "mpiSimulation.hpp"
20	#define TAKE_THIS_TAG_CHAR 0
21	#define TAKE_THIS_TAG_INT 1
22	#endif // is_mpi
23
24	namespace dumpRead{
25
26	#ifdef IS_MPI
27	void nodeZeroError( void );
28	void anonymousNodeDie( void );
29	#endif // is_mpi
30
31	class FilePos{
32
33	public:
34	FilePos(){ myPos = NULL; }
35	FilePos( fpos_t* thePos ) { myPos = thePos; }
36	~FilePos(){ if( myPos != NULL ) delete myPos; }
37
38	FilePos &operator=(fpos_t thePos){ myPos = thePos; return this; }
39
40	void setPos( fpos_t *thePos ){ myPos = thePos; }
41	fpos_t *getPos( void ){ return myPos; }
42
43	private:
44
45	fpos_t *myPos
46
47	};
48
49	bool operator<(FilePos a, FilePos b){
50	return (a.getPos())->__pos < (b.getPos())->__pos; }
51
52	bool operator==(FilePos a, FilePos b){
53	return (a.getPos())->__pos == (b.getPos())->__pos; }
54
55	vector<FilePos> frameStart;
56	}
57
58	using namespace dumpRead;
59
60
61	DumpReader :: DumpReader( char *in_name ){
62
63	isScanned = false;
64
65	#ifdef IS_MPI
66	if (worldRank == 0) {
67	#endif
68
69	inFile = fopen(in_name, "r");
70	if(inFile == NULL){
71	sprintf(painCave.errMsg,
72	"Cannot open file: %s\n", in_name);
73	painCave.isFatal = 1;
74	simError();
75	}
76
77	strcpy( inName, in_name);
78	#ifdef IS_MPI
79	}
80	strcpy( checkPointMsg, "Dump file opened for reading successfully." );
81	MPIcheckPoint();
82	#endif
83	return;
84	}
85
86	DumpReader :: ~DumpReader( ){
87	#ifdef IS_MPI
88	if (worldRank == 0) {
89	#endif
90	int error;
91	error = fclose( inFile );
92	if( error ){
93	sprintf( painCave.errMsg,
94	"Error closing %s\n", inName );
95	simError();
96	}
97	#ifdef IS_MPI
98	}
99	strcpy( checkPointMsg, "Dump file closed successfully." );
100	MPIcheckPoint();
101	#endif
102
103	return;
104	}
105
106	int DumpReader::getNframes( void ){
107
108	if( !isScanned ) scanFile();
109	return nFrames;
110	}
111
112	void DumpReader::scanFile( void ){
113
114	int vectorSize;
115	int i, j, k;
116	int lineNum = 0;
117	char readBuffer[2000];
118	char* foo;
119	fpos_t *currPos;
120
121	#ifdef IS_MPI
122	if( worldRank == 0 ){
123	#endif // is_mpi
124
125	rewind( inFile );
126
127	currPos = new fpos_t;
128	fgetpos( inFile, currPos );
129	fgets( readBuffer, sizeof( readBuffer ), in_file );
130	lineNum++;
131	if( feof( in_file ) ){
132	sprintf( painCave.errMsg,
133	"File \"%s\" ended unexpectedly at line %d\n",
134	inName,
135	lineNum );
136	painCave.isFatal = 1;
137	simError();
138	}
139
140	while( !feof( in_file ) ){
141
142	frameStart.push_back( FilePos( currPos ) );
143	i = atoi(readBuffer);
144
145	fgets( readBuffer, sizeof( readBuffer ), in_file );
146	lineNum++;
147	if( feof( in_file ) ){
148	sprintf( painCave.errMsg,
149	"File \"%s\" ended unexpectedly at line %d\n",
150	inName,
151	lineNum );
152	painCave.isFatal = 1;
153	simError();
154	}
155
156	if(outTime){
157	foo = strtok( readBuffer, " ,;\t" );
158	time = atof( foo );
159	}
160
161	for(j=0; j<i; j++){
162
163	fgets( readBuffer, sizeof( readBuffer ), in_file );
164	lineNum++;
165	if( feof( in_file ) ){
166	sprintf( painCave.errMsg,
167	"File \"%s\" ended unexpectedly at line %d,"
168	" with atom %d\n",
169	inName,
170	lineNum,
171	j );
172	painCave.isFatal = 1;
173	simError();
174	}
175
176	}
177
178	currPos = new fpos_t;
179	fgetpos( inFile, currPos );
180	fgets( readBuffer, sizeof( readBuffer ), in_file );
181	lineNum++;
182	}
183
184	delete currPos;
185	vectorSize = frameStart.size();
186	rewind( inFile );
187
188	#ifdef IS_MPI
189	}
190	strcpy( checkPointMsg, "Successfully scanned DumpFile\n" );
191	MPIcheckPoint();
192	#endif // is_mpi
193
194	nFrames = vectorSize;
195	isScanned = true;
196	}
197
198	void DumpReader :: readFrame( SimInfo* the_simnfo, int whichFrame){
199
200	simnfo = the_simnfo;
201
202	this->readSet( whichFrame );
203	}
204
205
206
207	void DumpReader :: readSet( int whichFrame ){
208
209	int i, j, done, which_node, which_atom; // loop counter
210
211	const int BUFFERSIZE = 2000; // size of the read buffer
212	int n_atoms; // the number of atoms
213	char read_buffer[BUFFERSIZE]; //the line buffer for reading
214	#ifdef IS_MPI
215	char send_buffer[BUFFERSIZE];
216	#endif
217
218	char *eof_test; // ptr to see when we reach the end of the file
219	char *parseErr;
220	int procIndex;
221	double boxMat[9];
222	double theBoxMat3[3][3];
223
224	fpos_t *framePos;
225
226	#ifndef IS_MPI
227
228	framePos = startFrame[whichFrame].getPos();
229
230
231
232
233	eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
234	if( eof_test == NULL ){
235	sprintf( painCave.errMsg,
236	"DumpReader error: error reading 1st line of \"%s\"\n",
237	inName );
238	painCave.isFatal = 1;
239	simError();
240	}
241
242	n_atoms = atoi( read_buffer );
243
244	Atom **atoms = simnfo->atoms;
245	DirectionalAtom* dAtom;
246
247	if( n_atoms != simnfo->n_atoms ){
248	sprintf( painCave.errMsg,
249	"DumpReader error. %s n_atoms, %d, "
250	"does not match the BASS file's n_atoms, %d.\n",
251	inName, n_atoms, simnfo->n_atoms );
252	painCave.isFatal = 1;
253	simError();
254	}
255
256	//read the box mat from the comment line
257
258	eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
259	if(eof_test == NULL){
260	sprintf( painCave.errMsg,
261	"error in reading commment in %s\n", inName);
262	painCave.isFatal = 1;
263	simError();
264	}
265
266	parseErr = parseCommentLine( read_buffer, time, boxMat );
267	if( parseErr != NULL ){
268	strcpy( painCave.errMsg, parseErr );
269	painCave.isFatal = 1;
270	simError();
271	}
272
273	simnfo->setTime( time );
274
275	for(i=0;i<3;i++)
276	for(j=0;j<3;j++) theBoxMat3[i][j] = boxMat[3*j+i];
277
278
279	simnfo->setBoxM( theBoxMat3 );
280
281
282	for( i=0; i < n_atoms; i++){
283
284	eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
285	if(eof_test == NULL){
286	sprintf(painCave.errMsg,
287	"error in reading file %s\n"
288	"natoms = %d; index = %d\n"
289	"error reading the line from the file.\n",
290	inName, n_atoms, i );
291	painCave.isFatal = 1;
292	simError();
293	}
294
295
296	parseErr = parseDumpLine( read_buffer, i );
297	if( parseErr != NULL ){
298	strcpy( painCave.errMsg, parseErr );
299	painCave.isFatal = 1;
300	simError();
301	}
302	}
303
304
305	// MPI Section of code..........
306	#else //IS_MPI
307
308	// first thing first, suspend fatalities.
309	painCave.isEventLoop = 1;
310
311	int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
312	int haveError;
313
314	MPI_Status istatus;
315	int *AtomToProcMap = mpiSim->getAtomToProcMap();
316
317
318	haveError = 0;
319	if (worldRank == 0) {
320
321	eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
322	if( eof_test == NULL ){
323	sprintf( painCave.errMsg,
324	"Error reading 1st line of %d \n ",inName);
325	haveError = 1;
326	simError();
327	}
328
329	n_atoms = atoi( read_buffer );
330
331	Atom **atoms = simnfo->atoms;
332	DirectionalAtom* dAtom;
333
334	// Check to see that the number of atoms in the intial configuration file is the
335	// same as declared in simBass.
336
337	if( n_atoms != mpiSim->getTotAtoms() ){
338	sprintf( painCave.errMsg,
339	"Initialize from File error. %s n_atoms, %d, "
340	"does not match the BASS file's n_atoms, %d.\n",
341	inName, n_atoms, simnfo->n_atoms );
342	haveError= 1;
343	simError();
344	}
345
346	//read the time and boxMat from the comment line
347
348	eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
349	if(eof_test == NULL){
350	sprintf( painCave.errMsg,
351	"error in reading commment in %s\n", inName);
352	haveError = 1;
353	simError();
354	}
355
356	parseErr = parseCommentLine( read_buffer, time, boxMat );
357	if( parseErr != NULL ){
358	strcpy( painCave.errMsg, parseErr );
359	haveError = 1;
360	simError();
361	}
362
363	MPI_Bcast(time, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD );
364
365	MPI_Bcast(boxMat, 9, MPI_DOUBLE, 0, MPI_COMM_WORLD );
366
367	if(haveError) nodeZeroError();
368
369	for (i=0 ; i < mpiSim->getTotAtoms(); i++) {
370
371	eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
372	if(eof_test == NULL){
373	sprintf(painCave.errMsg,
374	"error in reading file %s\n"
375	"natoms = %d; index = %d\n"
376	"error reading the line from the file.\n",
377	inName, n_atoms, i );
378	haveError= 1;
379	simError();
380	}
381
382	if(haveError) nodeZeroError();
383
384	// Get the Node number which wants this atom:
385	which_node = AtomToProcMap[i];
386	if (which_node == 0) {
387	parseErr = parseDumpLine( read_buffer, i );
388	if( parseErr != NULL ){
389	strcpy( painCave.errMsg, parseErr );
390	haveError = 1;
391	simError();
392	}
393	if(haveError) nodeZeroError();
394	}
395
396	else {
397
398	myStatus = 1;
399	MPI_Send(&myStatus, 1, MPI_INT, which_node,
400	TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
401	MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
402	TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
403	MPI_Send(&i, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
404	MPI_COMM_WORLD);
405	MPI_Recv(&myStatus, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
406	MPI_COMM_WORLD, &istatus);
407
408	if(!myStatus) nodeZeroError();
409	}
410	}
411	myStatus = -1;
412	for (j = 0; j < mpiSim->getNumberProcessors(); j++) {
413	MPI_Send( &myStatus, 1, MPI_INT, j,
414	TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
415	}
416
417	} else {
418
419	MPI_Bcast(time, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
420	MPI_Bcast(boxMat, 9, MPI_DOUBLE, 0, MPI_COMM_WORLD);
421
422	done = 0;
423	while (!done) {
424
425	MPI_Recv(&myStatus, 1, MPI_INT, 0,
426	TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
427
428	if(!myStatus) anonymousNodeDie();
429
430	if(myStatus < 0) break;
431
432	MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, 0,
433	TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
434	MPI_Recv(&which_atom, 1, MPI_INT, 0,
435	TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
436
437	myStatus = 1;
438	parseErr = parseDumpLine( read_buffer, which_atom );
439	if( parseErr != NULL ){
440	strcpy( painCave.errMsg, parseErr );
441	myStatus = 0;;
442	simError();
443	}
444
445	MPI_Send( &myStatus, 1, MPI_INT, 0,
446	TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
447
448	}
449	}
450
451	// last thing last, enable fatalities.
452	painCave.isEventLoop = 0;
453
454	simnfo->setTime( time );
455
456	for(i=0;i<3;i++)
457	for(j=0;j<3;j++) theBoxMat3[i][j] = boxMat[3*j+i];
458
459	simnfo->setBoxM( theBoxMat3 );
460
461
462	#endif
463	}
464
465	char* DumpReader::parseDumpLine(char* readLine, int globalIndex){
466
467	char *foo; // the pointer to the current string token
468
469	double rx, ry, rz; // position place holders
470	double vx, vy, vz; // velocity placeholders
471	double q[4]; // the quaternions
472	double jx, jy, jz; // angular velocity placeholders;
473	double qSqr, qLength; // needed to normalize the quaternion vector.
474
475	Atom **atoms = simnfo->atoms;
476	DirectionalAtom* dAtom;
477
478	int j, n_atoms, atomIndex;
479
480	#ifdef IS_MPI
481	n_atoms = mpiSim->getTotAtoms();
482	atomIndex=-1;
483	for (j=0; j < mpiSim->getMyNlocal(); j++) {
484	if (atoms[j]->getGlobalIndex() == globalIndex) atomIndex = j;
485	}
486	if (atomIndex == -1) {
487	sprintf( painCave.errMsg,
488	"Initialize from file error. Atom at index %d "
489	"in file %s does not exist on processor %d .\n",
490	globalIndex, inName, mpiSim->getMyNode() );
491	return strdup( painCave.errMsg );
492	}
493	#else
494	n_atoms = simnfo->n_atoms;
495	atomIndex = globalIndex;
496	#endif // is_mpi
497
498	// set the string tokenizer
499
500	foo = strtok(readLine, " ,;\t");
501
502	// check the atom name to the current atom
503
504	if( strcmp( foo, atoms[atomIndex]->getType() ) ){
505	sprintf( painCave.errMsg,
506	"Initialize from file error. Atom %s at index %d "
507	"in file %s does not"
508	" match the BASS atom %s.\n",
509	foo, atomIndex, inName, atoms[atomIndex]->getType() );
510	return strdup( painCave.errMsg );
511	}
512
513	// get the positions
514
515	foo = strtok(NULL, " ,;\t");
516	if(foo == NULL){
517	sprintf( painCave.errMsg,
518	"error in reading postition x from %s\n"
519	"natoms = %d, index = %d\n",
520	inName, n_atoms, atomIndex );
521	return strdup( painCave.errMsg );
522	}
523	rx = atof( foo );
524
525	foo = strtok(NULL, " ,;\t");
526	if(foo == NULL){
527	sprintf( painCave.errMsg,
528	"error in reading postition y from %s\n"
529	"natoms = %d, index = %d\n",
530	inName, n_atoms, atomIndex );
531	return strdup( painCave.errMsg );
532	}
533	ry = atof( foo );
534
535	foo = strtok(NULL, " ,;\t");
536	if(foo == NULL){
537	sprintf( painCave.errMsg,
538	"error in reading postition z from %s\n"
539	"natoms = %d, index = %d\n",
540	inName, n_atoms, atomIndex );
541	return strdup( painCave.errMsg );
542	}
543	rz = atof( foo );
544
545
546	// get the velocities
547
548	foo = strtok(NULL, " ,;\t");
549	if(foo == NULL){
550	sprintf( painCave.errMsg,
551	"error in reading velocity x from %s\n"
552	"natoms = %d, index = %d\n",
553	inName, n_atoms, atomIndex );
554	return strdup( painCave.errMsg );
555	}
556	vx = atof( foo );
557
558	foo = strtok(NULL, " ,;\t");
559	if(foo == NULL){
560	sprintf( painCave.errMsg,
561	"error in reading velocity y from %s\n"
562	"natoms = %d, index = %d\n",
563	inName, n_atoms, atomIndex );
564	return strdup( painCave.errMsg );
565	}
566	vy = atof( foo );
567
568	foo = strtok(NULL, " ,;\t");
569	if(foo == NULL){
570	sprintf( painCave.errMsg,
571	"error in reading velocity z from %s\n"
572	"natoms = %d, index = %d\n",
573	inName, n_atoms, atomIndex );
574	return strdup( painCave.errMsg );
575	}
576	vz = atof( foo );
577
578
579	// get the quaternions
580
581	if( atoms[atomIndex]->isDirectional() ){
582
583	foo = strtok(NULL, " ,;\t");
584	if(foo == NULL){
585	sprintf(painCave.errMsg,
586	"error in reading quaternion 0 from %s\n"
587	"natoms = %d, index = %d\n",
588	inName, n_atoms, atomIndex );
589	return strdup( painCave.errMsg );
590	}
591	q[0] = atof( foo );
592
593	foo = strtok(NULL, " ,;\t");
594	if(foo == NULL){
595	sprintf( painCave.errMsg,
596	"error in reading quaternion 1 from %s\n"
597	"natoms = %d, index = %d\n",
598	inName, n_atoms, atomIndex );
599	return strdup( painCave.errMsg );
600	}
601	q[1] = atof( foo );
602
603	foo = strtok(NULL, " ,;\t");
604	if(foo == NULL){
605	sprintf( painCave.errMsg,
606	"error in reading quaternion 2 from %s\n"
607	"natoms = %d, index = %d\n",
608	inName, n_atoms, atomIndex );
609	return strdup( painCave.errMsg );
610	}
611	q[2] = atof( foo );
612
613	foo = strtok(NULL, " ,;\t");
614	if(foo == NULL){
615	sprintf( painCave.errMsg,
616	"error in reading quaternion 3 from %s\n"
617	"natoms = %d, index = %d\n",
618	inName, n_atoms, atomIndex );
619	return strdup( painCave.errMsg );
620	}
621	q[3] = atof( foo );
622
623	// get the angular velocities
624
625	foo = strtok(NULL, " ,;\t");
626	if(foo == NULL){
627	sprintf( painCave.errMsg,
628	"error in reading angular momentum jx from %s\n"
629	"natoms = %d, index = %d\n",
630	inName, n_atoms, atomIndex );
631	return strdup( painCave.errMsg );
632	}
633	jx = atof( foo );
634
635	foo = strtok(NULL, " ,;\t");
636	if(foo == NULL){
637	sprintf( painCave.errMsg,
638	"error in reading angular momentum jy from %s\n"
639	"natoms = %d, index = %d\n",
640	inName, n_atoms, atomIndex );
641	return strdup( painCave.errMsg );
642	}
643	jy = atof(foo );
644
645	foo = strtok(NULL, " ,;\t");
646	if(foo == NULL){
647	sprintf( painCave.errMsg,
648	"error in reading angular momentum jz from %s\n"
649	"natoms = %d, index = %d\n",
650	inName, n_atoms, atomIndex );
651	return strdup( painCave.errMsg );
652	}
653	jz = atof( foo );
654
655	dAtom = ( DirectionalAtom* )atoms[atomIndex];
656
657	// check that the quaternion vector is normalized
658
659	qSqr = (q[0] * q[0]) + (q[1] * q[1]) + (q[2] * q[2]) + (q[3] * q[3]);
660
661	qLength = sqrt( qSqr );
662	q[0] = q[0] / qLength;
663	q[1] = q[1] / qLength;
664	q[2] = q[2] / qLength;
665	q[3] = q[3] / qLength;
666
667	dAtom->setQ( q );
668
669	// add the angular velocities
670
671	dAtom->setJx( jx );
672	dAtom->setJy( jy );
673	dAtom->setJz( jz );
674	}
675
676	// add the positions and velocities to the atom
677
678	atoms[atomIndex]->setX( rx );
679	atoms[atomIndex]->setY( ry );
680	atoms[atomIndex]->setZ( rz );
681
682	atoms[atomIndex]->set_vx( vx );
683	atoms[atomIndex]->set_vy( vy );
684	atoms[atomIndex]->set_vz( vz );
685
686	return NULL;
687	}
688
689
690	char* DumpReader::parseCommentLine(char* readLine, double time,
691	double boxMat[9]){
692
693	char *foo; // the pointer to the current string token
694	int j;
695
696	// set the string tokenizer
697
698	foo = strtok(readLine, " ,;\t");
699	if(foo == NULL){
700	sprintf( painCave.errMsg,
701	"error in reading time from %s\n",
702	inName );
703	return strdup( painCave.errMsg );
704	}
705	time = atof( foo );
706
707	// get the Hx vector
708
709	foo = strtok(NULL, " ,;\t");
710	if(foo == NULL){
711	sprintf( painCave.errMsg,
712	"error in reading Hx[0] from %s\n",
713	inName );
714	return strdup( painCave.errMsg );
715	}
716	boxMat[0] = atof( foo );
717
718	foo = strtok(NULL, " ,;\t");
719	if(foo == NULL){
720	sprintf( painCave.errMsg,
721	"error in reading Hx[1] from %s\n",
722	inName );
723	return strdup( painCave.errMsg );
724	}
725	boxMat[1] = atof( foo );
726
727	foo = strtok(NULL, " ,;\t");
728	if(foo == NULL){
729	sprintf( painCave.errMsg,
730	"error in reading Hx[2] from %s\n",
731	inName );
732	return strdup( painCave.errMsg );
733	}
734	boxMat[2] = atof( foo );
735
736	// get the Hy vector
737
738	foo = strtok(NULL, " ,;\t");
739	if(foo == NULL){
740	sprintf( painCave.errMsg,
741	"error in reading Hy[0] from %s\n",
742	inName );
743	return strdup( painCave.errMsg );
744	}
745	boxMat[3] = atof( foo );
746
747	foo = strtok(NULL, " ,;\t");
748	if(foo == NULL){
749	sprintf( painCave.errMsg,
750	"error in reading Hy[1] from %s\n",
751	inName );
752	return strdup( painCave.errMsg );
753	}
754	boxMat[4] = atof( foo );
755
756	foo = strtok(NULL, " ,;\t");
757	if(foo == NULL){
758	sprintf( painCave.errMsg,
759	"error in reading Hy[2] from %s\n",
760	inName );
761	return strdup( painCave.errMsg );
762	}
763	boxMat[5] = atof( foo );
764
765	// get the Hz vector
766
767	foo = strtok(NULL, " ,;\t");
768	if(foo == NULL){
769	sprintf( painCave.errMsg,
770	"error in reading Hz[0] from %s\n",
771	inName );
772	return strdup( painCave.errMsg );
773	}
774	boxMat[6] = atof( foo );
775
776	foo = strtok(NULL, " ,;\t");
777	if(foo == NULL){
778	sprintf( painCave.errMsg,
779	"error in reading Hz[1] from %s\n",
780	inName );
781	return strdup( painCave.errMsg );
782	}
783	boxMat[7] = atof( foo );
784
785	foo = strtok(NULL, " ,;\t");
786	if(foo == NULL){
787	sprintf( painCave.errMsg,
788	"error in reading Hz[2] from %s\n",
789	inName );
790	return strdup( painCave.errMsg );
791	}
792	boxMat[8] = atof( foo );
793
794	return NULL;
795	}
796
797
798	#ifdef IS_MPI
799
800	// a couple of functions to let us escape the read loop
801
802	void initFile::nodeZeroError( void ){
803	int j, myStatus;
804
805	myStatus = 0;
806	for (j = 0; j < mpiSim->getNumberProcessors(); j++) {
807	MPI_Send( &myStatus, 1, MPI_INT, j,
808	TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
809	}
810
811
812	MPI_Finalize();
813	exit (0);
814
815	}
816
817	void initFile::anonymousNodeDie( void ){
818
819	MPI_Finalize();
820	exit (0);
821	}
822
823	#endif //is_mpi