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;
char* inName;

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( void ){
  
  inFile = fopen(inName, "r");
  if(inFile ==NULL){
    fprintf(stderr,
            "Error opening file \"%s\"\n",
            inName);
    exit(0);
  }

  fileOpen = 1;
}  

void closeFile( void ){

  fclose( inFile );

  fileOpen = 0;
}

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

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

  while( !feof( inFile ) ){

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

    currPT = (fpos_t *)malloc(sizeof(fpos_t));
    fgetpos(inFile, currPT);
    
    fgets( readBuffer, sizeof( readBuffer ), inFile );
    lineNum++;
  }
  
  // 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 ){
  
  const int nLipAtoms = NL_ATOMS;
  const int nBonds    = NBONDS;
  const int nLipids   = NLIPIDS;
  const int nSSD      = NSSD;
  const int nAtoms    = nLipAtoms * nLipids + nSSD;

  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){
      fprintf( stderr,
               "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];
  int i,j,k;

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