src/tools/gmx_sgangle.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  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.
  14
  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.
  19  *
  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.
  26  *
  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.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <math.h>
  39
  40 #include "sysstuff.h"
  41 #include "string.h"
  42 #include "typedefs.h"
  43 #include "smalloc.h"
  44 #include "macros.h"
  45 #include "rmpbc.h"
  46 #include "vec.h"
  47 #include "xvgr.h"
  48 #include "pbc.h"
  49 #include "copyrite.h"
  50 #include "futil.h"
  51 #include "statutil.h"
  52 #include "index.h"
  53
  54 /* this version only works correctly if one of the entries in the index file
  55    is a plane (three atoms specified) and the other a vector. Distance
  56    is calculated from the center of the plane to both atoms of the vector */
  57
  58 static void print_types(atom_id index1[], int gnx1, char *group1,
  59                        atom_id index2[], int gnx2, char *group2,
  60                        t_topology *top)
  61 {
  62   int i,j;
  63
  64   fprintf(stderr,"\n");
  65   fprintf(stderr,"Group %s contains the following atoms: \n",group1);
  66   for(i=0;i<gnx1;i++)
  67     fprintf(stderr,"Atomname %d: %s\n",i,*(top->atoms.atomname[index1[i]]));
  68   fprintf(stderr,"\n");
  69
  70   fprintf(stderr,"Group %s contains the following atoms: \n",group2);
  71   for(j=0;j<gnx2;j++)
  72     fprintf(stderr,"Atomname %d: %s\n",j,*(top->atoms.atomname[index2[j]]));
  73   fprintf(stderr,"\n");
  74
  75   fprintf(stderr,"Careful: distance only makes sense in some situations.\n\n");
  76 }
  77
  78 static void calculate_normal(atom_id index[],rvec x[],rvec result,rvec center)
  79 {
  80   rvec c1,c2;
  81   int i;
  82
  83   /* calculate centroid of triangle spanned by the three points */
  84   for(i=0;i<3;i++)
  85     center[i] = (x[index[0]][i] + x[index[1]][i] + x[index[2]][i])/3;
  86
  87   /* use P1P2 x P1P3 to calculate normal, given three points P1-P3 */
  88   rvec_sub(x[index[1]],x[index[0]],c1);    /* find two vectors */
  89   rvec_sub(x[index[2]],x[index[0]],c2);
  90
  91   cprod(c1,c2,result);                    /* take crossproduct between these */
  92 }
  93
  94 /* calculate the angle and distance between the two groups */
  95 static void calc_angle(int ePBC,matrix box,rvec x[], atom_id index1[],
  96                        atom_id index2[], int gnx1, int gnx2,
  97                        real *angle,      real *distance,
  98                        real *distance1,  real *distance2)
  99
 100 /* distance is distance between centers, distance 1 between center of plane
 101    and atom one of vector, distance 2 same for atom two
 102 */
 103
 104 {
 105   rvec
 106     normal1,normal2,    /* normals on planes of interest */
 107     center1,center2,    /* center of triangle of points given to define plane,*/
 108                         /* or center of vector if a vector is given */
 109     h1,h2,h3,h4,h5;     /* temp. vectors */
 110   t_pbc pbc;
 111
 112   set_pbc(&pbc,ePBC,box);
 113
 114   switch(gnx1)
 115     {
 116     case 3:           /* group 1 defines plane */
 117       calculate_normal(index1,x,normal1,center1);
 118       break;
 119     case 2:           /* group 1 defines vector */
 120       rvec_sub(x[index1[0]],x[index1[1]],normal1);
 121       rvec_add(x[index1[0]],x[index1[1]],h1);
 122       svmul(0.5,h1,center1);  /* center is geometric mean */
 123       break;
 124     default:          /* group 1 does none of the above */
 125       gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
 126     }
 127
 128   switch(gnx2)
 129     {
 130     case 3:          /* group 2 defines plane */
 131       calculate_normal(index2,x,normal2,center2);
 132       break;
 133     case 2:          /* group 2 defines vector */
 134       rvec_sub(x[index2[0]],x[index2[1]],normal2);
 135       rvec_add(x[index2[0]],x[index2[1]],h2);
 136       svmul(0.5,h2,center2);  /* center is geometric mean */
 137       break;
 138     case 0:
 139       normal2[XX] = 0;
 140       normal2[YY] = 0;
 141       normal2[ZZ] = 1;
 142       center2[XX] = 0;
 143       center2[YY] = 0;
 144       center2[ZZ] = 0;
 145       break;
 146     default:         /* group 2 does none of the above */
 147       gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
 148     }
 149
 150   *angle = cos_angle(normal1,normal2);
 151
 152   if (box)
 153     pbc_dx(&pbc,center1,center2,h3);
 154   else
 155     rvec_sub(center1,center2,h3);
 156   *distance = norm(h3);
 157
 158   if (gnx1 == 3 && gnx2 == 2) {
 159     if (box) {
 160       pbc_dx(&pbc,center1,x[index2[0]],h4);
 161       pbc_dx(&pbc,center1,x[index2[1]],h5);
 162     } else {
 163       rvec_sub(center1,x[index2[0]],h4);
 164       rvec_sub(center1,x[index2[1]],h5);
 165     }
 166     *distance1 = norm(h4);
 167     *distance2 = norm(h5);
 168   }
 169   else if (gnx1 == 2 && gnx2 ==3) {
 170     rvec_sub(center1,x[index1[0]],h4);
 171     rvec_sub(center1,x[index1[1]],h5);
 172     *distance1 = norm(h4);
 173     *distance2 = norm(h5);
 174   }
 175   else {
 176     *distance1 = 0; *distance2 = 0;
 177   }
 178 }
 179
 180 void sgangle_plot(const char *fn,const char *afile,const char *dfile,
 181                   const char *d1file, const char *d2file,
 182                   atom_id index1[], int gnx1, char *grpn1,
 183                   atom_id index2[], int gnx2, char *grpn2,
 184                   t_topology *top,int ePBC,const output_env_t oenv)
 185 {
 186   FILE
 187     *sg_angle,           /* xvgr file with angles */
 188     *sg_distance = NULL, /* xvgr file with distances */
 189     *sg_distance1 = NULL,/* xvgr file with distance between plane and atom */
 190     *sg_distance2 = NULL;/* xvgr file with distance between plane and atom2 */
 191   real
 192     t,                   /* time */
 193     angle,               /* cosine of angle between two groups */
 194     distance,            /* distance between two groups. */
 195     distance1,           /* distance between plane and one of two atoms */
 196     distance2;           /* same for second of two atoms */
 197   int        status,natoms,teller=0;
 198   rvec       *x0;   /* coordinates, and coordinates corrected for pb */
 199   matrix     box;
 200   char       buf[256];   /* for xvgr title */
 201
 202   if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
 203     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 204
 205   sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
 206   sg_angle = xvgropen(afile,buf,"Time (ps)","Angle (degrees)",oenv);
 207
 208   if (dfile) {
 209     sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
 210     sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
 211   }
 212
 213   if (d1file) {
 214     sprintf(buf,"Distance between plane and first atom of vector");
 215     sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
 216   }
 217
 218   if (d2file) {
 219     sprintf(buf,"Distance between plane and second atom of vector");
 220     sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
 221   }
 222
 223   do
 224     {
 225       teller++;
 226
 227       rm_pbc(&(top->idef),ePBC,natoms,box,x0,x0);
 228
 229       calc_angle(ePBC,box,x0,index1,index2,gnx1,gnx2,&angle,
 230                  &distance,&distance1,&distance2);
 231
 232       fprintf(sg_angle,"%12g  %12g  %12g\n",t,angle,acos(angle)*180.0/M_PI);
 233       if (dfile)
 234         fprintf(sg_distance,"%12g  %12g\n",t,distance);
 235       if (d1file)
 236         fprintf(sg_distance1,"%12g  %12g\n",t,distance1);
 237       if (d2file)
 238         fprintf(sg_distance2,"%12g  %12g\n",t,distance1);
 239
 240     } while (read_next_x(oenv,status,&t,natoms,x0,box));
 241
 242   fprintf(stderr,"\n");
 243   close_trj(status);
 244   fclose(sg_angle);
 245   if (dfile)
 246     fclose(sg_distance);
 247   if (d1file)
 248     fclose(sg_distance1);
 249   if (d2file)
 250     fclose(sg_distance2);
 251
 252   sfree(x0);
 253 }
 254
 255 static void calc_angle_single(int ePBC,
 256                               matrix box,
 257                               rvec xzero[],
 258                               rvec x[],
 259                               atom_id index1[],
 260                               atom_id index2[],
 261                               int gnx1,
 262                               int gnx2,
 263                               real *angle,
 264                               real *distance,
 265                               real *distance1,
 266                               real *distance2)
 267 {
 268   t_pbc pbc;
 269
 270   /* distance is distance between centers, distance 1 between center of plane
 271      and atom one of vector, distance 2 same for atom two
 272   */
 273
 274   rvec  normal1,normal2;        /* normals on planes of interest */
 275   rvec  center1,center2;
 276     /* center of triangle of pts to define plane,
 277      * or center of vector if a vector is given
 278      */
 279   rvec  h1,h2,h3,h4,h5;         /* temp. vectors */
 280
 281   if (box)
 282     set_pbc(&pbc,ePBC,box);
 283
 284   switch(gnx1) {
 285   case 3:           /* group 1 defines plane */
 286     calculate_normal(index1,xzero,normal1,center1);
 287     break;
 288   case 2:           /* group 1 defines vector */
 289     rvec_sub(xzero[index1[0]],xzero[index1[1]],normal1);
 290     rvec_add(xzero[index1[0]],xzero[index1[1]],h1);
 291     svmul(0.5,h1,center1);  /* center is geometric mean */
 292     break;
 293   default:          /* group 1 does none of the above */
 294     gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
 295   }
 296
 297   switch(gnx2) {
 298   case 3:          /* group 2 defines plane */
 299     calculate_normal(index2,x,normal2,center2);
 300     break;
 301   case 2:          /* group 2 defines vector */
 302     rvec_sub(x[index2[0]],x[index2[1]],normal2);
 303     rvec_add(x[index2[0]],x[index2[1]],h2);
 304     svmul(0.5,h2,center2);  /* center is geometric mean */
 305     break;
 306   default:         /* group 2 does none of the above */
 307     gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
 308   }
 309
 310   *angle = cos_angle(normal1,normal2);
 311
 312   if (box)
 313     pbc_dx(&pbc,center1,center2,h3);
 314   else
 315     rvec_sub(center1,center2,h3);
 316   *distance = norm(h3);
 317
 318   if (gnx1 == 3 && gnx2 == 2) {
 319     if (box) {
 320       pbc_dx(&pbc,center1,x[index2[0]],h4);
 321       pbc_dx(&pbc,center1,x[index2[1]],h5);
 322     } else {
 323       rvec_sub(center1,x[index2[0]],h4);
 324       rvec_sub(center1,x[index2[1]],h5);
 325     }
 326     *distance1 = norm(h4);
 327     *distance2 = norm(h5);
 328   } else if (gnx1 == 2 && gnx2 ==3) {
 329     rvec_sub(center1,xzero[index1[0]],h4);
 330     rvec_sub(center1,xzero[index1[1]],h5);
 331     *distance1 = norm(h4);
 332     *distance2 = norm(h5);
 333   } else {
 334     *distance1 = 0; *distance2 = 0;
 335   }
 336 }
 337
 338
 339 void sgangle_plot_single(const char *fn,const char *afile,const char *dfile,
 340                          const char *d1file, const char *d2file,
 341                          atom_id index1[], int gnx1, char *grpn1,
 342                          atom_id index2[], int gnx2, char *grpn2,
 343                          t_topology *top,int ePBC, const output_env_t oenv)
 344 {
 345   FILE
 346     *sg_angle,           /* xvgr file with angles */
 347     *sg_distance = NULL, /* xvgr file with distances */
 348     *sg_distance1 = NULL,/* xvgr file with distance between plane and atom */
 349     *sg_distance2 = NULL;/* xvgr file with distance between plane and atom2 */
 350   real
 351     t,                   /* time */
 352     angle,               /* cosine of angle between two groups */
 353     distance,            /* distance between two groups. */
 354     distance1,           /* distance between plane and one of two atoms */
 355     distance2;           /* same for second of two atoms */
 356   int        status,natoms,teller=0;
 357   int        i;
 358   rvec       *x0;   /* coordinates, and coordinates corrected for pb */
 359   rvec       *xzero;
 360   matrix     box;
 361   char       buf[256];   /* for xvgr title */
 362
 363   if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
 364     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 365
 366   sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
 367   sg_angle = xvgropen(afile,buf,"Time (ps)","Cos(angle) ",oenv);
 368
 369   if (dfile) {
 370     sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
 371     sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
 372   }
 373
 374   if (d1file) {
 375     sprintf(buf,"Distance between plane and first atom of vector");
 376     sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
 377   }
 378
 379   if (d2file) {
 380     sprintf(buf,"Distance between plane and second atom of vector");
 381     sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
 382   }
 383
 384   snew(xzero,natoms);
 385
 386   do {
 387     teller++;
 388
 389     rm_pbc(&(top->idef),ePBC,natoms,box,x0,x0);
 390     if (teller==1) {
 391       for(i=0;i<natoms;i++)
 392         copy_rvec(x0[i],xzero[i]);
 393     }
 394
 395
 396     calc_angle_single(ePBC,box,
 397                       xzero,x0,index1,index2,gnx1,gnx2,&angle,
 398                       &distance,&distance1,&distance2);
 399
 400     fprintf(sg_angle,"%12g  %12g  %12g\n",t,angle,acos(angle)*180.0/M_PI);
 401     if (dfile)
 402       fprintf(sg_distance,"%12g  %12g\n",t,distance);
 403     if (d1file)
 404       fprintf(sg_distance1,"%12g  %12g\n",t,distance1);
 405     if (d2file)
 406       fprintf(sg_distance2,"%12g  %12g\n",t,distance1);
 407
 408   } while (read_next_x(oenv,status,&t,natoms,x0,box));
 409
 410   fprintf(stderr,"\n");
 411   close_trj(status);
 412   fclose(sg_angle);
 413   if (dfile)
 414     fclose(sg_distance);
 415   if (d1file)
 416     fclose(sg_distance1);
 417   if (d2file)
 418     fclose(sg_distance2);
 419
 420   sfree(x0);
 421 }
 422
 423
 424
 425 int gmx_sgangle(int argc,char *argv[])
 426 {
 427   const char *desc[] = {
 428     "Compute the angle and distance between two groups. ",
 429     "The groups are defined by a number of atoms given in an index file and",
 430     "may be two or three atoms in size.",
 431     "If -one is set, only one group should be specified in the index",
 432     "file and the angle between this group at time 0 and t will be computed.",
 433     "The angles calculated depend on the order in which the atoms are ",
 434     "given. Giving for instance 5 6 will rotate the vector 5-6 with ",
 435     "180 degrees compared to giving 6 5. [PAR]If three atoms are given, ",
 436     "the normal on the plane spanned by those three atoms will be",
 437     "calculated, using the formula  P1P2 x P1P3.",
 438     "The cos of the angle is calculated, using the inproduct of the two",
 439     "normalized vectors.[PAR]",
 440     "Here is what some of the file options do:[BR]",
 441     "-oa: 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]",
 442     "-od: Distance between two groups. Distance is taken from the center of one group to the center of the other group.[BR]",
 443     "-od1: If one plane and one vector is given, the distances for each of the atoms from the center of the plane is given seperately.[BR]",
 444     "-od2: For two planes this option has no meaning."
 445   };
 446
 447   output_env_t oenv;
 448   const char *fna, *fnd, *fnd1, *fnd2;
 449   char *      grpname[2];                       /* name of the two groups */
 450   int       gnx[2];                             /* size of the two groups */
 451   t_topology *top;                              /* topology             */
 452   int       ePBC;
 453   atom_id   *index[2];
 454   static bool bOne = FALSE, bZ=FALSE;
 455   t_pargs pa[] = {
 456     { "-one", FALSE, etBOOL, {&bOne},
 457       "Only one group compute angle between vector at time zero and time t"},
 458     { "-z", FALSE, etBOOL, {&bZ},
 459         "Use the Z-axis as reference" }
 460   };
 461 #define NPA asize(pa)
 462
 463   t_filenm  fnm[] = {                           /* files for g_sgangle  */
 464     { efTRX, "-f", NULL,  ffREAD },             /* trajectory file      */
 465     { efNDX, NULL, NULL,  ffREAD },             /* index file           */
 466     { efTPX, NULL, NULL,  ffREAD },             /* topology file        */
 467     { efXVG,"-oa","sg_angle",ffWRITE },         /* xvgr output file     */
 468     { efXVG, "-od","sg_dist",ffOPTWR },         /* xvgr output file     */
 469     { efXVG, "-od1", "sg_dist1",ffOPTWR },      /* xvgr output file     */
 470     { efXVG, "-od2", "sg_dist2",ffOPTWR }       /* xvgr output file     */
 471   };
 472
 473 #define NFILE asize(fnm)
 474
 475   CopyRight(stderr,argv[0]);
 476   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 477                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 478
 479
 480   top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);     /* read topology file */
 481
 482   fna = opt2fn("-oa",NFILE,fnm);
 483   fnd = opt2fn_null("-od",NFILE,fnm);
 484   fnd1 = opt2fn_null("-od1",NFILE,fnm);
 485   fnd2 = opt2fn_null("-od2",NFILE,fnm);
 486
 487   /* read index file. */
 488   if(bOne) {
 489     rd_index(ftp2fn(efNDX,NFILE,fnm),1,gnx,index,grpname);
 490     print_types(index[0],gnx[0],grpname[0],
 491                 index[0],gnx[0],grpname[0],top);
 492
 493     sgangle_plot_single(ftp2fn(efTRX,NFILE,fnm), fna, fnd, fnd1, fnd2,
 494                         index[0],gnx[0],grpname[0],
 495                         index[0],gnx[0],grpname[0],
 496                         top,ePBC,oenv);
 497   }  else {
 498     rd_index(ftp2fn(efNDX,NFILE,fnm),bZ ? 1 : 2,gnx,index,grpname);
 499     if (!bZ)
 500       print_types(index[0],gnx[0],grpname[0],
 501                 index[1],gnx[1],grpname[1],top);
 502     else {
 503       gnx[1] = 0;
 504       grpname[1] = strdup("Z-axis");
 505     }
 506     sgangle_plot(ftp2fn(efTRX,NFILE,fnm), fna, fnd, fnd1, fnd2,
 507                  index[0],gnx[0],grpname[0],
 508                  index[1],gnx[1],grpname[1],
 509                  top,ePBC,oenv);
 510   }
 511
 512   do_view(oenv,fna,"-nxy");     /* view xvgr file */
 513   do_view(oenv,fnd,"-nxy");     /* view xvgr file */
 514   do_view(oenv,fnd1,"-nxy");     /* view xvgr file */
 515   do_view(oenv,fnd2,"-nxy");     /* view xvgr file */
 516
 517   thanx(stderr);
 518   return 0;
 519 }