OOPSE/libmdtools/DumpReader.cpp

#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 "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


DumpReader :: DumpReader(const 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();
  }
  
  inFileName = 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
  vector<fpos_t*>::iterator i;

  int error;
  error = fclose( inFile );
  if( error ){
    sprintf( painCave.errMsg,
             "Error closing %s\n", inFileName.c_str());
    simError();
  }

  for(i = framePos.begin(); i != framePos.end(); ++i)
    delete *i;
  framePos.clear();
  
#ifdef IS_MPI
  }
  strcpy( checkPointMsg, "Dump file closed successfully." );
  MPIcheckPoint();
#endif

  return;
}

int DumpReader::getNframes( void ){

  if( !isScanned )
    scanFile();
  return framePos.size();
}

void DumpReader::scanFile( void ){

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


#ifdef IS_MPI
  if( worldRank == 0 ){
#endif // is_mpi
    
    rewind( inFile );
    
    currPos = new fpos_t;
    fgetpos( inFile, currPos );
    fgets( readBuffer, sizeof( readBuffer ), inFile );
    lineNum++;
    if( feof( inFile ) ){
      sprintf( painCave.errMsg,
               "File \"%s\" ended unexpectedly at line %d\n", 
               inFileName.c_str(),
               lineNum );
      painCave.isFatal = 1;
      simError();
    }

    while( !feof( inFile ) ){
      
      framePos.push_back(currPos);

      i = atoi(readBuffer);
      
      fgets( readBuffer, sizeof( readBuffer ), inFile );
      lineNum++;    
      if( feof( inFile ) ){
        sprintf( painCave.errMsg,
                 "File \"%s\" ended unexpectedly at line %d\n",
                 inFileName.c_str(),
                 lineNum );
        painCave.isFatal = 1;
        simError();
      }
            
      for(j=0; j<i; j++){
        
        fgets( readBuffer, sizeof( readBuffer ), inFile );
        lineNum++;    
        if( feof( inFile ) ){
          sprintf( painCave.errMsg,
                   "File \"%s\" ended unexpectedly at line %d,"
                   " with atom %d\n", 
                   inFileName.c_str(),
                   lineNum, 
                   j );
          painCave.isFatal = 1;
          simError();
        }
        
      }
      
      currPos = new fpos_t;
      fgetpos( inFile, currPos );
      fgets( readBuffer, sizeof( readBuffer ), inFile );
      lineNum++;
    }
    
    delete currPos;
    rewind( inFile );
     
    isScanned = true;

#ifdef IS_MPI
  }
  strcpy( checkPointMsg, "Successfully scanned DumpFile\n" );
  MPIcheckPoint();
#endif // is_mpi
}

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

  simnfo = the_simnfo;

  this->readSet( whichFrame );
}


