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;

  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();
    }

    nTotObjs = 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 != simnfo->getTotIntegrableObjects()){
      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[i]);
           
       }


      }
      else{
      //molecule belongs to slave nodes

        MPI_Recv(&nCurObj, 1, MPI_INT, 0,
               TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
      
       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();

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

      }
      
    }
    
  }
  else{
  //actions taken at slave nodes
    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_Recv(&nCurObj, 1, MPI_INT, 0,
                        TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);

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