OOPSE/libmdtools/InitializeFromFile.cpp

#include <iostream>
#include <cmath>

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

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

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

namespace initFile{
  void nodeZeroError( void );
  void anonymousNodeDie( void );
}

using namespace initFile;

#endif // is_mpi

InitializeFromFile :: InitializeFromFile( char *in_name ){
#ifdef IS_MPI
  if (worldRank == 0) {
#endif

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

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

  return;
}


void InitializeFromFile :: read_xyz( SimInfo* the_simnfo ){

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

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

  char *eof_test; // ptr to see when we reach the end of the file 
  char *parseErr;
  int procIndex;
  double boxMat[9];
  

  simnfo = the_simnfo;


#ifndef IS_MPI
  eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
  if( eof_test == NULL ){
    sprintf( painCave.errMsg,
             "InitializeFromFile error: error reading 1st line of \"%s\"\n",
             c_in_name );
    painCave.isFatal = 1;
    simError();
  }

  n_atoms = atoi( read_buffer );

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

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

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

  simnfo->setBoxM( boxMat );


  for( i=0; i < n_atoms; i++){
    
    eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
    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",
              c_in_name, n_atoms, i );
      painCave.isFatal = 1;
      simError();
    }

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


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

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

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

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

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

    MPI_Bcast(boxMat, 9, MPI_DOUBLE, 0, MPI_COMM_WORLD );
  
    if(haveError) nodeZeroError();

    for (i=0 ; i < mpiSim->getTotAtoms(); i++) {
      
      eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
      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",
                c_in_name, n_atoms, i );
        haveError= 1;
        simError();
      }
    
      if(haveError) nodeZeroError();

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

    done = 0;
    while (!done) {

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

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

   
#endif
}

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

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

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

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

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


  // get the velocities

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

    // check that the quaternion vector is normalized

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

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

  return NULL;
}


char* InitializeFromFile::parseBoxLine(char* readLine, double boxMat[9]){

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

  // set the string tokenizer
  
  foo = strtok(readLine, " ,;\t");
  // ignore the first token which is the time stamp.


  // get the Hx vector

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

  // get the Hy vector

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

  // get the Hz vector

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

  return NULL;
}


#ifdef IS_MPI

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

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

  MPI_Finalize();
  exit (0);
  
}

void initFile::anonymousNodeDie( void ){

  MPI_Finalize();
  exit (0);
}

