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 i, j;
  int lineNum = 0;
  char readBuffer[2000];
  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 ), 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;
  unsigned int 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->getNMolGlobal(); 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->getNMolGlobal(); 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]);

    if (fabs(qSqr) < 1e-6) {
      sprintf(painCave.errMsg,
          "initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0).\n");
       return strdup(painCave.errMsg);
    }

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