trunk/tcProps/readWrite.c

#define _FILE_OFFSET_BITS 64

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

#include "params.h"
#include "tcProps.h"
#include "readWrite.h"


#define BUFFER_SIZE 2000
int isScanned;
int fileOpen;
int nFrames;
double *frameTimes;
FILE* inFile;


struct linkedPos{
  
  fpos_t *myPos;
  double timeStamp;
  double Hmat[3][3];
  struct linkedPos* next;
};

struct staticPos{
    
  fpos_t *myPos;
  double Hmat[3][3];
};
struct staticPos* posArray;

void parseDumpLine( char readBuffer[BUFFER_SIZE], int index, 
                    struct atomCoord* atoms );
void setU( double q[4], double u[3] );
void normalizeQ( double q[4] );

void openFile( char* inName ){
  
  inFile = fopen(inName);
  if(inFile ==NULL){
    fprintf(stderr,
            "Error opening file \"%s\"\n",
            inName);
    exit(0);
  }

  fileOpen = 1;
}  

void closeFile( void ){

  fclose( inFile );

  fileOpen = 0;
}

int setFrames( void ){

  int i,j,k;
  struct linkedPos* headPos;
  struct linkedPos* currPos;
  fpos_t *currPT;
  char readBuffer[BUFFER_SIZE];
  char trash[BUFFER_SIZE];
  char* foo;
  int lineNum = 0;
  
  
  if( !fileOpen ){
    fprintf(stderr,
            "Attempt to scan the file without first opening the file.\n" );
    exit(0);
  }


  nFrames = 0;
  headPos = (struct linkedPos*)malloc(sizeof(struct linkedPos));
  currPos = headPos;
  while( !feof( inFile ) ){

    currPT = (fpos_t *)malloc(sizeof(fpos_t));
    fgetpos(inFile, currPT);
    
    fgets( readBuffer, sizeof( readBuffer ), inFile );
    lineNum++;
    if( feof( inFile ) ){
      fprintf( stderr,
               "File \"%s\" ended unexpectedly at line %d\n", 
               inName,
               lineNum );
      exit(0);
    }
    
    currPos->next = (struct linkedPos*)malloc(sizeof(struct linkedPos));
    currPos = currPos->next;
    currPos->myPos = currPT;      
    nFrames++;

    i = atoi(readBuffer);

    fgets( readBuffer, sizeof( readBuffer ), inFile );
    lineNum++;
    if( feof( inFile ) ){
      fprintf( stderr,
               "File \"%s\" ended unexpectedly at line %d\n", 
               inName,
               lineNum );
      exit(0);
    }
    
    
    // parse the comment line
    
    foo = strtok(readBuffer, " ,;\t");
      
    if(foo == NULL){
      fprintf(stderr,
              "error in reading time from %s at line %d\n",
              inName, lineNum );
      exit(0);
    }
    currPos->timeStamp = atof( foo );
      
    // get the Hx vector
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      fprintf(stderr,
              "error in reading Hx[0] from %s at line %d\n",
              inName, lineNum );
      exit(0);
    }
    currPos->Hmat[0][0] = atof( foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      fprintf(stderr,
              "error in reading Hx[1] from %s at line %d\n",
              inName, lineNum );
      exit(0);
    }
    currPos->Hmat[1][0] = atof( foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      fprintf(stderr,
              "error in reading Hx[2] from %s at line %d\n",
              inName, lineNum );
      exit(0);
    }
    currPos->Hmat[2][0] = atof( foo );    
      
    // get the Hy vector
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      fprintf(stderr,
              "error in reading Hy[0] from %s at line %d\n",
              inName, lineNum );
      exit(0);
    }
    currPos->Hmat[0][1] = atof( foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      fprintf(stderr,
              "error in reading Hy[1] from %s at line %d\n",
              inName, lineNum );
      exit(0);
    }
    currPos->Hmat[1][1] = atof( foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      fprintf(stderr,
              "error in reading Hy[2] from %s at line %d\n",
              inName, lineNum );
      exit(0);
    }
    currPos->Hmat[2][1] = atof( foo );    
      
    // get the Hz vector
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      fprintf(stderr,
              "error in reading Hz[0] from %s at line %d\n",
              inName, lineNum );
      exit(0);
    }
    currPos->Hmat[0][2] = atof( foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      fprintf(stderr,
              "error in reading Hz[1] from %s at line %d\n",
              inName, lineNum );
      exit(0);
    }
    currPos->Hmat[1][2] = atof( foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      fprintf(stderr,
              "error in reading Hz[2] from %s at line %d\n",
              inName, lineNum );
      exit(0);
    }
    currPos->Hmat[2][2] = atof( foo );
    
    // skip the atoms

    for(j=0; j<i; j++){
      
      fgets( readBuffer, sizeof( readBuffer ), inFile );
      lineNum++;    
      if( feof( inFile ) ){
        fprintf( stderr,
                 "File \"%s\" ended unexpectedly at line %d,"
                 " with atom %d\n", 
                 inName,
                 lineNum, 
                 j );
        exit(0);
      }
    }
  }
  
  // allocate the static position array

  posArray = (struct staticPos*)calloc(nFrames, sizeof(struct staticPos));
  frameTimes = (double *)calloc(nFrames, sizeof(double) );

  i=0;
  currPos = headPos->next;
  while( currPos != NULL ){
    
    if(i>=nFrames){
      
      fprintf( stderr,
               "The number of frame pointers exceeded the number of frames.\n"
               "  OOPS!\n" );
      exit(0);
    }
    
    frameTimes[i] = currPos->timeStamp;
    posArray[i].myPos = currPos->myPos;
    for(j=0;j<3;j++){
      for(k=0;k<3;k++){
        posArray[i].Hmat[j][k] = currPos->Hmat[j][k];
      }
    }

    i++;
    currPos = currPos->next;
  }  

  // clean up the linked list

  currPos = headPos;
  while( currPos != NULL ){
    headPos = currPos;
    currPos = headPos->next;
    free(headPos);
  }


  isScanned = 1;
}