void DumpReader :: readSet( int whichFrame ){

   int i, j;

#ifdef IS_MPI
  int done, which_node, which_atom; // loop counter
#endif //is_mpi

  const int BUFFERSIZE = 2000; // size of the read buffer
  int nTotObjs; // the number of atoms
  char read_buffer[BUFFERSIZE]; //the line buffer for reading

  char *eof_test; // ptr to see when we reach the end of the file
  char *parseErr;

  vector<StuntDouble*> integrableObjects;


#ifndef IS_MPI

  fsetpos(inFile, framePos[whichFrame]);
  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",
             inFileName.c_str() );
    painCave.isFatal = 1;
    simError();
  }

  nTotObjs = atoi( read_buffer );

  if( nTotObjs != simnfo->getTotIntegrableObjects() ){
    sprintf( painCave.errMsg,
             "DumpReader error. %s n_atoms, %d, "
             "does not match the BASS file's n_atoms, %d.\n",
             inFileName.c_str(), nTotObjs, simnfo->getTotIntegrableObjects());
    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", inFileName.c_str());
    painCave.isFatal = 1;
    simError();
  }

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

  //parse dump lines

  for( i=0; i < simnfo->n_mol; i++){

    integrableObjects = (simnfo->molecules[i]).getIntegrableObjects();

    for(j = 0; j < integrableObjects.size(); j++){

      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",
              inFileName.c_str(), nTotObjs, i );
        painCave.isFatal = 1;
        simError();
      }
      
      parseErr = parseDumpLine( read_buffer, integrableObjects[j]);
      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 *MolToProcMap = mpiSim->getMolToProcMap();
  int localIndex;
  int nCurObj;
  int nitems;

  nTotObjs = simnfo->getTotIntegrableObjects();
  haveError = 0;
  if (worldRank == 0) {
     fsetpos(inFile,  framePos[whichFrame]);

    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
    if( eof_test == NULL ){
      sprintf( painCave.errMsg,
               "Error reading 1st line of %s \n ",inFileName.c_str());
      haveError = 1;
      simError();
    }

    nitems = atoi( read_buffer );

    // Check to see that the number of integrable objects  in the intial configuration file is the
    // same as declared in simBass.

    if( nTotObjs != nitems){
      sprintf( painCave.errMsg,
               "DumpReadererror. %s n_atoms, %d, "
               "does not match the BASS file's n_atoms, %d.\n",
               inFileName.c_str(), nTotObjs, simnfo->getTotIntegrableObjects());
      haveError= 1;
      simError();
    }

    //read the 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", inFileName.c_str());
      haveError = 1;
      simError();
    }

    //Every single processor will parse the comment line by itself
    //By using this way, we might lose some efficiency, but if we want to add
    //more parameters into comment line, we only need to modify function
    //parseCommentLine

    MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, 0, MPI_COMM_WORLD);

    parseErr = parseCommentLine( read_buffer, simnfo);

    if( parseErr != NULL ){
      strcpy( painCave.errMsg, parseErr );
      haveError = 1;
      simError();
    }

    for (i=0 ; i < mpiSim->getTotNmol(); i++) {
      which_node = MolToProcMap[i];
      if(which_node == 0){
       //molecules belong to master node

      localIndex = mpiSim->getGlobalToLocalMol(i);

      if(localIndex == -1) {
        strcpy(painCave.errMsg, "Molecule not found on node 0!");
        haveError = 1;
        simError();
      }

       integrableObjects = (simnfo->molecules[localIndex]).getIntegrableObjects();
       for(j=0; j < integrableObjects.size(); j++){
        
          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",
                    inFileName.c_str(), nTotObjs, i );
                haveError= 1;
                simError();
          }
          
          if(haveError) nodeZeroError();

          parseDumpLine(read_buffer, integrableObjects[j]);
           
       }


      }
      else{
      //molecule belongs to slave nodes

        MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
               TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
      
       for(j=0; j < nCurObj; j++){
        
          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",
                    inFileName.c_str(), nTotObjs, i );
                haveError= 1;
                simError();
          }
          
          if(haveError) nodeZeroError();

            MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
                      TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
           
       }

      }
      
    }
    
  }
  else{
  //actions taken at slave nodes
    MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, 0, MPI_COMM_WORLD);

    parseErr = parseCommentLine( read_buffer, simnfo);

    if( parseErr != NULL ){
      strcpy( painCave.errMsg, parseErr );
      haveError = 1;
      simError();
    }
  
    for (i=0 ; i < mpiSim->getTotNmol(); i++) {
      which_node = MolToProcMap[i];
      
      if(which_node == worldRank){
      //molecule with global index i belongs to this processor
      
        localIndex = mpiSim->getGlobalToLocalMol(i);

        if(localIndex == -1) {
          sprintf(painCave.errMsg, "Molecule not found on node %d\n", worldRank);
          haveError = 1;
          simError();
        }

        integrableObjects = (simnfo->molecules[localIndex]).getIntegrableObjects();        

        nCurObj = integrableObjects.size();
        
        MPI_Send(&nCurObj, 1, MPI_INT, 0,
                        TAKE_THIS_TAG_INT, MPI_COMM_WORLD);

        for(j = 0; j < integrableObjects.size(); j++){

          MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, 0,
                              TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);

          parseErr = parseDumpLine(read_buffer, integrableObjects[j]);

          if( parseErr != NULL ){
                strcpy( painCave.errMsg, parseErr );
                simError();
          }

        }
          
      }
      
    }

  }

#endif
}

