1 |
#include <stdio.h> |
2 |
#include <stdlib.h> |
3 |
#include <string.h> |
4 |
#include <math.h> |
5 |
#include <unistd.h> |
6 |
#include <sys/types.h> |
7 |
#include <sys/stat.h> |
8 |
|
9 |
struct coords{ |
10 |
double x,y,z; // cartesian coords |
11 |
double q[4]; // the quanternions |
12 |
char name[30]; |
13 |
}; |
14 |
|
15 |
struct linked_xyz{ |
16 |
struct coords *r; |
17 |
double time; |
18 |
double boxX, boxY, boxZ; |
19 |
struct linked_xyz *next; |
20 |
}; |
21 |
|
22 |
char *program_name; /*the name of the program */ |
23 |
int printWater = 0; |
24 |
int printDipole = 0; |
25 |
|
26 |
void usage(void); |
27 |
|
28 |
void rotWrite( FILE* out_file, struct coords* theSSD ); |
29 |
|
30 |
void setRot( double A[3][3], double q[4] ); |
31 |
|
32 |
void rotMe( double A[3][3], double r[3] ); |
33 |
|
34 |
void map( double *x, double *y, double *z, double centerX, double centerY, |
35 |
double centerZ, double boxX, double boxY, double boxZ ); |
36 |
|
37 |
int main(argc, argv) |
38 |
int argc; |
39 |
char *argv[]; |
40 |
{ |
41 |
|
42 |
|
43 |
struct coords *out_coords; |
44 |
|
45 |
int i, j, k, l; /* loop counters */ |
46 |
mode_t dir_mode = S_IRWXU; |
47 |
|
48 |
int lineCount = 0; //the line number |
49 |
int newN; // the new number of atoms |
50 |
int nSSD=0; // the number of SSD per frame |
51 |
|
52 |
unsigned int nAtoms; /*the number of atoms in each time step */ |
53 |
char read_buffer[1000]; /*the line buffer for reading */ |
54 |
char *eof_test; /*ptr to see when we reach the end of the file */ |
55 |
char *foo; /*the pointer to the current string token */ |
56 |
FILE *in_file; /* the input file */ |
57 |
FILE *out_file; /*the output file */ |
58 |
char *out_prefix = NULL; /*the prefix of the output file */ |
59 |
char out_name[500]; /*the output name */ |
60 |
int have_outName =0; |
61 |
char in_name[500]; // the input file name |
62 |
char temp_name[500]; |
63 |
char *root_name = NULL; /*the root name of the input file */ |
64 |
unsigned int n_out = 0; /*keeps track of which output file is being written*/ |
65 |
|
66 |
|
67 |
int skipFrames = 1; |
68 |
|
69 |
struct linked_xyz *current_frame; |
70 |
struct linked_xyz *temp_frame; |
71 |
|
72 |
int nframes =0; |
73 |
|
74 |
int wrap = 0; |
75 |
|
76 |
double cX = 0.0; |
77 |
double cY = 0.0; |
78 |
double cZ = 0.0; |
79 |
|
80 |
double mcX, mcY, mcZ; |
81 |
double newCX, newCY, newCZ; |
82 |
double dx, dy, dz; |
83 |
int mcount; |
84 |
|
85 |
int repeatX = 0; |
86 |
int repeatY = 0; |
87 |
int repeatZ = 0; |
88 |
|
89 |
int nRepeatX = 0; |
90 |
int nRepeatY = 0; |
91 |
int nRepeatZ = 0; |
92 |
|
93 |
struct coords* rCopy; |
94 |
|
95 |
int done; |
96 |
char current_flag; |
97 |
|
98 |
program_name = argv[0]; /*save the program name in case we need it*/ |
99 |
|
100 |
for( i = 1; i < argc; i++){ |
101 |
|
102 |
if(argv[i][0] =='-'){ |
103 |
|
104 |
if(argv[i][1] == '-' ){ |
105 |
|
106 |
// parse long word options |
107 |
|
108 |
if( !strcmp( argv[i], "--repeatX" ) ){ |
109 |
repeatX = 1; |
110 |
i++; |
111 |
nRepeatX = atoi( argv[i] ); |
112 |
} |
113 |
|
114 |
else if( !strcmp( argv[i], "--repeatY" ) ){ |
115 |
repeatY = 1; |
116 |
i++; |
117 |
nRepeatY = atoi( argv[i] ); |
118 |
} |
119 |
|
120 |
else if( !strcmp( argv[i], "--repeatZ" ) ){ |
121 |
repeatZ = 1; |
122 |
i++; |
123 |
nRepeatZ = atoi( argv[i] ); |
124 |
} |
125 |
|
126 |
else{ |
127 |
fprintf( stderr, |
128 |
"Invalid option \"%s\"\n", argv[i] ); |
129 |
usage(); |
130 |
} |
131 |
} |
132 |
|
133 |
else{ |
134 |
|
135 |
// parse single character options |
136 |
|
137 |
done =0; |
138 |
j = 1; |
139 |
current_flag = argv[i][j]; |
140 |
while( (current_flag != '\0') && (!done) ){ |
141 |
|
142 |
switch(current_flag){ |
143 |
|
144 |
case 'o': |
145 |
// -o <prefix> => the output |
146 |
|
147 |
i++; |
148 |
strcpy( out_name, argv[i] ); |
149 |
have_outName = 1; |
150 |
done = 1; |
151 |
break; |
152 |
|
153 |
case 'n': |
154 |
// -n <#> => print every <#> frames |
155 |
|
156 |
i++; |
157 |
skipFrames = atoi(argv[i]); |
158 |
done = 1; |
159 |
break; |
160 |
|
161 |
case 'h': |
162 |
// -h => give the usage |
163 |
|
164 |
usage(); |
165 |
break; |
166 |
|
167 |
case 'd': |
168 |
// print the dipoles |
169 |
|
170 |
printDipole = 1; |
171 |
break; |
172 |
|
173 |
case 'w': |
174 |
// print the waters |
175 |
|
176 |
printWater = 1; |
177 |
break; |
178 |
|
179 |
case 'm': |
180 |
// map to a periodic box |
181 |
|
182 |
wrap = 1; |
183 |
break; |
184 |
|
185 |
default: |
186 |
|
187 |
(void)fprintf(stderr, "Bad option \"-%s\"\n", current_flag); |
188 |
usage(); |
189 |
} |
190 |
j++; |
191 |
current_flag = argv[i][j]; |
192 |
} |
193 |
} |
194 |
} |
195 |
|
196 |
else{ |
197 |
|
198 |
if( root_name != NULL ){ |
199 |
fprintf( stderr, |
200 |
"Error at \"%s\", program does not currently support\n" |
201 |
"more than one input file.\n" |
202 |
"\n", |
203 |
argv[i]); |
204 |
usage(); |
205 |
} |
206 |
|
207 |
root_name = argv[i]; |
208 |
} |
209 |
} |
210 |
|
211 |
if(root_name == NULL){ |
212 |
usage(); |
213 |
} |
214 |
|
215 |
sprintf( in_name, "%s", root_name ); |
216 |
if( !have_outName ) sprintf( out_name, "%s.xyz", root_name ); |
217 |
|
218 |
in_file = fopen(in_name, "r"); |
219 |
if(in_file == NULL){ |
220 |
printf("Cannot open file \"%s\" for reading.\n", in_name); |
221 |
exit(8); |
222 |
} |
223 |
|
224 |
out_file = fopen( out_name, "w" ); |
225 |
if( out_file == NULL ){ |
226 |
printf("Cannot open file \"%s\" for writing.\n", out_name); |
227 |
exit(8); |
228 |
} |
229 |
|
230 |
// start reading the first frame |
231 |
|
232 |
eof_test = fgets(read_buffer, sizeof(read_buffer), in_file); |
233 |
lineCount++; |
234 |
|
235 |
current_frame = (struct linked_xyz *)malloc(sizeof(struct linked_xyz)); |
236 |
current_frame->next = NULL; |
237 |
|
238 |
|
239 |
while(eof_test != NULL){ |
240 |
|
241 |
(void)sscanf(read_buffer, "%d", &nAtoms); |
242 |
current_frame->r = |
243 |
(struct coords *)calloc(nAtoms, sizeof(struct coords)); |
244 |
|
245 |
// read and the comment line and grab the time and box dimensions |
246 |
|
247 |
eof_test = fgets(read_buffer, sizeof(read_buffer), in_file); |
248 |
lineCount++; |
249 |
if(eof_test == NULL){ |
250 |
printf("error in reading file at line: %d\n", lineCount); |
251 |
exit(8); |
252 |
} |
253 |
|
254 |
foo = strtok( read_buffer, " ,;\t" ); |
255 |
sscanf( read_buffer, "%lf", ¤t_frame->time ); |
256 |
|
257 |
foo = strtok(NULL, " ,;\t"); |
258 |
if(foo == NULL){ |
259 |
printf("error in reading file at line: %d\n", lineCount); |
260 |
exit(8); |
261 |
} |
262 |
(void)sscanf(foo, "%lf",¤t_frame->boxX); |
263 |
|
264 |
foo = strtok(NULL, " ,;\t"); |
265 |
if(foo == NULL){ |
266 |
printf("error in reading file at line: %d\n", lineCount); |
267 |
exit(8); |
268 |
} |
269 |
(void)sscanf(foo, "%lf",¤t_frame->boxY); |
270 |
|
271 |
foo = strtok(NULL, " ,;\t"); |
272 |
if(foo == NULL){ |
273 |
printf("error in reading file at line: %d\n", lineCount); |
274 |
exit(8); |
275 |
} |
276 |
(void)sscanf(foo, "%lf",¤t_frame->boxZ); |
277 |
|
278 |
for( i=0; i < nAtoms; i++){ |
279 |
|
280 |
eof_test = fgets(read_buffer, sizeof(read_buffer), in_file); |
281 |
lineCount++; |
282 |
if(eof_test == NULL){ |
283 |
printf("error in reading file at line: %d\n", lineCount); |
284 |
exit(8); |
285 |
} |
286 |
|
287 |
foo = strtok(read_buffer, " ,;\t"); |
288 |
if ( !strcmp( "SSD", foo ) ) nSSD++; |
289 |
(void)strcpy(current_frame->r[i].name, foo); //copy the atom name |
290 |
|
291 |
// next we grab the positions |
292 |
|
293 |
foo = strtok(NULL, " ,;\t"); |
294 |
if(foo == NULL){ |
295 |
printf("error in reading postition x from %s\n" |
296 |
"natoms = %d, line = %d\n", |
297 |
in_name, nAtoms, lineCount ); |
298 |
exit(8); |
299 |
} |
300 |
(void)sscanf( foo, "%lf", ¤t_frame->r[i].x ); |
301 |
|
302 |
foo = strtok(NULL, " ,;\t"); |
303 |
if(foo == NULL){ |
304 |
printf("error in reading postition y from %s\n" |
305 |
"natoms = %d, line = %d\n", |
306 |
in_name, nAtoms, lineCount ); |
307 |
exit(8); |
308 |
} |
309 |
(void)sscanf( foo, "%lf", ¤t_frame->r[i].y ); |
310 |
|
311 |
foo = strtok(NULL, " ,;\t"); |
312 |
if(foo == NULL){ |
313 |
printf("error in reading postition z from %s\n" |
314 |
"natoms = %d, line = %d\n", |
315 |
in_name, nAtoms, lineCount ); |
316 |
exit(8); |
317 |
} |
318 |
(void)sscanf( foo, "%lf", ¤t_frame->r[i].z ); |
319 |
|
320 |
// get the velocities |
321 |
|
322 |
foo = strtok(NULL, " ,;\t"); |
323 |
if(foo == NULL){ |
324 |
printf("error in reading velocity x from %s\n" |
325 |
"natoms = %d, line = %d\n", |
326 |
in_name, nAtoms, lineCount ); |
327 |
exit(8); |
328 |
} |
329 |
|
330 |
|
331 |
foo = strtok(NULL, " ,;\t"); |
332 |
if(foo == NULL){ |
333 |
printf("error in reading velocity y from %s\n" |
334 |
"natoms = %d, line = %d\n", |
335 |
in_name, nAtoms, lineCount ); |
336 |
exit(8); |
337 |
} |
338 |
|
339 |
|
340 |
foo = strtok(NULL, " ,;\t"); |
341 |
if(foo == NULL){ |
342 |
printf("error in reading velocity z from %s\n" |
343 |
"natoms = %d, line = %d\n", |
344 |
in_name, nAtoms, lineCount ); |
345 |
exit(8); |
346 |
} |
347 |
|
348 |
// get the quaternions |
349 |
|
350 |
foo = strtok(NULL, " ,;\t"); |
351 |
if(foo == NULL){ |
352 |
printf("error in reading quaternion 0 from %s\n" |
353 |
"natoms = %d, line = %d\n", |
354 |
in_name, nAtoms, lineCount ); |
355 |
exit(8); |
356 |
} |
357 |
(void)sscanf( foo, "%lf", ¤t_frame->r[i].q[0] ); |
358 |
|
359 |
foo = strtok(NULL, " ,;\t"); |
360 |
if(foo == NULL){ |
361 |
printf("error in reading quaternion 1 from %s\n" |
362 |
"natoms = %d, line = %d\n", |
363 |
in_name, nAtoms, lineCount ); |
364 |
exit(8); |
365 |
} |
366 |
(void)sscanf( foo, "%lf", ¤t_frame->r[i].q[1] ); |
367 |
|
368 |
foo = strtok(NULL, " ,;\t"); |
369 |
if(foo == NULL){ |
370 |
printf("error in reading quaternion 2 from %s\n" |
371 |
"natoms = %d, line = %d\n", |
372 |
in_name, nAtoms, lineCount ); |
373 |
exit(8); |
374 |
} |
375 |
(void)sscanf( foo, "%lf", ¤t_frame->r[i].q[2] ); |
376 |
|
377 |
foo = strtok(NULL, " ,;\t"); |
378 |
if(foo == NULL){ |
379 |
printf("error in reading quaternion 3 from %s\n" |
380 |
"natoms = %d, line = %d\n", |
381 |
in_name, nAtoms, lineCount ); |
382 |
exit(8); |
383 |
} |
384 |
(void)sscanf( foo, "%lf", ¤t_frame->r[i].q[3] ); |
385 |
|
386 |
// get the angular velocities |
387 |
|
388 |
foo = strtok(NULL, " ,;\t"); |
389 |
if(foo == NULL){ |
390 |
printf("error in reading angular momentum jx from %s\n" |
391 |
"natoms = %d, line = %d\n", |
392 |
in_name, nAtoms, lineCount ); |
393 |
exit(8); |
394 |
} |
395 |
|
396 |
foo = strtok(NULL, " ,;\t"); |
397 |
if(foo == NULL){ |
398 |
printf("error in reading angular momentum jy from %s\n" |
399 |
"natoms = %d, line = %d\n", |
400 |
in_name, nAtoms, lineCount ); |
401 |
exit(8); |
402 |
} |
403 |
|
404 |
foo = strtok(NULL, " ,;\t"); |
405 |
if(foo == NULL){ |
406 |
printf("error in reading angular momentum jz from %s\n" |
407 |
"natoms = %d, line = %d\n", |
408 |
in_name, nAtoms, lineCount ); |
409 |
exit(8); |
410 |
} |
411 |
} |
412 |
|
413 |
// if necassary, wrap into the periodic image |
414 |
|
415 |
if(wrap){ |
416 |
|
417 |
for( i=0; i<nAtoms; i++ ){ |
418 |
|
419 |
/* |
420 |
if( !strcmp( current_frame->r[i].name, "HEAD" ) ){ |
421 |
|
422 |
// find the center of the lipid |
423 |
|
424 |
mcX = 0.0; |
425 |
mcY = 0.0; |
426 |
mcZ = 0.0; |
427 |
mcount = 0; |
428 |
|
429 |
mcX += current_frame->r[i].x; |
430 |
mcY += current_frame->r[i].y; |
431 |
mcZ += current_frame->r[i].z; |
432 |
mcount++; |
433 |
i++; |
434 |
|
435 |
while( strcmp( current_frame->r[i].name, "HEAD" ) && |
436 |
strcmp( current_frame->r[i].name, "SSD" ) ){ |
437 |
|
438 |
mcX += current_frame->r[i].x; |
439 |
mcY += current_frame->r[i].y; |
440 |
mcZ += current_frame->r[i].z; |
441 |
mcount++; |
442 |
i++; |
443 |
} |
444 |
|
445 |
mcX /= mcount; |
446 |
mcY /= mcount; |
447 |
mcZ /= mcount; |
448 |
|
449 |
newCX = mcX; |
450 |
newCY = mcY; |
451 |
newCZ = mcZ; |
452 |
|
453 |
map( &newCX, &newCY, &newCZ, |
454 |
cX, cY, cZ, |
455 |
current_frame->boxX, |
456 |
current_frame->boxY, |
457 |
current_frame->boxZ ); |
458 |
|
459 |
dx = mcX - newCX; |
460 |
dy = mcY - newCY; |
461 |
dz = mcZ - newCZ; |
462 |
|
463 |
for(j=(i-mcount); j<mcount; j++ ){ |
464 |
|
465 |
current_frame->r[j].x -= dx; |
466 |
current_frame->r[j].y -= dy; |
467 |
current_frame->r[j].z -= dz; |
468 |
} |
469 |
|
470 |
i--; // so we don't skip the next atom |
471 |
} |
472 |
*/ |
473 |
// else |
474 |
map( ¤t_frame->r[i].x, |
475 |
¤t_frame->r[i].y, |
476 |
¤t_frame->r[i].z, |
477 |
cX, cY, cZ, |
478 |
current_frame->boxX, |
479 |
current_frame->boxY, |
480 |
current_frame->boxZ ); |
481 |
|
482 |
} |
483 |
} |
484 |
|
485 |
|
486 |
|
487 |
|
488 |
nframes++; |
489 |
|
490 |
// write out here |
491 |
|
492 |
if(printWater) newN = nAtoms + ( nSSD * 3 ); |
493 |
else newN = nAtoms - nSSD; |
494 |
|
495 |
newN *= (nRepeatX+1); |
496 |
newN *= (nRepeatY+1); |
497 |
newN *= (nRepeatZ+1); |
498 |
|
499 |
if( !(nframes%skipFrames) ){ |
500 |
fprintf( out_file, |
501 |
"%d\n" |
502 |
"%lf\t%lf\t%lf\t%lf\n", |
503 |
newN, current_frame->time, |
504 |
current_frame->boxX, current_frame->boxY, |
505 |
current_frame->boxZ ); |
506 |
|
507 |
rCopy = (struct coords*) |
508 |
calloc( nAtoms, sizeof( struct coords )); |
509 |
|
510 |
for(i=0; i<nAtoms; i++){ |
511 |
rCopy[i].q[0] = current_frame->r[i].q[0]; |
512 |
rCopy[i].q[1] = current_frame->r[i].q[1]; |
513 |
rCopy[i].q[2] = current_frame->r[i].q[2]; |
514 |
rCopy[i].q[3] = current_frame->r[i].q[3]; |
515 |
|
516 |
strcpy(rCopy[i].name, current_frame->r[i].name); |
517 |
} |
518 |
|
519 |
for(j=0; j<(nRepeatX+1); j++){ |
520 |
for(k=0; k<(nRepeatY+1); k++){ |
521 |
for(l=0; l<(nRepeatZ+1); l++){ |
522 |
|
523 |
for(i=0; i<nAtoms; i++){ |
524 |
rCopy[i].x = current_frame->r[i].x + (j * current_frame->boxX); |
525 |
rCopy[i].y = current_frame->r[i].y + (k * current_frame->boxY); |
526 |
rCopy[i].z = current_frame->r[i].z + (l * current_frame->boxZ); |
527 |
|
528 |
if( !strcmp( "SSD", rCopy[i].name ) ){ |
529 |
|
530 |
rotWrite( out_file, &rCopy[i] ); |
531 |
} |
532 |
|
533 |
else if( !strcmp( "HEAD", rCopy[i].name ) ){ |
534 |
|
535 |
rotWrite( out_file, &rCopy[i] ); |
536 |
} |
537 |
else{ |
538 |
|
539 |
fprintf( out_file, |
540 |
"%s\t%lf\t%lf\t%lf\n", |
541 |
rCopy[i].name, |
542 |
rCopy[i].x, |
543 |
rCopy[i].y, |
544 |
rCopy[i].z ); |
545 |
} |
546 |
} |
547 |
} |
548 |
} |
549 |
} |
550 |
|
551 |
free( rCopy ); |
552 |
} |
553 |
|
554 |
/*free up memory */ |
555 |
|
556 |
temp_frame = current_frame->next; |
557 |
current_frame->next = NULL; |
558 |
|
559 |
if(temp_frame != NULL){ |
560 |
|
561 |
free(temp_frame->r); |
562 |
free(temp_frame); |
563 |
} |
564 |
|
565 |
/* make a new frame */ |
566 |
|
567 |
temp_frame = (struct linked_xyz *)malloc(sizeof(struct linked_xyz)); |
568 |
temp_frame->next = current_frame; |
569 |
current_frame = temp_frame; |
570 |
|
571 |
nSSD = 0; |
572 |
|
573 |
eof_test = fgets(read_buffer, sizeof(read_buffer), in_file); |
574 |
lineCount++; |
575 |
} |
576 |
|
577 |
(void)fclose(in_file); |
578 |
|
579 |
return 0; |
580 |
|
581 |
} |
582 |
|
583 |
|
584 |
|
585 |
void rotWrite( FILE* out_file, struct coords* theDipole ){ |
586 |
|
587 |
double u[3]; |
588 |
double h1[3]; |
589 |
double h2[3]; |
590 |
double ox[3]; |
591 |
|
592 |
double A[3][3]; |
593 |
|
594 |
|
595 |
|
596 |
u[0] = 0.0; u[1] = 0.0; u[2] = 1.0; |
597 |
|
598 |
h1[0] = 0.0; h1[1] = -0.75695; h1[2] = 0.5206; |
599 |
|
600 |
h2[0] = 0.0; h2[1] = 0.75695; h2[2] = 0.5206; |
601 |
|
602 |
ox[0] = 0.0; ox[1] = 0.0; ox[2] = -0.0654; |
603 |
|
604 |
|
605 |
setRot( A, theDipole->q ); |
606 |
rotMe( A, u ); |
607 |
|
608 |
if( !strcmp( "SSD", theDipole->name )){ |
609 |
|
610 |
rotMe( A, h1 ); |
611 |
rotMe( A, h2 ); |
612 |
rotMe( A, ox ); |
613 |
|
614 |
if( printWater ){ |
615 |
|
616 |
if(printDipole) { |
617 |
|
618 |
fprintf( out_file, |
619 |
"X\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", |
620 |
theDipole->x, |
621 |
theDipole->y, |
622 |
theDipole->z, |
623 |
u[0], u[1], u[2] ); |
624 |
} |
625 |
else{ |
626 |
|
627 |
fprintf( out_file, |
628 |
"X\t%lf\t%lf\t%lf\n", |
629 |
theDipole->x, |
630 |
theDipole->y, |
631 |
theDipole->z ); |
632 |
} |
633 |
|
634 |
fprintf( out_file, |
635 |
"O\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n", |
636 |
theDipole->x + ox[0], |
637 |
theDipole->y + ox[1], |
638 |
theDipole->z + ox[2] ); |
639 |
|
640 |
fprintf( out_file, |
641 |
"H\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n", |
642 |
theDipole->x + h1[0], |
643 |
theDipole->y + h1[1], |
644 |
theDipole->z + h1[2] ); |
645 |
|
646 |
fprintf( out_file, |
647 |
"H\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n", |
648 |
theDipole->x + h2[0], |
649 |
theDipole->y + h2[1], |
650 |
theDipole->z + h2[2] ); |
651 |
} |
652 |
} |
653 |
|
654 |
else if( !strcmp( "HEAD", theDipole->name )) { |
655 |
|
656 |
if( printDipole ){ |
657 |
|
658 |
fprintf( out_file, |
659 |
"HEAD\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", |
660 |
theDipole->x, |
661 |
theDipole->y, |
662 |
theDipole->z, |
663 |
u[0], u[1], u[2] ); |
664 |
} |
665 |
else{ |
666 |
|
667 |
fprintf( out_file, |
668 |
"HEAD\t%lf\t%lf\t%lf\n", |
669 |
theDipole->x, |
670 |
theDipole->y, |
671 |
theDipole->z ); |
672 |
} |
673 |
} |
674 |
} |
675 |
|
676 |
|
677 |
void setRot( double A[3][3], double the_q[4] ){ |
678 |
|
679 |
double q0Sqr, q1Sqr, q2Sqr, q3Sqr; |
680 |
|
681 |
q0Sqr = the_q[0] * the_q[0]; |
682 |
q1Sqr = the_q[1] * the_q[1]; |
683 |
q2Sqr = the_q[2] * the_q[2]; |
684 |
q3Sqr = the_q[3] * the_q[3]; |
685 |
|
686 |
A[0][0] = q0Sqr + q1Sqr - q2Sqr - q3Sqr; |
687 |
A[0][1] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] ); |
688 |
A[0][2] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] ); |
689 |
|
690 |
A[1][0] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] ); |
691 |
A[1][1] = q0Sqr - q1Sqr + q2Sqr - q3Sqr; |
692 |
A[1][2] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] ); |
693 |
|
694 |
A[2][0] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] ); |
695 |
A[2][1] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] ); |
696 |
A[2][2] = q0Sqr - q1Sqr -q2Sqr +q3Sqr; |
697 |
} |
698 |
|
699 |
|
700 |
|
701 |
void rotMe( double A[3][3], double r[3] ){ |
702 |
|
703 |
int i,j; |
704 |
double rb[3]; // the body frame vector |
705 |
|
706 |
for(i=0; i<3; i++ ){ |
707 |
rb[i] = r[i]; |
708 |
} |
709 |
|
710 |
for( i=0; i<3; i++ ){ |
711 |
r[i] = 0.0; |
712 |
|
713 |
for( j=0; j<3; j++ ){ |
714 |
r[i] += A[j][i] * rb[j]; |
715 |
} |
716 |
} |
717 |
} |
718 |
|
719 |
|
720 |
|
721 |
/*************************************************************************** |
722 |
* prints out the usage for the command line arguments, then exits. |
723 |
***************************************************************************/ |
724 |
|
725 |
void usage(){ |
726 |
(void)fprintf(stderr, |
727 |
"The proper usage is: %s [options] <dumpfile>\n" |
728 |
"\n" |
729 |
"Options:\n" |
730 |
"\n" |
731 |
" -h Display this message\n" |
732 |
" -o <out_name> The output file (Defaults to <dumpfile>.xyz)\n" |
733 |
" -d print the dipole moments\n" |
734 |
" -w print the waters\n" |
735 |
" -m map to the periodic box\n" |
736 |
" -n <#> print every <#> frames\n" |
737 |
"\n" |
738 |
" --repeatX <#> The number of images to repeat in the x direction\n" |
739 |
" --repeatY <#> The number of images to repeat in the y direction\n" |
740 |
" --repeatZ <#> The number of images to repeat in the z direction\n" |
741 |
|
742 |
"\n", |
743 |
|
744 |
program_name); |
745 |
exit(8); |
746 |
} |
747 |
|
748 |
void map( x, y, z, centerX, centerY, centerZ, boxX, boxY, boxZ ) |
749 |
double *x, *y, *z; |
750 |
double centerX, centerY, centerZ; |
751 |
double boxX, boxY, boxZ; |
752 |
{ |
753 |
|
754 |
*x -= centerX; |
755 |
*y -= centerY; |
756 |
*z -= centerZ; |
757 |
|
758 |
if(*x < 0) *x -= boxX * (double)( (int)( (*x / boxX) - 0.5 ) ); |
759 |
else *x -= boxX * (double)( (int)( (*x / boxX ) + 0.5)); |
760 |
|
761 |
if(*y < 0) *y -= boxY * (double)( (int)( (*y / boxY) - 0.5 ) ); |
762 |
else *y -= boxY * (double)( (int)( (*y / boxY ) + 0.5)); |
763 |
|
764 |
if(*z < 0) *z -= boxZ * (double)( (int)( (*z / boxZ) - 0.5 ) ); |
765 |
else *z -= boxZ * (double)( (int)( (*z / boxZ ) + 0.5)); |
766 |
} |