void readFrame(int theFrame, struct atomCoord* atoms, double Hmat[3][3] ){

  // list of 'a priori' constants

  const int nLipAtoms = NL_ATOMS;
  const int nBonds    = NBONDS;
  const int nLipids   = NLIPIDS;
  const int nSSD      = NSSD;
  const int nAtoms    = nLipAtoms * nLipids + nSSD;

  fpos_t *framePos;
  char readBuffer[BUFFER_SIZE];

  int i, j, k;


  // quick checks

  if (theFrame >= nFrames){
    
    fprintf(stderr,
            "frame %d, is out of range\n",
            theFrame );
    exit(0);
  }

  if( !fileOpen ){
    
    fprintf( stderr,
             "File is closed, cannot read frame %d.\n"
             theFrame );
    exit(0);
  }

  // access the frame
  
  framePos = posArray[theFrame].myPos;
  fsetPos( inFile, framePos );

  for(j=0;j<3;j++){
    for(k=0;k<3;k++){
      Hmat[j][k] = posArray[theFrame].Hmat[j][k];
    }
  }
  
  fgets( readBuffer, sizeof( readBuffer ), inFile );
  if( feof( inFile ) ){
    fprintf( stderr,
             "File \"%s\" ended unexpectedly in Frame %d\n", 
             inName,
             theFrame );
    exit(0);
  }
  
  i = atoi( readBuffer );
  if( i != nAtoms ){
    fprintf( stderr,
             "Error in frame %d: frame atoms, %d, does not params atoms, %d\n",
             theFrame,
             i,
             nAtoms );
    exit(0);
  }

  // read and toss the comment line

  fgets( readBuffer, sizeof( readBuffer ), inFile );
  if( feof( inFile ) ){
    fprintf( stderr,
             "File \"%s\" ended unexpectedly in Frame %d\n", 
             inName,
             theFrame );
    exit(0);
  }

  // read the atoms

  for( i=0; i < nAtoms; i++){
    
    fgets( readBuffer, sizeof( readBuffer ), inFile );
    if( feof( inFile ) ){
      fprintf( stderr,
               "File \"%s\" ended unexpectedly in Frame %d at atom %d\n", 
               inName,
               theFrame,
               i );
      exit(0);
    }

    
    parseDumpLine( readBuffer, i, atoms );
  }
}