#endif //is_mpi
Revision:	586
Committed:	Wed Jul 9 22:14:06 2003 UTC (21 years ago) by mmeineke
File size:	16075 byte(s)
Log Message:	Bug fixing NPTi and NPTf. there is some error in the caclulation of HmatInverse.
#	User	Rev	Content
1	mmeineke	377	#include <iostream>
2			#include <cmath>
3
4			#include <stdio.h>
5			#include <stdlib.h>
6			#include <string.h>
7			#include <unistd.h>
8			#include <sys/types.h>
9			#include <sys/stat.h>
10
11			#include "ReadWrite.hpp"
12			#include "simError.h"
13
14			#ifdef IS_MPI
15	gezelter	417	#include <mpi.h>
16	mmeineke	377	#include "mpiSimulation.hpp"
17	chuckv	436	#define TAKE_THIS_TAG_CHAR 0
18			#define TAKE_THIS_TAG_INT 1
19
20	mmeineke	440	namespace initFile{
21			void nodeZeroError( void );
22			void anonymousNodeDie( void );
23			}
24	chuckv	436
25	mmeineke	440	using namespace initFile;
26
27	mmeineke	377	#endif // is_mpi
28
29			InitializeFromFile :: InitializeFromFile( char *in_name ){
30			#ifdef IS_MPI
31			if (worldRank == 0) {
32			#endif
33
34			c_in_file = fopen(in_name, "r");
35			if(c_in_file == NULL){
36			sprintf(painCave.errMsg,
37			"Cannot open file: %s\n", in_name);
38			painCave.isFatal = 1;
39			simError();
40			}
41
42			strcpy( c_in_name, in_name);
43			#ifdef IS_MPI
44			}
45			strcpy( checkPointMsg, "Infile opened for reading successfully." );
46			MPIcheckPoint();
47			#endif
48			return;
49			}
50
51			InitializeFromFile :: ~InitializeFromFile( ){
52			#ifdef IS_MPI
53			if (worldRank == 0) {
54			#endif
55			int error;
56			error = fclose( c_in_file );
57			if( error ){
58			sprintf( painCave.errMsg,
59			"Error closing %s\n", c_in_name );
60			simError();
61			}
62			#ifdef IS_MPI
63			}
64			strcpy( checkPointMsg, "Infile closed successfully." );
65			MPIcheckPoint();
66			#endif
67
68			return;
69			}
70
71
72	mmeineke	586	void InitializeFromFile :: read_xyz( SimInfo* the_simnfo ){
73	mmeineke	377
74	gezelter	417	int i, j, done, which_node, which_atom; // loop counter
75	mmeineke	377
76			const int BUFFERSIZE = 2000; // size of the read buffer
77			int n_atoms; // the number of atoms
78			char read_buffer[BUFFERSIZE]; //the line buffer for reading
79			#ifdef IS_MPI
80			char send_buffer[BUFFERSIZE];
81			#endif
82
83			char *eof_test; // ptr to see when we reach the end of the file
84			char *parseErr;
85			int procIndex;
86	mmeineke	586	double boxMat[9];
87
88	mmeineke	377
89	mmeineke	586	simnfo = the_simnfo;
90	mmeineke	377
91
92			#ifndef IS_MPI
93			eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
94			if( eof_test == NULL ){
95			sprintf( painCave.errMsg,
96			"InitializeFromFile error: error reading 1st line of \"%s\"\n",
97			c_in_name );
98			painCave.isFatal = 1;
99			simError();
100			}
101
102			n_atoms = atoi( read_buffer );
103
104	mmeineke	586	Atom **atoms = simnfo->atoms;
105	mmeineke	377	DirectionalAtom* dAtom;
106
107	mmeineke	586	if( n_atoms != simnfo->n_atoms ){
108	mmeineke	377	sprintf( painCave.errMsg,
109			"Initialize from File error. %s n_atoms, %d, "
110			"does not match the BASS file's n_atoms, %d.\n",
111	mmeineke	586	c_in_name, n_atoms, simnfo->n_atoms );
112	mmeineke	377	painCave.isFatal = 1;
113			simError();
114			}
115
116	mmeineke	586	//read the box mat from the comment line
117	mmeineke	377
118			eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
119			if(eof_test == NULL){
120			sprintf( painCave.errMsg,
121			"error in reading commment in %s\n", c_in_name);
122			painCave.isFatal = 1;
123			simError();
124			}
125
126	mmeineke	586	parseErr = parseBoxLine( read_buffer, boxMat );
127			if( parseErr != NULL ){
128			strcpy( painCave.errMsg, parseErr );
129			painCave.isFatal = 1;
130			simError();
131			}
132
133			simnfo->setBoxM( boxMat );
134
135
136	mmeineke	377	for( i=0; i < n_atoms; i++){
137
138			eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
139			if(eof_test == NULL){
140			sprintf(painCave.errMsg,
141			"error in reading file %s\n"
142			"natoms = %d; index = %d\n"
143			"error reading the line from the file.\n",
144			c_in_name, n_atoms, i );
145			painCave.isFatal = 1;
146			simError();
147			}
148
149
150			parseErr = parseDumpLine( read_buffer, i );
151			if( parseErr != NULL ){
152			strcpy( painCave.errMsg, parseErr );
153			painCave.isFatal = 1;
154			simError();
155	mmeineke	586	}
156	mmeineke	377	}
157
158
159			// MPI Section of code..........
160			#else //IS_MPI
161
162	chuckv	436	// first thing first, suspend fatalities.
163			painCave.isEventLoop = 1;
164
165			int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
166			int haveError;
167
168	mmeineke	447	MPI_Status istatus;
169	gezelter	417	int *AtomToProcMap = mpiSim->getAtomToProcMap();
170	mmeineke	377
171	chuckv	436
172			haveError = 0;
173	mmeineke	377	if (worldRank == 0) {
174	chuckv	436
175	mmeineke	377	eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
176			if( eof_test == NULL ){
177			sprintf( painCave.errMsg,
178			"Error reading 1st line of %d \n ",c_in_name);
179	chuckv	436	haveError = 1;
180	mmeineke	377	simError();
181			}
182
183			n_atoms = atoi( read_buffer );
184
185	mmeineke	586	Atom **atoms = simnfo->atoms;
186	mmeineke	377	DirectionalAtom* dAtom;
187
188			// Check to see that the number of atoms in the intial configuration file is the
189			// same as declared in simBass.
190
191			if( n_atoms != mpiSim->getTotAtoms() ){
192			sprintf( painCave.errMsg,
193			"Initialize from File error. %s n_atoms, %d, "
194			"does not match the BASS file's n_atoms, %d.\n",
195	mmeineke	586	c_in_name, n_atoms, simnfo->n_atoms );
196	chuckv	436	haveError= 1;
197	mmeineke	377	simError();
198			}
199
200	mmeineke	586	//read the boxMat from the comment line
201	mmeineke	377
202			eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
203			if(eof_test == NULL){
204			sprintf( painCave.errMsg,
205			"error in reading commment in %s\n", c_in_name);
206	mmeineke	586	haveError = 1;
207	mmeineke	377	simError();
208			}
209	mmeineke	586
210			parseErr = parseBoxLine( read_buffer, boxMat );
211			if( parseErr != NULL ){
212			strcpy( painCave.errMsg, parseErr );
213			haveError = 1;
214			simError();
215			}
216
217			MPI_Bcast(boxMat, 9, MPI_DOUBLE, 0, MPI_COMM_WORLD );
218	mmeineke	377
219	chuckv	436	if(haveError) nodeZeroError();
220	gezelter	417
221			for (i=0 ; i < mpiSim->getTotAtoms(); i++) {
222
223	mmeineke	377	eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
224			if(eof_test == NULL){
225			sprintf(painCave.errMsg,
226			"error in reading file %s\n"
227			"natoms = %d; index = %d\n"
228			"error reading the line from the file.\n",
229			c_in_name, n_atoms, i );
230	chuckv	436	haveError= 1;
231	mmeineke	377	simError();
232			}
233
234	chuckv	436	if(haveError) nodeZeroError();
235
236	gezelter	417	// Get the Node number which wants this atom:
237			which_node = AtomToProcMap[i];
238	chuckv	436	if (which_node == 0) {
239	gezelter	417	parseErr = parseDumpLine( read_buffer, i );
240			if( parseErr != NULL ){
241			strcpy( painCave.errMsg, parseErr );
242	chuckv	436	haveError = 1;
243	gezelter	417	simError();
244			}
245	chuckv	436	if(haveError) nodeZeroError();
246			}
247
248			else {
249
250			myStatus = 1;
251	mmeineke	447	MPI_Send(&myStatus, 1, MPI_INT, which_node,
252			TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
253			MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
254			TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
255			MPI_Send(&i, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
256			MPI_COMM_WORLD);
257			MPI_Recv(&myStatus, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
258			MPI_COMM_WORLD, &istatus);
259	chuckv	436
260			if(!myStatus) nodeZeroError();
261	mmeineke	377	}
262			}
263	chuckv	436	myStatus = -1;
264	gezelter	417	for (j = 0; j < mpiSim->getNumberProcessors(); j++) {
265	mmeineke	447	MPI_Send( &myStatus, 1, MPI_INT, j,
266			TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
267	gezelter	417	}
268	gezelter	419
269	gezelter	417	} else {
270	gezelter	419
271	mmeineke	586	MPI_Bcast(boxMat, 9, MPI_DOUBLE, 0, MPI_COMM_WORLD);
272
273	gezelter	417	done = 0;
274			while (!done) {
275	chuckv	436
276	mmeineke	447	MPI_Recv(&myStatus, 1, MPI_INT, 0,
277			TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
278	chuckv	436
279			if(!myStatus) anonymousNodeDie();
280
281			if(myStatus < 0) break;
282
283	mmeineke	447	MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, 0,
284			TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
285			MPI_Recv(&which_atom, 1, MPI_INT, 0,
286			TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
287	chuckv	436
288			myStatus = 1;
289			parseErr = parseDumpLine( read_buffer, which_atom );
290			if( parseErr != NULL ){
291			strcpy( painCave.errMsg, parseErr );
292			myStatus = 0;;
293			simError();
294	mmeineke	377	}
295	chuckv	436
296	mmeineke	447	MPI_Send( &myStatus, 1, MPI_INT, 0,
297			TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
298	chuckv	436
299	mmeineke	377	}
300			}
301	chuckv	436
302			// last thing last, enable fatalities.
303			painCave.isEventLoop = 0;
304	mmeineke	586	simnfo->setBoxM( boxMat );
305
306	chuckv	436
307	mmeineke	377	#endif
308			}
309
310	gezelter	419	char* InitializeFromFile::parseDumpLine(char* readLine, int globalIndex){
311	mmeineke	377
312			char *foo; // the pointer to the current string token
313
314			double rx, ry, rz; // position place holders
315			double vx, vy, vz; // velocity placeholders
316			double q[4]; // the quaternions
317			double jx, jy, jz; // angular velocity placeholders;
318			double qSqr, qLength; // needed to normalize the quaternion vector.
319
320	mmeineke	586	Atom **atoms = simnfo->atoms;
321	mmeineke	377	DirectionalAtom* dAtom;
322
323	gezelter	419	int j, n_atoms, atomIndex;
324	mmeineke	377
325			#ifdef IS_MPI
326			n_atoms = mpiSim->getTotAtoms();
327	gezelter	419	atomIndex=-1;
328			for (j=0; j < mpiSim->getMyNlocal(); j++) {
329			if (atoms[j]->getGlobalIndex() == globalIndex) atomIndex = j;
330			}
331			if (atomIndex == -1) {
332			sprintf( painCave.errMsg,
333			"Initialize from file error. Atom at index %d "
334			"in file %s does not exist on processor %d .\n",
335			globalIndex, c_in_name, mpiSim->getMyNode() );
336			return strdup( painCave.errMsg );
337			}
338	mmeineke	377	#else
339	mmeineke	586	n_atoms = simnfo->n_atoms;
340	gezelter	419	atomIndex = globalIndex;
341	mmeineke	377	#endif // is_mpi
342
343			// set the string tokenizer
344
345			foo = strtok(readLine, " ,;\t");
346
347			// check the atom name to the current atom
348
349			if( strcmp( foo, atoms[atomIndex]->getType() ) ){
350			sprintf( painCave.errMsg,
351			"Initialize from file error. Atom %s at index %d "
352			"in file %s does not"
353			" match the BASS atom %s.\n",
354			foo, atomIndex, c_in_name, atoms[atomIndex]->getType() );
355			return strdup( painCave.errMsg );
356			}
357
358			// get the positions
359
360			foo = strtok(NULL, " ,;\t");
361			if(foo == NULL){
362			sprintf( painCave.errMsg,
363			"error in reading postition x from %s\n"
364			"natoms = %d, index = %d\n",
365			c_in_name, n_atoms, atomIndex );
366			return strdup( painCave.errMsg );
367			}
368			rx = atof( foo );
369
370			foo = strtok(NULL, " ,;\t");
371			if(foo == NULL){
372			sprintf( painCave.errMsg,
373			"error in reading postition y from %s\n"
374			"natoms = %d, index = %d\n",
375			c_in_name, n_atoms, atomIndex );
376			return strdup( painCave.errMsg );
377			}
378			ry = atof( foo );
379
380			foo = strtok(NULL, " ,;\t");
381			if(foo == NULL){
382			sprintf( painCave.errMsg,
383			"error in reading postition z from %s\n"
384			"natoms = %d, index = %d\n",
385			c_in_name, n_atoms, atomIndex );
386			return strdup( painCave.errMsg );
387			}
388			rz = atof( foo );
389
390
391			// get the velocities
392
393			foo = strtok(NULL, " ,;\t");
394			if(foo == NULL){
395			sprintf( painCave.errMsg,
396			"error in reading velocity x from %s\n"
397			"natoms = %d, index = %d\n",
398			c_in_name, n_atoms, atomIndex );
399			return strdup( painCave.errMsg );
400			}
401			vx = atof( foo );
402
403			foo = strtok(NULL, " ,;\t");
404			if(foo == NULL){
405			sprintf( painCave.errMsg,
406			"error in reading velocity y from %s\n"
407			"natoms = %d, index = %d\n",
408			c_in_name, n_atoms, atomIndex );
409			return strdup( painCave.errMsg );
410			}
411			vy = atof( foo );
412
413			foo = strtok(NULL, " ,;\t");
414			if(foo == NULL){
415			sprintf( painCave.errMsg,
416			"error in reading velocity z from %s\n"
417			"natoms = %d, index = %d\n",
418			c_in_name, n_atoms, atomIndex );
419			return strdup( painCave.errMsg );
420			}
421			vz = atof( foo );
422
423
424			// get the quaternions
425
426			if( atoms[atomIndex]->isDirectional() ){
427
428			foo = strtok(NULL, " ,;\t");
429			if(foo == NULL){
430			sprintf(painCave.errMsg,
431			"error in reading quaternion 0 from %s\n"
432			"natoms = %d, index = %d\n",
433			c_in_name, n_atoms, atomIndex );
434			return strdup( painCave.errMsg );
435			}
436			q[0] = atof( foo );
437
438			foo = strtok(NULL, " ,;\t");
439			if(foo == NULL){
440			sprintf( painCave.errMsg,
441			"error in reading quaternion 1 from %s\n"
442			"natoms = %d, index = %d\n",
443			c_in_name, n_atoms, atomIndex );
444			return strdup( painCave.errMsg );
445			}
446			q[1] = atof( foo );
447
448			foo = strtok(NULL, " ,;\t");
449			if(foo == NULL){
450			sprintf( painCave.errMsg,
451			"error in reading quaternion 2 from %s\n"
452			"natoms = %d, index = %d\n",
453			c_in_name, n_atoms, atomIndex );
454			return strdup( painCave.errMsg );
455			}
456			q[2] = atof( foo );
457
458			foo = strtok(NULL, " ,;\t");
459			if(foo == NULL){
460			sprintf( painCave.errMsg,
461			"error in reading quaternion 3 from %s\n"
462			"natoms = %d, index = %d\n",
463			c_in_name, n_atoms, atomIndex );
464			return strdup( painCave.errMsg );
465			}
466			q[3] = atof( foo );
467
468			// get the angular velocities
469
470			foo = strtok(NULL, " ,;\t");
471			if(foo == NULL){
472			sprintf( painCave.errMsg,
473			"error in reading angular momentum jx from %s\n"
474			"natoms = %d, index = %d\n",
475			c_in_name, n_atoms, atomIndex );
476			return strdup( painCave.errMsg );
477			}
478			jx = atof( foo );
479
480			foo = strtok(NULL, " ,;\t");
481			if(foo == NULL){
482			sprintf( painCave.errMsg,
483			"error in reading angular momentum jy from %s\n"
484			"natoms = %d, index = %d\n",
485			c_in_name, n_atoms, atomIndex );
486			return strdup( painCave.errMsg );
487			}
488			jy = atof(foo );
489
490			foo = strtok(NULL, " ,;\t");
491			if(foo == NULL){
492			sprintf( painCave.errMsg,
493			"error in reading angular momentum jz from %s\n"
494			"natoms = %d, index = %d\n",
495			c_in_name, n_atoms, atomIndex );
496			return strdup( painCave.errMsg );
497			}
498			jz = atof( foo );
499
500			dAtom = ( DirectionalAtom* )atoms[atomIndex];
501
502			// check that the quaternion vector is normalized
503
504			qSqr = (q[0] * q[0]) + (q[1] * q[1]) + (q[2] * q[2]) + (q[3] * q[3]);
505
506			qLength = sqrt( qSqr );
507			q[0] = q[0] / qLength;
508			q[1] = q[1] / qLength;
509			q[2] = q[2] / qLength;
510			q[3] = q[3] / qLength;
511
512			dAtom->setQ( q );
513
514			// add the angular velocities
515
516			dAtom->setJx( jx );
517			dAtom->setJy( jy );
518			dAtom->setJz( jz );
519			}
520
521			// add the positions and velocities to the atom
522
523			atoms[atomIndex]->setX( rx );
524			atoms[atomIndex]->setY( ry );
525			atoms[atomIndex]->setZ( rz );
526
527			atoms[atomIndex]->set_vx( vx );
528			atoms[atomIndex]->set_vy( vy );
529			atoms[atomIndex]->set_vz( vz );
530
531			return NULL;
532			}
533	chuckv	436
534
535	mmeineke	586	char* InitializeFromFile::parseBoxLine(char* readLine, double boxMat[9]){
536
537			char *foo; // the pointer to the current string token
538			int j;
539
540			// set the string tokenizer
541
542			foo = strtok(readLine, " ,;\t");
543			// ignore the first token which is the time stamp.
544
545
546			// get the Hx vector
547
548			foo = strtok(NULL, " ,;\t");
549			if(foo == NULL){
550			sprintf( painCave.errMsg,
551			"error in reading Hx[0] from %s\n",
552			c_in_name );
553			return strdup( painCave.errMsg );
554			}
555			boxMat[0] = atof( foo );
556
557			foo = strtok(NULL, " ,;\t");
558			if(foo == NULL){
559			sprintf( painCave.errMsg,
560			"error in reading Hx[1] from %s\n",
561			c_in_name );
562			return strdup( painCave.errMsg );
563			}
564			boxMat[1] = atof( foo );
565
566			foo = strtok(NULL, " ,;\t");
567			if(foo == NULL){
568			sprintf( painCave.errMsg,
569			"error in reading Hx[2] from %s\n",
570			c_in_name );
571			return strdup( painCave.errMsg );
572			}
573			boxMat[2] = atof( foo );
574
575			// get the Hy vector
576
577			foo = strtok(NULL, " ,;\t");
578			if(foo == NULL){
579			sprintf( painCave.errMsg,
580			"error in reading Hy[0] from %s\n",
581			c_in_name );
582			return strdup( painCave.errMsg );
583			}
584			boxMat[3] = atof( foo );
585
586			foo = strtok(NULL, " ,;\t");
587			if(foo == NULL){
588			sprintf( painCave.errMsg,
589			"error in reading Hy[1] from %s\n",
590			c_in_name );
591			return strdup( painCave.errMsg );
592			}
593			boxMat[4] = atof( foo );
594
595			foo = strtok(NULL, " ,;\t");
596			if(foo == NULL){
597			sprintf( painCave.errMsg,
598			"error in reading Hy[2] from %s\n",
599			c_in_name );
600			return strdup( painCave.errMsg );
601			}
602			boxMat[5] = atof( foo );
603
604			// get the Hz vector
605
606			foo = strtok(NULL, " ,;\t");
607			if(foo == NULL){
608			sprintf( painCave.errMsg,
609			"error in reading Hz[0] from %s\n",
610			c_in_name );
611			return strdup( painCave.errMsg );
612			}
613			boxMat[6] = atof( foo );
614
615			foo = strtok(NULL, " ,;\t");
616			if(foo == NULL){
617			sprintf( painCave.errMsg,
618			"error in reading Hz[1] from %s\n",
619			c_in_name );
620			return strdup( painCave.errMsg );
621			}
622			boxMat[7] = atof( foo );
623
624			foo = strtok(NULL, " ,;\t");
625			if(foo == NULL){
626			sprintf( painCave.errMsg,
627			"error in reading Hz[2] from %s\n",
628			c_in_name );
629			return strdup( painCave.errMsg );
630			}
631			boxMat[8] = atof( foo );
632
633			return NULL;
634			}
635
636
637	chuckv	436	#ifdef IS_MPI
638
639			// a couple of functions to let us escape the read loop
640
641	mmeineke	440	void initFile::nodeZeroError( void ){
642	chuckv	436	int j, myStatus;
643
644			myStatus = 0;
645			for (j = 0; j < mpiSim->getNumberProcessors(); j++) {
646	mmeineke	447	MPI_Send( &myStatus, 1, MPI_INT, j,
647			TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
648	chuckv	436	}
649
650
651			MPI_Finalize();
652			exit (0);
653
654			}
655
656	mmeineke	440	void initFile::anonymousNodeDie( void ){
657	chuckv	436
658			MPI_Finalize();
659			exit (0);
660			}
661
662			#endif //is_mpi