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