src/tools/gmx_sgangle.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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,2013, 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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41 #include <math.h>
  42
  43 #include "sysstuff.h"
  44 #include <string.h>
  45 #include "typedefs.h"
  46 #include "smalloc.h"
  47 #include "macros.h"
  48 #include "rmpbc.h"
  49 #include "vec.h"
  50 #include "xvgr.h"
  51 #include "pbc.h"
  52 #include "copyrite.h"
  53 #include "futil.h"
  54 #include "statutil.h"
  55 #include "index.h"
  56 #include "gmx_ana.h"
  57
  58
  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 */
  62
  63 static void print_types(atom_id index1[], int gnx1, char *group1,
  64                        atom_id index2[], int gnx2, char *group2,
  65                        t_topology *top)
  66 {
  67   int i,j;
  68
  69   fprintf(stderr,"\n");
  70   fprintf(stderr,"Group %s contains the following atoms: \n",group1);
  71   for(i=0;i<gnx1;i++)
  72     fprintf(stderr,"Atomname %d: %s\n",i,*(top->atoms.atomname[index1[i]]));
  73   fprintf(stderr,"\n");
  74
  75   fprintf(stderr,"Group %s contains the following atoms: \n",group2);
  76   for(j=0;j<gnx2;j++)
  77     fprintf(stderr,"Atomname %d: %s\n",j,*(top->atoms.atomname[index2[j]]));
  78   fprintf(stderr,"\n");
  79
  80   fprintf(stderr,"Careful: distance only makes sense in some situations.\n\n");
  81 }
  82
  83 static void calculate_normal(atom_id index[],rvec x[],rvec result,rvec center)
  84 {
  85   rvec c1,c2;
  86   int i;
  87
  88   /* calculate centroid of triangle spanned by the three points */
  89   for(i=0;i<3;i++)
  90     center[i] = (x[index[0]][i] + x[index[1]][i] + x[index[2]][i])/3;
  91
  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);
  95
  96   cprod(c1,c2,result);                    /* take crossproduct between these */
  97 }
  98
  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)
 104
 105 /* distance is distance between centers, distance 1 between center of plane
 106    and atom one of vector, distance 2 same for atom two
 107 */
 108
 109 {
 110   rvec
 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 */
 115   t_pbc pbc;
 116
 117   set_pbc(&pbc,ePBC,box);
 118
 119   switch(gnx1)
 120     {
 121     case 3:           /* group 1 defines plane */
 122       calculate_normal(index1,x,normal1,center1);
 123       break;
 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 */
 128       break;
 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");
 131     }
 132
 133   switch(gnx2)
 134     {
 135     case 3:          /* group 2 defines plane */
 136       calculate_normal(index2,x,normal2,center2);
 137       break;
 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 */
 142       break;
 143     case 0:
 144       normal2[XX] = 0;
 145       normal2[YY] = 0;
 146       normal2[ZZ] = 1;
 147       center2[XX] = 0;
 148       center2[YY] = 0;
 149       center2[ZZ] = 0;
 150       break;
 151     default:         /* group 2 does none of the above */
 152       gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
 153     }
 154
 155   *angle = cos_angle(normal1,normal2);
 156
 157   if (box)
 158     pbc_dx(&pbc,center1,center2,h3);
 159   else
 160     rvec_sub(center1,center2,h3);
 161   *distance = norm(h3);
 162
 163   if (gnx1 == 3 && gnx2 == 2) {
 164     if (box) {
 165       pbc_dx(&pbc,center1,x[index2[0]],h4);
 166       pbc_dx(&pbc,center1,x[index2[1]],h5);
 167     } else {
 168       rvec_sub(center1,x[index2[0]],h4);
 169       rvec_sub(center1,x[index2[1]],h5);
 170     }
 171     *distance1 = norm(h4);
 172     *distance2 = norm(h5);
 173   }
 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);
 179   }
 180   else {
 181     *distance1 = 0; *distance2 = 0;
 182   }
 183 }
 184
 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)
 190 {
 191   FILE
 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 */
 196   real
 197     t,                   /* time */
 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 */
 202   t_trxstatus *status;
 203   int        natoms,teller=0;
 204   rvec       *x0;   /* coordinates, and coordinates corrected for pb */
 205   matrix     box;
 206   char       buf[256];   /* for xvgr title */
 207   gmx_rmpbc_t  gpbc=NULL;
 208
 209
 210   if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
 211     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 212
 213   sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
 214   sg_angle = xvgropen(afile,buf,"Time (ps)","Angle (degrees)",oenv);
 215
 216   if (dfile) {
 217     sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
 218     sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
 219   }
 220
 221   if (d1file) {
 222     sprintf(buf,"Distance between plane and first atom of vector");
 223     sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
 224   }
 225
 226   if (d2file) {
 227     sprintf(buf,"Distance between plane and second atom of vector");
 228     sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
 229   }
 230
 231   gpbc = gmx_rmpbc_init(&(top->idef),ePBC,natoms,box);
 232
 233   do
 234     {
 235       teller++;
 236
 237       gmx_rmpbc(gpbc,natoms,box,x0);
 238
 239       calc_angle(ePBC,box,x0,index1,index2,gnx1,gnx2,&angle,
 240                  &distance,&distance1,&distance2);
 241
 242       fprintf(sg_angle,"%12g  %12g  %12g\n",t,angle,acos(angle)*180.0/M_PI);
 243       if (dfile)
 244         fprintf(sg_distance,"%12g  %12g\n",t,distance);
 245       if (d1file)
 246         fprintf(sg_distance1,"%12g  %12g\n",t,distance1);
 247       if (d2file)
 248         fprintf(sg_distance2,"%12g  %12g\n",t,distance1);
 249
 250     } while (read_next_x(oenv,status,&t,natoms,x0,box));
 251
 252   gmx_rmpbc_done(gpbc);
 253
 254   fprintf(stderr,"\n");
 255   close_trj(status);
 256   ffclose(sg_angle);
 257   if (dfile)
 258     ffclose(sg_distance);
 259   if (d1file)
 260     ffclose(sg_distance1);
 261   if (d2file)
 262     ffclose(sg_distance2);
 263
 264   sfree(x0);
 265 }
 266
 267 static void calc_angle_single(int ePBC,
 268                               matrix box,
 269                               rvec xzero[],
 270                               rvec x[],
 271                               atom_id index1[],
 272                               atom_id index2[],
 273                               int gnx1,
 274                               int gnx2,
 275                               real *angle,
 276                               real *distance,
 277                               real *distance1,
 278                               real *distance2)
 279 {
 280   t_pbc pbc;
 281
 282   /* distance is distance between centers, distance 1 between center of plane
 283      and atom one of vector, distance 2 same for atom two
 284   */
 285
 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
 290      */
 291   rvec  h1,h2,h3,h4,h5;         /* temp. vectors */
 292
 293   if (box)
 294     set_pbc(&pbc,ePBC,box);
 295
 296   switch(gnx1) {
 297   case 3:           /* group 1 defines plane */
 298     calculate_normal(index1,xzero,normal1,center1);
 299     break;
 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 */
 304     break;
 305   default:          /* group 1 does none of the above */
 306     gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
 307   }
 308
 309   switch(gnx2) {
 310   case 3:          /* group 2 defines plane */
 311     calculate_normal(index2,x,normal2,center2);
 312     break;
 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 */
 317     break;
 318   default:         /* group 2 does none of the above */
 319     gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
 320   }
 321
 322   *angle = cos_angle(normal1,normal2);
 323
 324   if (box)
 325     pbc_dx(&pbc,center1,center2,h3);
 326   else
 327     rvec_sub(center1,center2,h3);
 328   *distance = norm(h3);
 329
 330   if (gnx1 == 3 && gnx2 == 2) {
 331     if (box) {
 332       pbc_dx(&pbc,center1,x[index2[0]],h4);
 333       pbc_dx(&pbc,center1,x[index2[1]],h5);
 334     } else {
 335       rvec_sub(center1,x[index2[0]],h4);
 336       rvec_sub(center1,x[index2[1]],h5);
 337     }
 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);
 345   } else {
 346     *distance1 = 0; *distance2 = 0;
 347   }
 348 }
 349
 350
 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)
 356 {
 357   FILE
 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 */
 362   real
 363     t,                   /* time */
 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 */
 368   t_trxstatus *status;
 369   int        natoms,teller=0;
 370   int        i;
 371   rvec       *x0;   /* coordinates, and coordinates corrected for pb */
 372   rvec       *xzero;
 373   matrix     box;
 374   char       buf[256];   /* for xvgr title */
 375   gmx_rmpbc_t  gpbc=NULL;
 376
 377
 378   if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
 379     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 380
 381   sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
 382   sg_angle = xvgropen(afile,buf,"Time (ps)","Cos(angle) ",oenv);
 383
 384   if (dfile) {
 385     sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
 386     sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
 387   }
 388
 389   if (d1file) {
 390     sprintf(buf,"Distance between plane and first atom of vector");
 391     sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
 392   }
 393
 394   if (d2file) {
 395     sprintf(buf,"Distance between plane and second atom of vector");
 396     sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
 397   }
 398
 399   snew(xzero,natoms);
 400   gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
 401
 402   do {
 403     teller++;
 404
 405     gmx_rmpbc(gpbc,natoms,box,x0);
 406     if (teller==1) {
 407       for(i=0;i<natoms;i++)
 408         copy_rvec(x0[i],xzero[i]);
 409     }
 410
 411
 412     calc_angle_single(ePBC,box,
 413                       xzero,x0,index1,index2,gnx1,gnx2,&angle,
 414                       &distance,&distance1,&distance2);
 415
 416     fprintf(sg_angle,"%12g  %12g  %12g\n",t,angle,acos(angle)*180.0/M_PI);
 417     if (dfile)
 418       fprintf(sg_distance,"%12g  %12g\n",t,distance);
 419     if (d1file)
 420       fprintf(sg_distance1,"%12g  %12g\n",t,distance1);
 421     if (d2file)
 422       fprintf(sg_distance2,"%12g  %12g\n",t,distance1);
 423
 424   } while (read_next_x(oenv,status,&t,natoms,x0,box));
 425   gmx_rmpbc_done(gpbc);
 426
 427   fprintf(stderr,"\n");
 428   close_trj(status);
 429   ffclose(sg_angle);
 430   if (dfile)
 431     ffclose(sg_distance);
 432   if (d1file)
 433     ffclose(sg_distance1);
 434   if (d2file)
 435     ffclose(sg_distance2);
 436
 437   sfree(x0);
 438 }
 439
 440
 441
 442 int gmx_sgangle(int argc,char *argv[])
 443 {
 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."
 462   };
 463
 464   output_env_t oenv;
 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             */
 469   int       ePBC;
 470   atom_id   *index[2];
 471   static gmx_bool bOne = FALSE, bZ=FALSE;
 472   t_pargs pa[] = {
 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" }
 477   };
 478 #define NPA asize(pa)
 479
 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     */
 488   };
 489
 490 #define NFILE asize(fnm)
 491
 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);
 495
 496
 497   top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);     /* read topology file */
 498
 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);
 503
 504   /* read index file. */
 505   if(bOne) {
 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);
 509
 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],
 513                         top,ePBC,oenv);
 514   }  else {
 515     rd_index(ftp2fn(efNDX,NFILE,fnm),bZ ? 1 : 2,gnx,index,grpname);
 516     if (!bZ)
 517       print_types(index[0],gnx[0],grpname[0],
 518                 index[1],gnx[1],grpname[1],top);
 519     else {
 520       gnx[1] = 0;
 521       grpname[1] = strdup("Z-axis");
 522     }
 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],
 526                  top,ePBC,oenv);
 527   }
 528
 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 */
 533
 534   thanx(stderr);
 535   return 0;
 536 }