void parseDumpLine( char readLine[BUFFER_SIZE], int i, 
                    struct atomCoord* atoms ){
  char* foo;
  double q[4];

  // set the string tokenizer
  
  foo = strtok(readLine, " ,;\t");
  
  // check the atom ID

  switch( atoms[i].type ){
    
  case HEAD:
    if( strcmp(foo, "HEAD") ){
      fprintf( stderr,
               "Atom %s does not match master array type \"HEAD\".\n"
               foo );
      exit(0);
    }
    break;

  case CH:
    if( strcmp(foo, "CH") ){
      fprintf( stderr,
               "Atom %s does not match master array type \"CH\".\n"
               foo );
      exit(0);
    }
    break;

  case CH2:
    if( strcmp(foo, "CH2") ){
      fprintf( stderr,
               "Atom %s does not match master array type \"CH2\".\n"
               foo );
      exit(0);
    }
    break;

  case CH3:
    if( strcmp(foo, "CH3") ){
      fprintf( stderr,
               "Atom %s does not match master array type \"CH3\".\n"
               foo );
      exit(0);
    }
    break;

  case SSD:
    if( strcmp(foo, "SSD") ){
      fprintf( stderr,
               "Atom %s does not match master array type \"SSD\".\n"
               foo );
      exit(0);
    }
    break;
      
  default:
    fprintf(stderr,
            "Atom %d is an unrecognized enumeration\n",
            i );
    exit(0);
  }
  
  // get the positions
  
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    fprintf( stderr,
             "error in reading postition x from %s\n"
             "natoms  = %d, i = %d\n",
             inName, nAtoms, i );
    exit(0);
  }
  atoms[i].pos[0] = atof( foo );
   
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    fprintf( stderr,
             "error in reading postition y from %s\n"
             "natoms  = %d, i = %d\n",
             inName, nAtoms, i );
    exit(0);
  }
  atoms[i].pos[1] = atof( foo );
     
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    fprintf( stderr,
             "error in reading postition z from %s\n"
             "natoms  = %d, i = %d\n",
             inName, nAtoms, i );
    exit(0);
  }
  atoms[i].pos[2] = atof( foo );    


  // get the velocities

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    fprintf( stderr,
             "error in reading velocity x from %s\n"
             "natoms  = %d, i = %d\n",
             inName, nAtoms, i );
    exit(0);
  }
  atoms[i].vel[0] = atof( foo );
    
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    fprintf( stderr,
             "error in reading velocity y from %s\n"
             "natoms  = %d, i = %d\n",
             inName, nAtoms, i );
    exit(0);
  }
  atoms[i].vel[1] = atof( foo );
    
  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    fprintf( stderr,
             "error in reading velocity z from %s\n"
             "natoms  = %d, i = %d\n",
             inName, nAtoms, i );
    exit(0);
  }
  atoms[i].vel[2] = atof( foo );
    
    
  // get the quaternions
    
  if( (atoms[i].type == HEAD) || (atoms[i].type == SSD) ){
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf(painCave.errMsg,
              "error in reading quaternion 0 from %s\n"
              "natoms  = %d, i = %d\n",
              inName, nAtoms, i );
      exit(0);
    }
    q[0] = atof( foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      fprintf( stderr,
               "error in reading quaternion 1 from %s\n"
               "natoms  = %d, i = %d\n",
               inName, nAtoms, i );
      exit(0);
    }
    q[1] = atof( foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      fprintf( stderr,
               "error in reading quaternion 2 from %s\n"
               "natoms  = %d, i = %d\n",
               inName, nAtoms, i );
      exit(0);
    }
    q[2] = atof( foo );
      
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      fprintf( stderr,
               "error in reading quaternion 3 from %s\n"
               "natoms  = %d, i = %d\n",
               inName, nAtoms, i );
      exit(0);
    }
    q[3] = atof( foo );
  
    normalizeQ( q );

    setU( q, atoms[i].u );
  }
}


