3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
56 /* this version only works correctly if one of the entries in the index file
57 is a plane (three atoms specified) and the other a vector. Distance
58 is calculated from the center of the plane to both atoms of the vector */
60 static void print_types(atom_id index1[], int gnx1, char *group1,
61 atom_id index2[], int gnx2, char *group2,
67 fprintf(stderr,"Group %s contains the following atoms: \n",group1);
69 fprintf(stderr,"Atomname %d: %s\n",i,*(top->atoms.atomname[index1[i]]));
72 fprintf(stderr,"Group %s contains the following atoms: \n",group2);
74 fprintf(stderr,"Atomname %d: %s\n",j,*(top->atoms.atomname[index2[j]]));
77 fprintf(stderr,"Careful: distance only makes sense in some situations.\n\n");
80 static void calculate_normal(atom_id index[],rvec x[],rvec result,rvec center)
85 /* calculate centroid of triangle spanned by the three points */
87 center[i] = (x[index[0]][i] + x[index[1]][i] + x[index[2]][i])/3;
89 /* use P1P2 x P1P3 to calculate normal, given three points P1-P3 */
90 rvec_sub(x[index[1]],x[index[0]],c1); /* find two vectors */
91 rvec_sub(x[index[2]],x[index[0]],c2);
93 cprod(c1,c2,result); /* take crossproduct between these */
96 /* calculate the angle and distance between the two groups */
97 static void calc_angle(int ePBC,matrix box,rvec x[], atom_id index1[],
98 atom_id index2[], int gnx1, int gnx2,
99 real *angle, real *distance,
100 real *distance1, real *distance2)
102 /* distance is distance between centers, distance 1 between center of plane
103 and atom one of vector, distance 2 same for atom two
108 normal1,normal2, /* normals on planes of interest */
109 center1,center2, /* center of triangle of points given to define plane,*/
110 /* or center of vector if a vector is given */
111 h1,h2,h3,h4,h5; /* temp. vectors */
114 set_pbc(&pbc,ePBC,box);
118 case 3: /* group 1 defines plane */
119 calculate_normal(index1,x,normal1,center1);
121 case 2: /* group 1 defines vector */
122 rvec_sub(x[index1[0]],x[index1[1]],normal1);
123 rvec_add(x[index1[0]],x[index1[1]],h1);
124 svmul(0.5,h1,center1); /* center is geometric mean */
126 default: /* group 1 does none of the above */
127 gmx_fatal(FARGS,"Something wrong with contents of index file. Groups should contain 2 or 3 atoms.\n");
132 case 3: /* group 2 defines plane */
133 calculate_normal(index2,x,normal2,center2);
135 case 2: /* group 2 defines vector */
136 rvec_sub(x[index2[0]],x[index2[1]],normal2);
137 rvec_add(x[index2[0]],x[index2[1]],h2);
138 svmul(0.5,h2,center2); /* center is geometric mean */
148 default: /* group 2 does none of the above */
149 gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
152 *angle = cos_angle(normal1,normal2);
155 pbc_dx(&pbc,center1,center2,h3);
157 rvec_sub(center1,center2,h3);
158 *distance = norm(h3);
160 if (gnx1 == 3 && gnx2 == 2) {
162 pbc_dx(&pbc,center1,x[index2[0]],h4);
163 pbc_dx(&pbc,center1,x[index2[1]],h5);
165 rvec_sub(center1,x[index2[0]],h4);
166 rvec_sub(center1,x[index2[1]],h5);
168 *distance1 = norm(h4);
169 *distance2 = norm(h5);
171 else if (gnx1 == 2 && gnx2 ==3) {
172 rvec_sub(center1,x[index1[0]],h4);
173 rvec_sub(center1,x[index1[1]],h5);
174 *distance1 = norm(h4);
175 *distance2 = norm(h5);
178 *distance1 = 0; *distance2 = 0;
182 void sgangle_plot(const char *fn,const char *afile,const char *dfile,
183 const char *d1file, const char *d2file,
184 atom_id index1[], int gnx1, char *grpn1,
185 atom_id index2[], int gnx2, char *grpn2,
186 t_topology *top,int ePBC,const output_env_t oenv)
189 *sg_angle, /* xvgr file with angles */
190 *sg_distance = NULL, /* xvgr file with distances */
191 *sg_distance1 = NULL,/* xvgr file with distance between plane and atom */
192 *sg_distance2 = NULL;/* xvgr file with distance between plane and atom2 */
195 angle, /* cosine of angle between two groups */
196 distance, /* distance between two groups. */
197 distance1, /* distance between plane and one of two atoms */
198 distance2; /* same for second of two atoms */
201 rvec *x0; /* coordinates, and coordinates corrected for pb */
203 char buf[256]; /* for xvgr title */
204 gmx_rmpbc_t gpbc=NULL;
207 if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
208 gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
210 sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
211 sg_angle = xvgropen(afile,buf,"Time (ps)","Angle (degrees)",oenv);
214 sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
215 sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
219 sprintf(buf,"Distance between plane and first atom of vector");
220 sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
224 sprintf(buf,"Distance between plane and second atom of vector");
225 sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
228 gpbc = gmx_rmpbc_init(&(top->idef),ePBC,natoms,box);
234 gmx_rmpbc(gpbc,natoms,box,x0);
236 calc_angle(ePBC,box,x0,index1,index2,gnx1,gnx2,&angle,
237 &distance,&distance1,&distance2);
239 fprintf(sg_angle,"%12g %12g %12g\n",t,angle,acos(angle)*180.0/M_PI);
241 fprintf(sg_distance,"%12g %12g\n",t,distance);
243 fprintf(sg_distance1,"%12g %12g\n",t,distance1);
245 fprintf(sg_distance2,"%12g %12g\n",t,distance1);
247 } while (read_next_x(oenv,status,&t,natoms,x0,box));
249 gmx_rmpbc_done(gpbc);
251 fprintf(stderr,"\n");
255 ffclose(sg_distance);
257 ffclose(sg_distance1);
259 ffclose(sg_distance2);
264 static void calc_angle_single(int ePBC,
279 /* distance is distance between centers, distance 1 between center of plane
280 and atom one of vector, distance 2 same for atom two
283 rvec normal1,normal2; /* normals on planes of interest */
284 rvec center1,center2;
285 /* center of triangle of pts to define plane,
286 * or center of vector if a vector is given
288 rvec h1,h2,h3,h4,h5; /* temp. vectors */
291 set_pbc(&pbc,ePBC,box);
294 case 3: /* group 1 defines plane */
295 calculate_normal(index1,xzero,normal1,center1);
297 case 2: /* group 1 defines vector */
298 rvec_sub(xzero[index1[0]],xzero[index1[1]],normal1);
299 rvec_add(xzero[index1[0]],xzero[index1[1]],h1);
300 svmul(0.5,h1,center1); /* center is geometric mean */
302 default: /* group 1 does none of the above */
303 gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
307 case 3: /* group 2 defines plane */
308 calculate_normal(index2,x,normal2,center2);
310 case 2: /* group 2 defines vector */
311 rvec_sub(x[index2[0]],x[index2[1]],normal2);
312 rvec_add(x[index2[0]],x[index2[1]],h2);
313 svmul(0.5,h2,center2); /* center is geometric mean */
315 default: /* group 2 does none of the above */
316 gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
319 *angle = cos_angle(normal1,normal2);
322 pbc_dx(&pbc,center1,center2,h3);
324 rvec_sub(center1,center2,h3);
325 *distance = norm(h3);
327 if (gnx1 == 3 && gnx2 == 2) {
329 pbc_dx(&pbc,center1,x[index2[0]],h4);
330 pbc_dx(&pbc,center1,x[index2[1]],h5);
332 rvec_sub(center1,x[index2[0]],h4);
333 rvec_sub(center1,x[index2[1]],h5);
335 *distance1 = norm(h4);
336 *distance2 = norm(h5);
337 } else if (gnx1 == 2 && gnx2 ==3) {
338 rvec_sub(center1,xzero[index1[0]],h4);
339 rvec_sub(center1,xzero[index1[1]],h5);
340 *distance1 = norm(h4);
341 *distance2 = norm(h5);
343 *distance1 = 0; *distance2 = 0;
348 void sgangle_plot_single(const char *fn,const char *afile,const char *dfile,
349 const char *d1file, const char *d2file,
350 atom_id index1[], int gnx1, char *grpn1,
351 atom_id index2[], int gnx2, char *grpn2,
352 t_topology *top,int ePBC, const output_env_t oenv)
355 *sg_angle, /* xvgr file with angles */
356 *sg_distance = NULL, /* xvgr file with distances */
357 *sg_distance1 = NULL,/* xvgr file with distance between plane and atom */
358 *sg_distance2 = NULL;/* xvgr file with distance between plane and atom2 */
361 angle, /* cosine of angle between two groups */
362 distance, /* distance between two groups. */
363 distance1, /* distance between plane and one of two atoms */
364 distance2; /* same for second of two atoms */
368 rvec *x0; /* coordinates, and coordinates corrected for pb */
371 char buf[256]; /* for xvgr title */
372 gmx_rmpbc_t gpbc=NULL;
375 if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
376 gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
378 sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
379 sg_angle = xvgropen(afile,buf,"Time (ps)","Cos(angle) ",oenv);
382 sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
383 sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
387 sprintf(buf,"Distance between plane and first atom of vector");
388 sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
392 sprintf(buf,"Distance between plane and second atom of vector");
393 sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
397 gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
402 gmx_rmpbc(gpbc,natoms,box,x0);
404 for(i=0;i<natoms;i++)
405 copy_rvec(x0[i],xzero[i]);
409 calc_angle_single(ePBC,box,
410 xzero,x0,index1,index2,gnx1,gnx2,&angle,
411 &distance,&distance1,&distance2);
413 fprintf(sg_angle,"%12g %12g %12g\n",t,angle,acos(angle)*180.0/M_PI);
415 fprintf(sg_distance,"%12g %12g\n",t,distance);
417 fprintf(sg_distance1,"%12g %12g\n",t,distance1);
419 fprintf(sg_distance2,"%12g %12g\n",t,distance1);
421 } while (read_next_x(oenv,status,&t,natoms,x0,box));
422 gmx_rmpbc_done(gpbc);
424 fprintf(stderr,"\n");
428 ffclose(sg_distance);
430 ffclose(sg_distance1);
432 ffclose(sg_distance2);
439 int gmx_sgangle(int argc,char *argv[])
441 const char *desc[] = {
442 "Compute the angle and distance between two groups. ",
443 "The groups are defined by a number of atoms given in an index file and",
444 "may be two or three atoms in size.",
445 "If [TT]-one[tt] is set, only one group should be specified in the index",
446 "file and the angle between this group at time 0 and t will be computed.",
447 "The angles calculated depend on the order in which the atoms are ",
448 "given. Giving, for instance, 5 6 will rotate the vector 5-6 with ",
449 "180 degrees compared to giving 6 5. [PAR]If three atoms are given, ",
450 "the normal on the plane spanned by those three atoms will be",
451 "calculated, using the formula P1P2 x P1P3.",
452 "The cos of the angle is calculated, using the inproduct of the two",
453 "normalized vectors.[PAR]",
454 "Here is what some of the file options do:[BR]",
455 "[TT]-oa[tt]: Angle between the two groups specified in the index file. If a group contains three atoms the normal to the plane defined by those three atoms will be used. If a group contains two atoms, the vector defined by those two atoms will be used.[BR]",
456 "[TT]-od[tt]: Distance between two groups. Distance is taken from the center of one group to the center of the other group.[BR]",
457 "[TT]-od1[tt]: If one plane and one vector is given, the distances for each of the atoms from the center of the plane is given separately.[BR]",
458 "[TT]-od2[tt]: For two planes this option has no meaning."
462 const char *fna, *fnd, *fnd1, *fnd2;
463 char * grpname[2]; /* name of the two groups */
464 int gnx[2]; /* size of the two groups */
465 t_topology *top; /* topology */
468 static gmx_bool bOne = FALSE, bZ=FALSE;
470 { "-one", FALSE, etBOOL, {&bOne},
471 "Only one group compute angle between vector at time zero and time t"},
472 { "-z", FALSE, etBOOL, {&bZ},
473 "Use the [IT]z[it]-axis as reference" }
475 #define NPA asize(pa)
477 t_filenm fnm[] = { /* files for g_sgangle */
478 { efTRX, "-f", NULL, ffREAD }, /* trajectory file */
479 { efNDX, NULL, NULL, ffREAD }, /* index file */
480 { efTPX, NULL, NULL, ffREAD }, /* topology file */
481 { efXVG,"-oa","sg_angle",ffWRITE }, /* xvgr output file */
482 { efXVG, "-od","sg_dist",ffOPTWR }, /* xvgr output file */
483 { efXVG, "-od1", "sg_dist1",ffOPTWR }, /* xvgr output file */
484 { efXVG, "-od2", "sg_dist2",ffOPTWR } /* xvgr output file */
487 #define NFILE asize(fnm)
489 CopyRight(stderr,argv[0]);
490 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
491 NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
494 top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC); /* read topology file */
496 fna = opt2fn("-oa",NFILE,fnm);
497 fnd = opt2fn_null("-od",NFILE,fnm);
498 fnd1 = opt2fn_null("-od1",NFILE,fnm);
499 fnd2 = opt2fn_null("-od2",NFILE,fnm);
501 /* read index file. */
503 rd_index(ftp2fn(efNDX,NFILE,fnm),1,gnx,index,grpname);
504 print_types(index[0],gnx[0],grpname[0],
505 index[0],gnx[0],grpname[0],top);
507 sgangle_plot_single(ftp2fn(efTRX,NFILE,fnm), fna, fnd, fnd1, fnd2,
508 index[0],gnx[0],grpname[0],
509 index[0],gnx[0],grpname[0],
512 rd_index(ftp2fn(efNDX,NFILE,fnm),bZ ? 1 : 2,gnx,index,grpname);
514 print_types(index[0],gnx[0],grpname[0],
515 index[1],gnx[1],grpname[1],top);
518 grpname[1] = strdup("Z-axis");
520 sgangle_plot(ftp2fn(efTRX,NFILE,fnm), fna, fnd, fnd1, fnd2,
521 index[0],gnx[0],grpname[0],
522 index[1],gnx[1],grpname[1],
526 do_view(oenv,fna,"-nxy"); /* view xvgr file */
527 do_view(oenv,fnd,"-nxy"); /* view xvgr file */
528 do_view(oenv,fnd1,"-nxy"); /* view xvgr file */
529 do_view(oenv,fnd2,"-nxy"); /* view xvgr file */