ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/tcProps/readWrite.c
Revision: 1058
Committed: Wed Feb 18 21:20:51 2004 UTC (20 years, 4 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

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
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( 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 void closeFile( void ){
56
57 fclose( inFile );
58
59 fileOpen = 0;
60 }
61
62 int 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 while( !feof( inFile ) ){
85
86 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 currPos->next = (struct linkedPos*)malloc(sizeof(struct linkedPos));
100 currPos = currPos->next;
101 currPos->myPos = currPT;
102 nFrames++;
103
104 i = atoi(readBuffer);
105
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
116
117 // 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 }
233
234 // allocate the static position array
235
236 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 isScanned = 1;
274 }
275
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 void parseDumpLine( char readLine[BUFFER_SIZE], int i,
372 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 switch( atoms[i].type ){
383
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 i );
433 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 "natoms = %d, i = %d\n",
443 inName, nAtoms, i );
444 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 "natoms = %d, i = %d\n",
453 inName, nAtoms, i );
454 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 "natoms = %d, i = %d\n",
463 inName, nAtoms, i );
464 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 "natoms = %d, i = %d\n",
476 inName, nAtoms, i );
477 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 "natoms = %d, i = %d\n",
486 inName, nAtoms, i );
487 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 "natoms = %d, i = %d\n",
496 inName, nAtoms, i );
497 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 "natoms = %d, i = %d\n",
511 inName, nAtoms, i );
512 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 "natoms = %d, i = %d\n",
521 inName, nAtoms, i );
522 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 "natoms = %d, i = %d\n",
531 inName, nAtoms, i );
532 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 "natoms = %d, i = %d\n",
541 inName, nAtoms, i );
542 exit(0);
543 }
544 q[3] = atof( foo );
545
546 normalizeQ( q );
547
548 setU( q, atoms[i].u );
549 }
550 }
551
552
553 void setU( double the_q[4], double u[3] ){
554
555 double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
556 double A[3][3];
557 double rb[3];
558
559 // initialize the axis
560
561 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 }
594
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