char* DumpReader::parseDumpLine(char* readLine, StuntDouble* sd){

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

  double pos[3]; // position place holders
  double vel[3]; // velocity placeholders
  double q[4]; // the quaternions
  double ji[3]; // angular velocity placeholders;
  double qSqr, qLength; // needed to normalize the quaternion vector.


  // set the string tokenizer

  foo = strtok(readLine, " ,;\t");

  // check the atom name to the current atom

  if( strcmp( foo, sd->getType() ) ){
    sprintf( painCave.errMsg,
             "DumpReader error.  Does not"
             " match the BASS atom %s.\n",
             sd->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",
             inFileName.c_str());
    return strdup( painCave.errMsg );
  }
  pos[0] = atof( foo );

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading postition y from %s\n",
             inFileName.c_str());
    return strdup( painCave.errMsg );
  }
  pos[1] = atof( foo );

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading postition z from %s\n",
             inFileName.c_str());
    return strdup( painCave.errMsg );
  }
  pos[2] = atof( foo );


  // get the velocities

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading velocity x from %s\n",
             inFileName.c_str() );
    return strdup( painCave.errMsg );
  }
  vel[0] = atof( foo );

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading velocity x from %s\n",
             inFileName.c_str() );
    return strdup( painCave.errMsg );
  }
  vel[1] = atof( foo );

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading velocity x from %s\n",
             inFileName.c_str() );
    return strdup( painCave.errMsg );
  }
  vel[2] = atof( foo );


  // add the positions and velocities to the atom

  sd->setPos( pos );
  sd->setVel( vel );

  if (!sd->isDirectional())
    return NULL;

  // get the quaternions

  if( sd->isDirectional() ){

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
                     "error in reading velocity x from %s\n",
                      inFileName.c_str() );
      return strdup( painCave.errMsg );
    }
    q[0] = atof( foo );

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
                     "error in reading velocity x from %s\n",
                      inFileName.c_str() );
      return strdup( painCave.errMsg );
    }
    q[1] = atof( foo );

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
                     "error in reading velocity x from %s\n",
                      inFileName.c_str() );
      return strdup( painCave.errMsg );
    }
    q[2] = atof( foo );

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
                     "error in reading velocity x from %s\n",
                      inFileName.c_str() );
      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 velocity x from %s\n",
                      inFileName.c_str() );
      return strdup( painCave.errMsg );
    }
    ji[0] = atof( foo );

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
                     "error in reading velocity x from %s\n",
                      inFileName.c_str() );
      return strdup( painCave.errMsg );
    }
    ji[1] = atof(foo );

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
                     "error in reading velocity x from %s\n",
                      inFileName.c_str() );
      return strdup( painCave.errMsg );
    }
    ji[2] = atof( foo );


    // 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;

    // add quaternion and angular velocities

    sd->setQ( q );
    sd->setJ( ji );
  }


  return NULL;
}


char* DumpReader::parseCommentLine(char* readLine, SimInfo* entry_plug){

  double currTime;
  double boxMat[9];
  double theBoxMat3[3][3];
  double chi;
  double integralOfChidt;
  double eta[9];

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

  // set the string tokenizer

  foo = strtok(readLine, " ,;\t");
  // set the timeToken.

  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading Time from %s\n",
             inFileName.c_str() );
    return strdup( painCave.errMsg );
  }

  currTime = atof( foo );
  entry_plug->setTime( currTime );

  //get H-Matrix

  for(int i = 0 ; i < 9; i++){
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
               "error in reading H[%d] from %s\n", i, inFileName.c_str() );
      return strdup( painCave.errMsg );
    }
    boxMat[i] = atof( foo );
  }

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

  //set H-Matrix
  entry_plug->setBoxM( theBoxMat3 );

  //get chi and integralOfChidt, they should appear by pair

  if( entry_plug->useInitXSstate ){
    foo = strtok(NULL, " ,;\t\n");
    if(foo != NULL){
      chi = atof(foo);
      
      foo = strtok(NULL, " ,;\t\n");
      if(foo == NULL){
        sprintf( painCave.errMsg,
                 "chi and integralOfChidt should appear by pair in %s\n", inFileName.c_str() );
        return strdup( painCave.errMsg );
      }
      integralOfChidt = atof( foo );
      
      //push chi and integralOfChidt into SimInfo::properties which can be
      //retrieved by integrator later
      DoubleData* chiValue = new DoubleData();
      chiValue->setID(CHIVALUE_ID);
      chiValue->setData(chi);
      entry_plug->addProperty(chiValue);
      
      DoubleData* integralOfChidtValue = new DoubleData();
      integralOfChidtValue->setID(INTEGRALOFCHIDT_ID);
      integralOfChidtValue->setData(integralOfChidt);
      entry_plug->addProperty(integralOfChidtValue);
      
    }
    else
      return NULL;
    
    //get eta
    foo = strtok(NULL, " ,;\t\n");
    if(foo != NULL ){
  
      for(int i = 0 ; i < 9; i++){
        
        if(foo == NULL){
          sprintf( painCave.errMsg,
                   "error in reading eta[%d] from %s\n", i, inFileName.c_str() );
          return strdup( painCave.errMsg );
        }
        eta[i] = atof( foo );
        foo = strtok(NULL, " ,;\t\n");
      }
    }
    else
      return NULL;
    
    //push eta into SimInfo::properties which can be
    //retrieved by integrator later
    //entry_plug->setBoxM( theBoxMat3 );
    DoubleArrayData* etaValue = new DoubleArrayData();
    etaValue->setID(ETAVALUE_ID);
    etaValue->setData(eta, 9);
    entry_plug->addProperty(etaValue);
  }

  return NULL;
}

#ifdef IS_MPI
void DumpReader::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 DumpReader::anonymousNodeDie( void ){

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