1 |
gezelter |
444 |
#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 |
|
|
} |