ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/tcProps/readWrite.c
Revision: 1085
Committed: Fri Mar 5 03:10:29 2004 UTC (20 years, 3 months ago) by mmeineke
Content type: text/plain
File size: 13671 byte(s)
Log Message:
who knows?, but it works for the most part

File Contents

# Content
1 #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 #define BUFFER_SIZE 2000
15 int isScanned;
16 int fileOpen;
17 int nFrames;
18 double *frameTimes;
19 FILE* inFile;
20 char* inName;
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 void parseDumpLine( char readBuffer[BUFFER_SIZE], int index,
38 struct atomCoord* atoms );
39 void setU( double q[4], double u[3] );
40 void normalizeQ( double q[4] );
41
42 void openFile( void ){
43
44 inFile = fopen(inName, "r");
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 void closeFile( void ){
56
57 fclose( inFile );
58
59 fileOpen = 0;
60 }
61
62 void setFrames( void ){
63
64 int i,j,k;
65 struct linkedPos* headPos;
66 struct linkedPos* currPos;
67 fpos_t *currPT;
68 char readBuffer[BUFFER_SIZE];
69 char trash[BUFFER_SIZE];
70 char* foo;
71 int lineNum = 0;
72
73
74 if( !fileOpen ){
75 fprintf(stderr,
76 "Attempt to scan the file without first opening the file.\n" );
77 exit(0);
78 }
79
80
81 nFrames = 0;
82 headPos = (struct linkedPos*)malloc(sizeof(struct linkedPos));
83 currPos = headPos;
84
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 while( !feof( inFile ) ){
99
100 currPos->next = (struct linkedPos*)malloc(sizeof(struct linkedPos));
101 currPos = currPos->next;
102 currPos->myPos = currPT;
103 nFrames++;
104
105 i = atoi(readBuffer);
106
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
117
118 // 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
234 currPT = (fpos_t *)malloc(sizeof(fpos_t));
235 fgetpos(inFile, currPT);
236
237 fgets( readBuffer, sizeof( readBuffer ), inFile );
238 lineNum++;
239 }
240
241 // allocate the static position array
242
243 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 isScanned = 1;
281 }
282
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 "File is closed, cannot read frame %d.\n",
314 theFrame );
315 exit(0);
316 }
317
318 // access the frame
319
320 framePos = posArray[theFrame].myPos;
321 fsetpos( inFile, framePos );
322
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 void parseDumpLine( char readLine[BUFFER_SIZE], int i,
379 struct atomCoord* atoms ){
380
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 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 switch( atoms[i].type ){
397
398 case HEAD:
399 if( strcmp(foo, "HEAD") ){
400 fprintf( stderr,
401 "Atom %s does not match master array type \"HEAD\".\n",
402 foo );
403 exit(0);
404 }
405 break;
406
407 case CH:
408 if( strcmp(foo, "CH") ){
409 fprintf( stderr,
410 "Atom %s does not match master array type \"CH\".\n",
411 foo );
412 exit(0);
413 }
414 break;
415
416 case CH2:
417 if( strcmp(foo, "CH2") ){
418 fprintf( stderr,
419 "Atom %s does not match master array type \"CH2\".\n",
420 foo );
421 exit(0);
422 }
423 break;
424
425 case CH3:
426 if( strcmp(foo, "CH3") ){
427 fprintf( stderr,
428 "Atom %s does not match master array type \"CH3\".\n",
429 foo );
430 exit(0);
431 }
432 break;
433
434 case SSD:
435 if( strcmp(foo, "SSD") ){
436 fprintf( stderr,
437 "Atom %s does not match master array type \"SSD\".\n",
438 foo );
439 exit(0);
440 }
441 break;
442
443 default:
444 fprintf(stderr,
445 "Atom %d is an unrecognized enumeration\n",
446 i );
447 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 "natoms = %d, i = %d\n",
457 inName, nAtoms, i );
458 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 "natoms = %d, i = %d\n",
467 inName, nAtoms, i );
468 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 "natoms = %d, i = %d\n",
477 inName, nAtoms, i );
478 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 "natoms = %d, i = %d\n",
490 inName, nAtoms, i );
491 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 "natoms = %d, i = %d\n",
500 inName, nAtoms, i );
501 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 "natoms = %d, i = %d\n",
510 inName, nAtoms, i );
511 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 fprintf( stderr,
523 "error in reading quaternion 0 from %s\n"
524 "natoms = %d, i = %d\n",
525 inName, nAtoms, i );
526 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 "natoms = %d, i = %d\n",
535 inName, nAtoms, i );
536 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 "natoms = %d, i = %d\n",
545 inName, nAtoms, i );
546 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 "natoms = %d, i = %d\n",
555 inName, nAtoms, i );
556 exit(0);
557 }
558 q[3] = atof( foo );
559
560 normalizeQ( q );
561
562 setU( q, atoms[i].u );
563 }
564 }
565
566
567 void setU( double the_q[4], double u[3] ){
568
569 double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
570 double A[3][3];
571 double rb[3];
572 int i,j,k;
573
574 // initialize the axis
575
576 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 }
609
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
623 double scanSmallestZ(int startFrame){
624 double smallest;
625 int i;
626
627 smallest = posArray[startFrame].Hmat[2][2];
628 for(i=startFrame;i<nFrames;i++)
629 smallest =
630 (smallest<posArray[i].Hmat[2][2])? smallest : posArray[i].Hmat[2][2];
631
632 return smallest;
633 }
634
635 double calcAverageXY( double startTime ){
636
637 int i, count, startFrame, startFound;
638 double avg;
639
640 startFound = 0;
641 startFrame = -1;
642 while( !startFound ){
643
644 startFrame++;
645
646 if(startFrame >= nFrames){
647
648 fprintf( stderr,
649 "Start Time, %G, was not found in the dump file.\n",
650 startTime );
651 exit(0);
652 }
653
654 if(startTime <= frameTimes[startFrame])
655 startFound = 1;
656
657
658 }
659
660 avg = 0.0;
661 count = 0;
662 for(i=startFrame;i<nFrames;i++){
663 avg += posArray[i].Hmat[0][0] * posArray[i].Hmat[1][1];
664 count++;
665 }
666 avg /= (double)count;
667
668 return avg;
669 }