void setU( double the_q[4], double u[3] ){

  double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
  double A[3][3];
  double rb[3];

  // initialize the axis

  rb[0] = 0.0;
  rb[1] = 0.0;
  rb[2] = 1.0;

  // set the rotation matrix
  
  q0Sqr = the_q[0] * the_q[0];
  q1Sqr = the_q[1] * the_q[1];
  q2Sqr = the_q[2] * the_q[2];
  q3Sqr = the_q[3] * the_q[3];
  
  A[0][0] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
  A[0][1] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
  A[0][2] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
  
  A[1][0]  = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
  A[1][1] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
  A[1][2] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
  
  A[2][0] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
  A[2][1] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
  A[2][2] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;

  // perform the rotation
  
  for( i=0; i<3; i++ ){
    u[i] = 0.0;
    
    for( j=0; j<3; j++ ){
      u[i] += A[j][i] * rb[j];
    }
  }
}

void normalizeQ( double q[4] ){
  
  double qSqr, qL;
  int i;
  
  qSqr = 0.0;
  for(i=0;i<4;i++)
    qSqr += q[i] * q[i];
  
  for(i=0;i<4;i++)
    q[i] /= qSqr;
}
      
Revision:	1058
Committed:	Wed Feb 18 21:20:51 2004 UTC (20 years, 5 months ago) by mmeineke
Content type:	text/plain
File size:	12474 byte(s)
Log Message:	finished up the rough draft of tcProps with Scd corr
#	User	Rev	Content
1	mmeineke	1052	#define _FILE_OFFSET_BITS 64
2
3			#include <stdio.h>
4			#include <stdlib.h>
5			#include <sys/types.h>
6			#include <sys/stat.h>
7			#include <string.h>
8
9			#include "params.h"
10			#include "tcProps.h"
11			#include "readWrite.h"
12
13
14	mmeineke	1053	#define BUFFER_SIZE 2000
15	mmeineke	1052	int isScanned;
16	mmeineke	1056	int fileOpen;
17			int nFrames;
18			double *frameTimes;
19	mmeineke	1052	FILE* inFile;
20
21
22			struct linkedPos{
23
24			fpos_t *myPos;
25			double timeStamp;
26			double Hmat[3][3];
27			struct linkedPos* next;
28			};
29
30			struct staticPos{
31
32			fpos_t *myPos;
33			double Hmat[3][3];
34			};
35			struct staticPos* posArray;
36
37	mmeineke	1056	void parseDumpLine( char readBuffer[BUFFER_SIZE], int index,
38			struct atomCoord* atoms );
39	mmeineke	1058	void setU( double q[4], double u[3] );
40			void normalizeQ( double q[4] );
41	mmeineke	1052
42	mmeineke	1056	void openFile( char* inName ){
43
44			inFile = fopen(inName);
45			if(inFile ==NULL){
46			fprintf(stderr,
47			"Error opening file \"%s\"\n",
48			inName);
49			exit(0);
50			}
51
52			fileOpen = 1;
53			}
54
55	mmeineke	1052	void closeFile( void ){
56
57			fclose( inFile );
58	mmeineke	1056
59			fileOpen = 0;
60	mmeineke	1052	}
61
62	mmeineke	1056	int setFrames( void ){
63	mmeineke	1052
64			int i,j,k;
65			struct linkedPos* headPos;
66			struct linkedPos* currPos;
67			fpos_t *currPT;
68	mmeineke	1053	char readBuffer[BUFFER_SIZE];
69	mmeineke	1055	char trash[BUFFER_SIZE];
70	mmeineke	1053	char* foo;
71			int lineNum = 0;
72
73
74	mmeineke	1056	if( !fileOpen ){
75	mmeineke	1052	fprintf(stderr,
76	mmeineke	1056	"Attempt to scan the file without first opening the file.\n" );
77	mmeineke	1052	exit(0);
78			}
79
80	mmeineke	1056
81			nFrames = 0;
82	mmeineke	1053	headPos = (struct linkedPos*)malloc(sizeof(struct linkedPos));
83	mmeineke	1055	currPos = headPos;
84	mmeineke	1052	while( !feof( inFile ) ){
85
86	mmeineke	1053	currPT = (fpos_t *)malloc(sizeof(fpos_t));
87			fgetpos(inFile, currPT);
88
89			fgets( readBuffer, sizeof( readBuffer ), inFile );
90			lineNum++;
91			if( feof( inFile ) ){
92			fprintf( stderr,
93			"File \"%s\" ended unexpectedly at line %d\n",
94			inName,
95			lineNum );
96			exit(0);
97			}
98
99	mmeineke	1055	currPos->next = (struct linkedPos*)malloc(sizeof(struct linkedPos));
100			currPos = currPos->next;
101			currPos->myPos = currPT;
102			nFrames++;
103
104	mmeineke	1053	i = atoi(readBuffer);
105	mmeineke	1055
106			fgets( readBuffer, sizeof( readBuffer ), inFile );
107			lineNum++;
108			if( feof( inFile ) ){
109			fprintf( stderr,
110			"File \"%s\" ended unexpectedly at line %d\n",
111			inName,
112			lineNum );
113			exit(0);
114			}
115	mmeineke	1053
116
117	mmeineke	1055	// parse the comment line
118
119			foo = strtok(readBuffer, " ,;\t");
120
121			if(foo == NULL){
122			fprintf(stderr,
123			"error in reading time from %s at line %d\n",
124			inName, lineNum );
125			exit(0);
126			}
127			currPos->timeStamp = atof( foo );
128
129			// get the Hx vector
130
131			foo = strtok(NULL, " ,;\t");
132			if(foo == NULL){
133			fprintf(stderr,
134			"error in reading Hx[0] from %s at line %d\n",
135			inName, lineNum );
136			exit(0);
137			}
138			currPos->Hmat[0][0] = atof( foo );
139
140			foo = strtok(NULL, " ,;\t");
141			if(foo == NULL){
142			fprintf(stderr,
143			"error in reading Hx[1] from %s at line %d\n",
144			inName, lineNum );
145			exit(0);
146			}
147			currPos->Hmat[1][0] = atof( foo );
148
149			foo = strtok(NULL, " ,;\t");
150			if(foo == NULL){
151			fprintf(stderr,
152			"error in reading Hx[2] from %s at line %d\n",
153			inName, lineNum );
154			exit(0);
155			}
156			currPos->Hmat[2][0] = atof( foo );
157
158			// get the Hy vector
159
160			foo = strtok(NULL, " ,;\t");
161			if(foo == NULL){
162			fprintf(stderr,
163			"error in reading Hy[0] from %s at line %d\n",
164			inName, lineNum );
165			exit(0);
166			}
167			currPos->Hmat[0][1] = atof( foo );
168
169			foo = strtok(NULL, " ,;\t");
170			if(foo == NULL){
171			fprintf(stderr,
172			"error in reading Hy[1] from %s at line %d\n",
173			inName, lineNum );
174			exit(0);
175			}
176			currPos->Hmat[1][1] = atof( foo );
177
178			foo = strtok(NULL, " ,;\t");
179			if(foo == NULL){
180			fprintf(stderr,
181			"error in reading Hy[2] from %s at line %d\n",
182			inName, lineNum );
183			exit(0);
184			}
185			currPos->Hmat[2][1] = atof( foo );
186
187			// get the Hz vector
188
189			foo = strtok(NULL, " ,;\t");
190			if(foo == NULL){
191			fprintf(stderr,
192			"error in reading Hz[0] from %s at line %d\n",
193			inName, lineNum );
194			exit(0);
195			}
196			currPos->Hmat[0][2] = atof( foo );
197
198			foo = strtok(NULL, " ,;\t");
199			if(foo == NULL){
200			fprintf(stderr,
201			"error in reading Hz[1] from %s at line %d\n",
202			inName, lineNum );
203			exit(0);
204			}
205			currPos->Hmat[1][2] = atof( foo );
206
207			foo = strtok(NULL, " ,;\t");
208			if(foo == NULL){
209			fprintf(stderr,
210			"error in reading Hz[2] from %s at line %d\n",
211			inName, lineNum );
212			exit(0);
213			}
214			currPos->Hmat[2][2] = atof( foo );
215
216			// skip the atoms
217
218			for(j=0; j<i; j++){
219
220			fgets( readBuffer, sizeof( readBuffer ), inFile );
221			lineNum++;
222			if( feof( inFile ) ){
223			fprintf( stderr,
224			"File \"%s\" ended unexpectedly at line %d,"
225			" with atom %d\n",
226			inName,
227			lineNum,
228			j );
229			exit(0);
230			}
231			}
232	mmeineke	1053	}
233	mmeineke	1056
234			// allocate the static position array
235	mmeineke	1052
236	mmeineke	1056	posArray = (struct staticPos*)calloc(nFrames, sizeof(struct staticPos));
237			frameTimes = (double *)calloc(nFrames, sizeof(double) );
238
239			i=0;
240			currPos = headPos->next;
241			while( currPos != NULL ){
242
243			if(i>=nFrames){
244
245			fprintf( stderr,
246			"The number of frame pointers exceeded the number of frames.\n"
247			" OOPS!\n" );
248			exit(0);
249			}
250
251			frameTimes[i] = currPos->timeStamp;
252			posArray[i].myPos = currPos->myPos;
253			for(j=0;j<3;j++){
254			for(k=0;k<3;k++){
255			posArray[i].Hmat[j][k] = currPos->Hmat[j][k];
256			}
257			}
258
259			i++;
260			currPos = currPos->next;
261			}
262
263			// clean up the linked list
264
265			currPos = headPos;
266			while( currPos != NULL ){
267			headPos = currPos;
268			currPos = headPos->next;
269			free(headPos);
270			}
271
272
273	mmeineke	1052	isScanned = 1;
274			}
275	mmeineke	1056
276
277			void readFrame(int theFrame, struct atomCoord* atoms, double Hmat[3][3] ){
278
279			// list of 'a priori' constants
280
281			const int nLipAtoms = NL_ATOMS;
282			const int nBonds = NBONDS;
283			const int nLipids = NLIPIDS;
284			const int nSSD = NSSD;
285			const int nAtoms = nLipAtoms * nLipids + nSSD;
286
287			fpos_t *framePos;
288			char readBuffer[BUFFER_SIZE];
289
290			int i, j, k;
291
292
293			// quick checks
294
295			if (theFrame >= nFrames){
296
297			fprintf(stderr,
298			"frame %d, is out of range\n",
299			theFrame );
300			exit(0);
301			}
302
303			if( !fileOpen ){
304
305			fprintf( stderr,
306			"File is closed, cannot read frame %d.\n"
307			theFrame );
308			exit(0);
309			}
310
311			// access the frame
312
313			framePos = posArray[theFrame].myPos;
314			fsetPos( inFile, framePos );
315
316			for(j=0;j<3;j++){
317			for(k=0;k<3;k++){
318			Hmat[j][k] = posArray[theFrame].Hmat[j][k];
319			}
320			}
321
322			fgets( readBuffer, sizeof( readBuffer ), inFile );
323			if( feof( inFile ) ){
324			fprintf( stderr,
325			"File \"%s\" ended unexpectedly in Frame %d\n",
326			inName,
327			theFrame );
328			exit(0);
329			}
330
331			i = atoi( readBuffer );
332			if( i != nAtoms ){
333			fprintf( stderr,
334			"Error in frame %d: frame atoms, %d, does not params atoms, %d\n",
335			theFrame,
336			i,
337			nAtoms );
338			exit(0);
339			}
340
341			// read and toss the comment line
342
343			fgets( readBuffer, sizeof( readBuffer ), inFile );
344			if( feof( inFile ) ){
345			fprintf( stderr,
346			"File \"%s\" ended unexpectedly in Frame %d\n",
347			inName,
348			theFrame );
349			exit(0);
350			}
351
352			// read the atoms
353
354			for( i=0; i < nAtoms; i++){
355
356			fgets( readBuffer, sizeof( readBuffer ), inFile );
357			if( feof( inFile ) ){
358			fprintf( stderr,
359			"File \"%s\" ended unexpectedly in Frame %d at atom %d\n",
360			inName,
361			theFrame,
362			i );
363			exit(0);
364			}
365
366
367			parseDumpLine( readBuffer, i, atoms );
368			}
369			}
370
371	mmeineke	1058	void parseDumpLine( char readLine[BUFFER_SIZE], int i,
372	mmeineke	1056	struct atomCoord* atoms ){
373			char* foo;
374			double q[4];
375
376			// set the string tokenizer
377
378			foo = strtok(readLine, " ,;\t");
379
380			// check the atom ID
381
382	mmeineke	1058	switch( atoms[i].type ){
383	mmeineke	1056
384			case HEAD:
385			if( strcmp(foo, "HEAD") ){
386			fprintf( stderr,
387			"Atom %s does not match master array type \"HEAD\".\n"
388			foo );
389			exit(0);
390			}
391			break;
392
393			case CH:
394			if( strcmp(foo, "CH") ){
395			fprintf( stderr,
396			"Atom %s does not match master array type \"CH\".\n"
397			foo );
398			exit(0);
399			}
400			break;
401
402			case CH2:
403			if( strcmp(foo, "CH2") ){
404			fprintf( stderr,
405			"Atom %s does not match master array type \"CH2\".\n"
406			foo );
407			exit(0);
408			}
409			break;
410
411			case CH3:
412			if( strcmp(foo, "CH3") ){
413			fprintf( stderr,
414			"Atom %s does not match master array type \"CH3\".\n"
415			foo );
416			exit(0);
417			}
418			break;
419
420			case SSD:
421			if( strcmp(foo, "SSD") ){
422			fprintf( stderr,
423			"Atom %s does not match master array type \"SSD\".\n"
424			foo );
425			exit(0);
426			}
427			break;
428
429			default:
430			fprintf(stderr,
431			"Atom %d is an unrecognized enumeration\n",
432	mmeineke	1058	i );
433	mmeineke	1056	exit(0);
434			}
435
436			// get the positions
437
438			foo = strtok(NULL, " ,;\t");
439			if(foo == NULL){
440			fprintf( stderr,
441			"error in reading postition x from %s\n"
442	mmeineke	1058	"natoms = %d, i = %d\n",
443			inName, nAtoms, i );
444	mmeineke	1056	exit(0);
445			}
446			atoms[i].pos[0] = atof( foo );
447
448			foo = strtok(NULL, " ,;\t");
449			if(foo == NULL){
450			fprintf( stderr,
451			"error in reading postition y from %s\n"
452	mmeineke	1058	"natoms = %d, i = %d\n",
453			inName, nAtoms, i );
454	mmeineke	1056	exit(0);
455			}
456			atoms[i].pos[1] = atof( foo );
457
458			foo = strtok(NULL, " ,;\t");
459			if(foo == NULL){
460			fprintf( stderr,
461			"error in reading postition z from %s\n"
462	mmeineke	1058	"natoms = %d, i = %d\n",
463			inName, nAtoms, i );
464	mmeineke	1056	exit(0);
465			}
466			atoms[i].pos[2] = atof( foo );
467
468
469			// get the velocities
470
471			foo = strtok(NULL, " ,;\t");
472			if(foo == NULL){
473			fprintf( stderr,
474			"error in reading velocity x from %s\n"
475	mmeineke	1058	"natoms = %d, i = %d\n",
476			inName, nAtoms, i );
477	mmeineke	1056	exit(0);
478			}
479			atoms[i].vel[0] = atof( foo );
480
481			foo = strtok(NULL, " ,;\t");
482			if(foo == NULL){
483			fprintf( stderr,
484			"error in reading velocity y from %s\n"
485	mmeineke	1058	"natoms = %d, i = %d\n",
486			inName, nAtoms, i );
487	mmeineke	1056	exit(0);
488			}
489			atoms[i].vel[1] = atof( foo );
490
491			foo = strtok(NULL, " ,;\t");
492			if(foo == NULL){
493			fprintf( stderr,
494			"error in reading velocity z from %s\n"
495	mmeineke	1058	"natoms = %d, i = %d\n",
496			inName, nAtoms, i );
497	mmeineke	1056	exit(0);
498			}
499			atoms[i].vel[2] = atof( foo );
500
501
502			// get the quaternions
503
504			if( (atoms[i].type == HEAD) \|\| (atoms[i].type == SSD) ){
505
506			foo = strtok(NULL, " ,;\t");
507			if(foo == NULL){
508			sprintf(painCave.errMsg,
509			"error in reading quaternion 0 from %s\n"
510	mmeineke	1058	"natoms = %d, i = %d\n",
511			inName, nAtoms, i );
512	mmeineke	1056	exit(0);
513			}
514			q[0] = atof( foo );
515
516			foo = strtok(NULL, " ,;\t");
517			if(foo == NULL){
518			fprintf( stderr,
519			"error in reading quaternion 1 from %s\n"
520	mmeineke	1058	"natoms = %d, i = %d\n",
521			inName, nAtoms, i );
522	mmeineke	1056	exit(0);
523			}
524			q[1] = atof( foo );
525
526			foo = strtok(NULL, " ,;\t");
527			if(foo == NULL){
528			fprintf( stderr,
529			"error in reading quaternion 2 from %s\n"
530	mmeineke	1058	"natoms = %d, i = %d\n",
531			inName, nAtoms, i );
532	mmeineke	1056	exit(0);
533			}
534			q[2] = atof( foo );
535
536			foo = strtok(NULL, " ,;\t");
537			if(foo == NULL){
538			fprintf( stderr,
539			"error in reading quaternion 3 from %s\n"
540	mmeineke	1058	"natoms = %d, i = %d\n",
541			inName, nAtoms, i );
542	mmeineke	1056	exit(0);
543			}
544			q[3] = atof( foo );
545	mmeineke	1058
546			normalizeQ( q );
547
548			setU( q, atoms[i].u );
549	mmeineke	1056	}
550	mmeineke	1058	}
551	mmeineke	1056
552
553	mmeineke	1058	void setU( double the_q[4], double u[3] ){
554	mmeineke	1056
555	mmeineke	1058	double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
556			double A[3][3];
557			double rb[3];
558	mmeineke	1056
559	mmeineke	1058	// initialize the axis
560	mmeineke	1056
561	mmeineke	1058	rb[0] = 0.0;
562			rb[1] = 0.0;
563			rb[2] = 1.0;
564
565			// set the rotation matrix
566
567			q0Sqr = the_q[0] * the_q[0];
568			q1Sqr = the_q[1] * the_q[1];
569			q2Sqr = the_q[2] * the_q[2];
570			q3Sqr = the_q[3] * the_q[3];
571
572			A[0][0] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
573			A[0][1] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
574			A[0][2] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
575
576			A[1][0] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
577			A[1][1] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
578			A[1][2] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
579
580			A[2][0] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
581			A[2][1] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
582			A[2][2] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
583
584			// perform the rotation
585
586			for( i=0; i<3; i++ ){
587			u[i] = 0.0;
588
589			for( j=0; j<3; j++ ){
590			u[i] += A[j][i] * rb[j];
591			}
592			}
593	mmeineke	1056	}
594	mmeineke	1058
595			void normalizeQ( double q[4] ){
596
597			double qSqr, qL;
598			int i;
599
600			qSqr = 0.0;
601			for(i=0;i<4;i++)
602			qSqr += q[i] * q[i];
603
604			for(i=0;i<4;i++)
605			q[i] /= qSqr;
606			}
607