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 #include "gmx_ana.h"
  54
  55
  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 */
  59
  60 static void print_types(atom_id index1[], int gnx1, char *group1,
  61                        atom_id index2[], int gnx2, char *group2,
  62                        t_topology *top)
  63 {
  64   int i,j;
  65
  66   fprintf(stderr,"\n");
  67   fprintf(stderr,"Group %s contains the following atoms: \n",group1);
  68   for(i=0;i<gnx1;i++)
  69     fprintf(stderr,"Atomname %d: %s\n",i,*(top->atoms.atomname[index1[i]]));
  70   fprintf(stderr,"\n");
  71
  72   fprintf(stderr,"Group %s contains the following atoms: \n",group2);
  73   for(j=0;j<gnx2;j++)
  74     fprintf(stderr,"Atomname %d: %s\n",j,*(top->atoms.atomname[index2[j]]));
  75   fprintf(stderr,"\n");
  76
  77   fprintf(stderr,"Careful: distance only makes sense in some situations.\n\n");
  78 }
  79
  80 static void calculate_normal(atom_id index[],rvec x[],rvec result,rvec center)
  81 {
  82   rvec c1,c2;
  83   int i;
  84
  85   /* calculate centroid of triangle spanned by the three points */
  86   for(i=0;i<3;i++)
  87     center[i] = (x[index[0]][i] + x[index[1]][i] + x[index[2]][i])/3;
  88
  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);
  92
  93   cprod(c1,c2,result);                    /* take crossproduct between these */
  94 }
  95
  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)
 101
 102 /* distance is distance between centers, distance 1 between center of plane
 103    and atom one of vector, distance 2 same for atom two
 104 */
 105
 106 {
 107   rvec
 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 */
 112   t_pbc pbc;
 113
 114   set_pbc(&pbc,ePBC,box);
 115
 116   switch(gnx1)
 117     {
 118     case 3:           /* group 1 defines plane */
 119       calculate_normal(index1,x,normal1,center1);
 120       break;
 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 */
 125       break;
 126     default:          /* group 1 does none of the above */
 127       gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
 128     }
 129
 130   switch(gnx2)
 131     {
 132     case 3:          /* group 2 defines plane */
 133       calculate_normal(index2,x,normal2,center2);
 134       break;
 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 */
 139       break;
 140     case 0:
 141       normal2[XX] = 0;
 142       normal2[YY] = 0;
 143       normal2[ZZ] = 1;
 144       center2[XX] = 0;
 145       center2[YY] = 0;
 146       center2[ZZ] = 0;
 147       break;
 148     default:         /* group 2 does none of the above */
 149       gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
 150     }
 151
 152   *angle = cos_angle(normal1,normal2);
 153
 154   if (box)
 155     pbc_dx(&pbc,center1,center2,h3);
 156   else
 157     rvec_sub(center1,center2,h3);
 158   *distance = norm(h3);
 159
 160   if (gnx1 == 3 && gnx2 == 2) {
 161     if (box) {
 162       pbc_dx(&pbc,center1,x[index2[0]],h4);
 163       pbc_dx(&pbc,center1,x[index2[1]],h5);
 164     } else {
 165       rvec_sub(center1,x[index2[0]],h4);
 166       rvec_sub(center1,x[index2[1]],h5);
 167     }
 168     *distance1 = norm(h4);
 169     *distance2 = norm(h5);
 170   }
 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);
 176   }
 177   else {
 178     *distance1 = 0; *distance2 = 0;
 179   }
 180 }
 181
 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)
 187 {
 188   FILE
 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 */
 193   real
 194     t,                   /* time */
 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 */
 199   t_trxstatus *status;
 200   int        natoms,teller=0;
 201   rvec       *x0;   /* coordinates, and coordinates corrected for pb */
 202   matrix     box;
 203   char       buf[256];   /* for xvgr title */
 204   gmx_rmpbc_t  gpbc=NULL;
 205
 206
 207   if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
 208     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 209
 210   sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
 211   sg_angle = xvgropen(afile,buf,"Time (ps)","Angle (degrees)",oenv);
 212
 213   if (dfile) {
 214     sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
 215     sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
 216   }
 217
 218   if (d1file) {
 219     sprintf(buf,"Distance between plane and first atom of vector");
 220     sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
 221   }
 222
 223   if (d2file) {
 224     sprintf(buf,"Distance between plane and second atom of vector");
 225     sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
 226   }
 227
 228   gpbc = gmx_rmpbc_init(&(top->idef),ePBC,natoms,box);
 229
 230   do
 231     {
 232       teller++;
 233
 234       gmx_rmpbc(gpbc,natoms,box,x0);
 235
 236       calc_angle(ePBC,box,x0,index1,index2,gnx1,gnx2,&angle,
 237                  &distance,&distance1,&distance2);
 238
 239       fprintf(sg_angle,"%12g  %12g  %12g\n",t,angle,acos(angle)*180.0/M_PI);
 240       if (dfile)
 241         fprintf(sg_distance,"%12g  %12g\n",t,distance);
 242       if (d1file)
 243         fprintf(sg_distance1,"%12g  %12g\n",t,distance1);
 244       if (d2file)
 245         fprintf(sg_distance2,"%12g  %12g\n",t,distance1);
 246
 247     } while (read_next_x(oenv,status,&t,natoms,x0,box));
 248
 249   gmx_rmpbc_done(gpbc);
 250
 251   fprintf(stderr,"\n");
 252   close_trj(status);
 253   ffclose(sg_angle);
 254   if (dfile)
 255     ffclose(sg_distance);
 256   if (d1file)
 257     ffclose(sg_distance1);
 258   if (d2file)
 259     ffclose(sg_distance2);
 260
 261   sfree(x0);
 262 }
 263
 264 static void calc_angle_single(int ePBC,
 265                               matrix box,
 266                               rvec xzero[],
 267                               rvec x[],
 268                               atom_id index1[],
 269                               atom_id index2[],
 270                               int gnx1,
 271                               int gnx2,
 272                               real *angle,
 273                               real *distance,
 274                               real *distance1,
 275                               real *distance2)
 276 {
 277   t_pbc pbc;
 278
 279   /* distance is distance between centers, distance 1 between center of plane
 280      and atom one of vector, distance 2 same for atom two
 281   */
 282
 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
 287      */
 288   rvec  h1,h2,h3,h4,h5;         /* temp. vectors */
 289
 290   if (box)
 291     set_pbc(&pbc,ePBC,box);
 292
 293   switch(gnx1) {
 294   case 3:           /* group 1 defines plane */
 295     calculate_normal(index1,xzero,normal1,center1);
 296     break;
 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 */
 301     break;
 302   default:          /* group 1 does none of the above */
 303     gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
 304   }
 305
 306   switch(gnx2) {
 307   case 3:          /* group 2 defines plane */
 308     calculate_normal(index2,x,normal2,center2);
 309     break;
 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 */
 314     break;
 315   default:         /* group 2 does none of the above */
 316     gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
 317   }
 318
 319   *angle = cos_angle(normal1,normal2);
 320
 321   if (box)
 322     pbc_dx(&pbc,center1,center2,h3);
 323   else
 324     rvec_sub(center1,center2,h3);
 325   *distance = norm(h3);
 326
 327   if (gnx1 == 3 && gnx2 == 2) {
 328     if (box) {
 329       pbc_dx(&pbc,center1,x[index2[0]],h4);
 330       pbc_dx(&pbc,center1,x[index2[1]],h5);
 331     } else {
 332       rvec_sub(center1,x[index2[0]],h4);
 333       rvec_sub(center1,x[index2[1]],h5);
 334     }
 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);
 342   } else {
 343     *distance1 = 0; *distance2 = 0;
 344   }
 345 }
 346
 347
 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)
 353 {
 354   FILE
 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 */
 359   real
 360     t,                   /* time */
 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 */
 365   t_trxstatus *status;
 366   int        natoms,teller=0;
 367   int        i;
 368   rvec       *x0;   /* coordinates, and coordinates corrected for pb */
 369   rvec       *xzero;
 370   matrix     box;
 371   char       buf[256];   /* for xvgr title */
 372   gmx_rmpbc_t  gpbc=NULL;
 373
 374
 375   if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
 376     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 377
 378   sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
 379   sg_angle = xvgropen(afile,buf,"Time (ps)","Cos(angle) ",oenv);
 380
 381   if (dfile) {
 382     sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
 383     sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
 384   }
 385
 386   if (d1file) {
 387     sprintf(buf,"Distance between plane and first atom of vector");
 388     sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
 389   }
 390
 391   if (d2file) {
 392     sprintf(buf,"Distance between plane and second atom of vector");
 393     sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
 394   }
 395
 396   snew(xzero,natoms);
 397   gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
 398
 399   do {
 400     teller++;
 401
 402     gmx_rmpbc(gpbc,natoms,box,x0);
 403     if (teller==1) {
 404       for(i=0;i<natoms;i++)
 405         copy_rvec(x0[i],xzero[i]);
 406     }
 407
 408
 409     calc_angle_single(ePBC,box,
 410                       xzero,x0,index1,index2,gnx1,gnx2,&angle,
 411                       &distance,&distance1,&distance2);
 412
 413     fprintf(sg_angle,"%12g  %12g  %12g\n",t,angle,acos(angle)*180.0/M_PI);
 414     if (dfile)
 415       fprintf(sg_distance,"%12g  %12g\n",t,distance);
 416     if (d1file)
 417       fprintf(sg_distance1,"%12g  %12g\n",t,distance1);
 418     if (d2file)
 419       fprintf(sg_distance2,"%12g  %12g\n",t,distance1);
 420
 421   } while (read_next_x(oenv,status,&t,natoms,x0,box));
 422   gmx_rmpbc_done(gpbc);
 423
 424   fprintf(stderr,"\n");
 425   close_trj(status);
 426   ffclose(sg_angle);
 427   if (dfile)
 428     ffclose(sg_distance);
 429   if (d1file)
 430     ffclose(sg_distance1);
 431   if (d2file)
 432     ffclose(sg_distance2);
 433
 434   sfree(x0);
 435 }
 436
 437
 438
 439 int gmx_sgangle(int argc,char *argv[])
 440 {
 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."
 459   };
 460
 461   output_env_t oenv;
 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             */
 466   int       ePBC;
 467   atom_id   *index[2];
 468   static gmx_bool bOne = FALSE, bZ=FALSE;
 469   t_pargs pa[] = {
 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 Z-axis as reference" }
 474   };
 475 #define NPA asize(pa)
 476
 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     */
 485   };
 486
 487 #define NFILE asize(fnm)
 488
 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);
 492
 493
 494   top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);     /* read topology file */
 495
 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);
 500
 501   /* read index file. */
 502   if(bOne) {
 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);
 506
 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],
 510                         top,ePBC,oenv);
 511   }  else {
 512     rd_index(ftp2fn(efNDX,NFILE,fnm),bZ ? 1 : 2,gnx,index,grpname);
 513     if (!bZ)
 514       print_types(index[0],gnx[0],grpname[0],
 515                 index[1],gnx[1],grpname[1],top);
 516     else {
 517       gnx[1] = 0;
 518       grpname[1] = strdup("Z-axis");
 519     }
 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],
 523                  top,ePBC,oenv);
 524   }
 525
 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 */
 530
 531   thanx(stderr);
 532   return 0;
 533 }