2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
59 /* this version only works correctly if one of the entries in the index file
60 is a plane (three atoms specified) and the other a vector. Distance
61 is calculated from the center of the plane to both atoms of the vector */
63 static void print_types(atom_id index1[], int gnx1, char *group1,
64 atom_id index2[], int gnx2, char *group2,
70 fprintf(stderr,"Group %s contains the following atoms: \n",group1);
72 fprintf(stderr,"Atomname %d: %s\n",i,*(top->atoms.atomname[index1[i]]));
75 fprintf(stderr,"Group %s contains the following atoms: \n",group2);
77 fprintf(stderr,"Atomname %d: %s\n",j,*(top->atoms.atomname[index2[j]]));
80 fprintf(stderr,"Careful: distance only makes sense in some situations.\n\n");
83 static void calculate_normal(atom_id index[],rvec x[],rvec result,rvec center)
88 /* calculate centroid of triangle spanned by the three points */
90 center[i] = (x[index[0]][i] + x[index[1]][i] + x[index[2]][i])/3;
92 /* use P1P2 x P1P3 to calculate normal, given three points P1-P3 */
93 rvec_sub(x[index[1]],x[index[0]],c1); /* find two vectors */
94 rvec_sub(x[index[2]],x[index[0]],c2);
96 cprod(c1,c2,result); /* take crossproduct between these */
99 /* calculate the angle and distance between the two groups */
100 static void calc_angle(int ePBC,matrix box,rvec x[], atom_id index1[],
101 atom_id index2[], int gnx1, int gnx2,
102 real *angle, real *distance,
103 real *distance1, real *distance2)
105 /* distance is distance between centers, distance 1 between center of plane
106 and atom one of vector, distance 2 same for atom two
111 normal1,normal2, /* normals on planes of interest */
112 center1,center2, /* center of triangle of points given to define plane,*/
113 /* or center of vector if a vector is given */
114 h1,h2,h3,h4,h5; /* temp. vectors */
117 set_pbc(&pbc,ePBC,box);
121 case 3: /* group 1 defines plane */
122 calculate_normal(index1,x,normal1,center1);
124 case 2: /* group 1 defines vector */
125 rvec_sub(x[index1[0]],x[index1[1]],normal1);
126 rvec_add(x[index1[0]],x[index1[1]],h1);
127 svmul(0.5,h1,center1); /* center is geometric mean */
129 default: /* group 1 does none of the above */
130 gmx_fatal(FARGS,"Something wrong with contents of index file. Groups should contain 2 or 3 atoms.\n");
135 case 3: /* group 2 defines plane */
136 calculate_normal(index2,x,normal2,center2);
138 case 2: /* group 2 defines vector */
139 rvec_sub(x[index2[0]],x[index2[1]],normal2);
140 rvec_add(x[index2[0]],x[index2[1]],h2);
141 svmul(0.5,h2,center2); /* center is geometric mean */
151 default: /* group 2 does none of the above */
152 gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
155 *angle = cos_angle(normal1,normal2);
158 pbc_dx(&pbc,center1,center2,h3);
160 rvec_sub(center1,center2,h3);
161 *distance = norm(h3);
163 if (gnx1 == 3 && gnx2 == 2) {
165 pbc_dx(&pbc,center1,x[index2[0]],h4);
166 pbc_dx(&pbc,center1,x[index2[1]],h5);
168 rvec_sub(center1,x[index2[0]],h4);
169 rvec_sub(center1,x[index2[1]],h5);
171 *distance1 = norm(h4);
172 *distance2 = norm(h5);
174 else if (gnx1 == 2 && gnx2 ==3) {
175 rvec_sub(center1,x[index1[0]],h4);
176 rvec_sub(center1,x[index1[1]],h5);
177 *distance1 = norm(h4);
178 *distance2 = norm(h5);
181 *distance1 = 0; *distance2 = 0;
185 void sgangle_plot(const char *fn,const char *afile,const char *dfile,
186 const char *d1file, const char *d2file,
187 atom_id index1[], int gnx1, char *grpn1,
188 atom_id index2[], int gnx2, char *grpn2,
189 t_topology *top,int ePBC,const output_env_t oenv)
192 *sg_angle, /* xvgr file with angles */
193 *sg_distance = NULL, /* xvgr file with distances */
194 *sg_distance1 = NULL,/* xvgr file with distance between plane and atom */
195 *sg_distance2 = NULL;/* xvgr file with distance between plane and atom2 */
198 angle, /* cosine of angle between two groups */
199 distance, /* distance between two groups. */
200 distance1, /* distance between plane and one of two atoms */
201 distance2; /* same for second of two atoms */
204 rvec *x0; /* coordinates, and coordinates corrected for pb */
206 char buf[256]; /* for xvgr title */
207 gmx_rmpbc_t gpbc=NULL;
210 if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
211 gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
213 sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
214 sg_angle = xvgropen(afile,buf,"Time (ps)","Angle (degrees)",oenv);
217 sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
218 sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
222 sprintf(buf,"Distance between plane and first atom of vector");
223 sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
227 sprintf(buf,"Distance between plane and second atom of vector");
228 sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
231 gpbc = gmx_rmpbc_init(&(top->idef),ePBC,natoms,box);
237 gmx_rmpbc(gpbc,natoms,box,x0);
239 calc_angle(ePBC,box,x0,index1,index2,gnx1,gnx2,&angle,
240 &distance,&distance1,&distance2);
242 fprintf(sg_angle,"%12g %12g %12g\n",t,angle,acos(angle)*180.0/M_PI);
244 fprintf(sg_distance,"%12g %12g\n",t,distance);
246 fprintf(sg_distance1,"%12g %12g\n",t,distance1);
248 fprintf(sg_distance2,"%12g %12g\n",t,distance1);
250 } while (read_next_x(oenv,status,&t,natoms,x0,box));
252 gmx_rmpbc_done(gpbc);
254 fprintf(stderr,"\n");
258 ffclose(sg_distance);
260 ffclose(sg_distance1);
262 ffclose(sg_distance2);
267 static void calc_angle_single(int ePBC,
282 /* distance is distance between centers, distance 1 between center of plane
283 and atom one of vector, distance 2 same for atom two
286 rvec normal1,normal2; /* normals on planes of interest */
287 rvec center1,center2;
288 /* center of triangle of pts to define plane,
289 * or center of vector if a vector is given
291 rvec h1,h2,h3,h4,h5; /* temp. vectors */
294 set_pbc(&pbc,ePBC,box);
297 case 3: /* group 1 defines plane */
298 calculate_normal(index1,xzero,normal1,center1);
300 case 2: /* group 1 defines vector */
301 rvec_sub(xzero[index1[0]],xzero[index1[1]],normal1);
302 rvec_add(xzero[index1[0]],xzero[index1[1]],h1);
303 svmul(0.5,h1,center1); /* center is geometric mean */
305 default: /* group 1 does none of the above */
306 gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
310 case 3: /* group 2 defines plane */
311 calculate_normal(index2,x,normal2,center2);
313 case 2: /* group 2 defines vector */
314 rvec_sub(x[index2[0]],x[index2[1]],normal2);
315 rvec_add(x[index2[0]],x[index2[1]],h2);
316 svmul(0.5,h2,center2); /* center is geometric mean */
318 default: /* group 2 does none of the above */
319 gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
322 *angle = cos_angle(normal1,normal2);
325 pbc_dx(&pbc,center1,center2,h3);
327 rvec_sub(center1,center2,h3);
328 *distance = norm(h3);
330 if (gnx1 == 3 && gnx2 == 2) {
332 pbc_dx(&pbc,center1,x[index2[0]],h4);
333 pbc_dx(&pbc,center1,x[index2[1]],h5);
335 rvec_sub(center1,x[index2[0]],h4);
336 rvec_sub(center1,x[index2[1]],h5);
338 *distance1 = norm(h4);
339 *distance2 = norm(h5);
340 } else if (gnx1 == 2 && gnx2 ==3) {
341 rvec_sub(center1,xzero[index1[0]],h4);
342 rvec_sub(center1,xzero[index1[1]],h5);
343 *distance1 = norm(h4);
344 *distance2 = norm(h5);
346 *distance1 = 0; *distance2 = 0;
351 void sgangle_plot_single(const char *fn,const char *afile,const char *dfile,
352 const char *d1file, const char *d2file,
353 atom_id index1[], int gnx1, char *grpn1,
354 atom_id index2[], int gnx2, char *grpn2,
355 t_topology *top,int ePBC, const output_env_t oenv)
358 *sg_angle, /* xvgr file with angles */
359 *sg_distance = NULL, /* xvgr file with distances */
360 *sg_distance1 = NULL,/* xvgr file with distance between plane and atom */
361 *sg_distance2 = NULL;/* xvgr file with distance between plane and atom2 */
364 angle, /* cosine of angle between two groups */
365 distance, /* distance between two groups. */
366 distance1, /* distance between plane and one of two atoms */
367 distance2; /* same for second of two atoms */
371 rvec *x0; /* coordinates, and coordinates corrected for pb */
374 char buf[256]; /* for xvgr title */
375 gmx_rmpbc_t gpbc=NULL;
378 if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
379 gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
381 sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
382 sg_angle = xvgropen(afile,buf,"Time (ps)","Cos(angle) ",oenv);
385 sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
386 sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
390 sprintf(buf,"Distance between plane and first atom of vector");
391 sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
395 sprintf(buf,"Distance between plane and second atom of vector");
396 sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
400 gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
405 gmx_rmpbc(gpbc,natoms,box,x0);
407 for(i=0;i<natoms;i++)
408 copy_rvec(x0[i],xzero[i]);
412 calc_angle_single(ePBC,box,
413 xzero,x0,index1,index2,gnx1,gnx2,&angle,
414 &distance,&distance1,&distance2);
416 fprintf(sg_angle,"%12g %12g %12g\n",t,angle,acos(angle)*180.0/M_PI);
418 fprintf(sg_distance,"%12g %12g\n",t,distance);
420 fprintf(sg_distance1,"%12g %12g\n",t,distance1);
422 fprintf(sg_distance2,"%12g %12g\n",t,distance1);
424 } while (read_next_x(oenv,status,&t,natoms,x0,box));
425 gmx_rmpbc_done(gpbc);
427 fprintf(stderr,"\n");
431 ffclose(sg_distance);
433 ffclose(sg_distance1);
435 ffclose(sg_distance2);
442 int gmx_sgangle(int argc,char *argv[])
444 const char *desc[] = {
445 "Compute the angle and distance between two groups. ",
446 "The groups are defined by a number of atoms given in an index file and",
447 "may be two or three atoms in size.",
448 "If [TT]-one[tt] is set, only one group should be specified in the index",
449 "file and the angle between this group at time 0 and t will be computed.",
450 "The angles calculated depend on the order in which the atoms are ",
451 "given. Giving, for instance, 5 6 will rotate the vector 5-6 with ",
452 "180 degrees compared to giving 6 5. [PAR]If three atoms are given, ",
453 "the normal on the plane spanned by those three atoms will be",
454 "calculated, using the formula P1P2 x P1P3.",
455 "The cos of the angle is calculated, using the inproduct of the two",
456 "normalized vectors.[PAR]",
457 "Here is what some of the file options do:[BR]",
458 "[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]",
459 "[TT]-od[tt]: Distance between two groups. Distance is taken from the center of one group to the center of the other group.[BR]",
460 "[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]",
461 "[TT]-od2[tt]: For two planes this option has no meaning."
465 const char *fna, *fnd, *fnd1, *fnd2;
466 char * grpname[2]; /* name of the two groups */
467 int gnx[2]; /* size of the two groups */
468 t_topology *top; /* topology */
471 static gmx_bool bOne = FALSE, bZ=FALSE;
473 { "-one", FALSE, etBOOL, {&bOne},
474 "Only one group compute angle between vector at time zero and time t"},
475 { "-z", FALSE, etBOOL, {&bZ},
476 "Use the [IT]z[it]-axis as reference" }
478 #define NPA asize(pa)
480 t_filenm fnm[] = { /* files for g_sgangle */
481 { efTRX, "-f", NULL, ffREAD }, /* trajectory file */
482 { efNDX, NULL, NULL, ffREAD }, /* index file */
483 { efTPX, NULL, NULL, ffREAD }, /* topology file */
484 { efXVG,"-oa","sg_angle",ffWRITE }, /* xvgr output file */
485 { efXVG, "-od","sg_dist",ffOPTWR }, /* xvgr output file */
486 { efXVG, "-od1", "sg_dist1",ffOPTWR }, /* xvgr output file */
487 { efXVG, "-od2", "sg_dist2",ffOPTWR } /* xvgr output file */
490 #define NFILE asize(fnm)
492 CopyRight(stderr,argv[0]);
493 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
494 NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
497 top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC); /* read topology file */
499 fna = opt2fn("-oa",NFILE,fnm);
500 fnd = opt2fn_null("-od",NFILE,fnm);
501 fnd1 = opt2fn_null("-od1",NFILE,fnm);
502 fnd2 = opt2fn_null("-od2",NFILE,fnm);
504 /* read index file. */
506 rd_index(ftp2fn(efNDX,NFILE,fnm),1,gnx,index,grpname);
507 print_types(index[0],gnx[0],grpname[0],
508 index[0],gnx[0],grpname[0],top);
510 sgangle_plot_single(ftp2fn(efTRX,NFILE,fnm), fna, fnd, fnd1, fnd2,
511 index[0],gnx[0],grpname[0],
512 index[0],gnx[0],grpname[0],
515 rd_index(ftp2fn(efNDX,NFILE,fnm),bZ ? 1 : 2,gnx,index,grpname);
517 print_types(index[0],gnx[0],grpname[0],
518 index[1],gnx[1],grpname[1],top);
521 grpname[1] = strdup("Z-axis");
523 sgangle_plot(ftp2fn(efTRX,NFILE,fnm), fna, fnd, fnd1, fnd2,
524 index[0],gnx[0],grpname[0],
525 index[1],gnx[1],grpname[1],
529 do_view(oenv,fna,"-nxy"); /* view xvgr file */
530 do_view(oenv,fnd,"-nxy"); /* view xvgr file */
531 do_view(oenv,fnd1,"-nxy"); /* view xvgr file */
532 do_view(oenv,fnd2,"-nxy"); /* view xvgr file */