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   int        status,natoms,teller=0;
 200   rvec       *x0;   /* coordinates, and coordinates corrected for pb */
 201   matrix     box;
 202   char       buf[256];   /* for xvgr title */
 203
 204   if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
 205     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 206
 207   sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
 208   sg_angle = xvgropen(afile,buf,"Time (ps)","Angle (degrees)",oenv);
 209
 210   if (dfile) {
 211     sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
 212     sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
 213   }
 214
 215   if (d1file) {
 216     sprintf(buf,"Distance between plane and first atom of vector");
 217     sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
 218   }
 219
 220   if (d2file) {
 221     sprintf(buf,"Distance between plane and second atom of vector");
 222     sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
 223   }
 224
 225   do
 226     {
 227       teller++;
 228
 229       rm_pbc(&(top->idef),ePBC,natoms,box,x0,x0);
 230
 231       calc_angle(ePBC,box,x0,index1,index2,gnx1,gnx2,&angle,
 232                  &distance,&distance1,&distance2);
 233
 234       fprintf(sg_angle,"%12g  %12g  %12g\n",t,angle,acos(angle)*180.0/M_PI);
 235       if (dfile)
 236         fprintf(sg_distance,"%12g  %12g\n",t,distance);
 237       if (d1file)
 238         fprintf(sg_distance1,"%12g  %12g\n",t,distance1);
 239       if (d2file)
 240         fprintf(sg_distance2,"%12g  %12g\n",t,distance1);
 241
 242     } while (read_next_x(oenv,status,&t,natoms,x0,box));
 243
 244   fprintf(stderr,"\n");
 245   close_trj(status);
 246   ffclose(sg_angle);
 247   if (dfile)
 248     ffclose(sg_distance);
 249   if (d1file)
 250     ffclose(sg_distance1);
 251   if (d2file)
 252     ffclose(sg_distance2);
 253
 254   sfree(x0);
 255 }
 256
 257 static void calc_angle_single(int ePBC,
 258                               matrix box,
 259                               rvec xzero[],
 260                               rvec x[],
 261                               atom_id index1[],
 262                               atom_id index2[],
 263                               int gnx1,
 264                               int gnx2,
 265                               real *angle,
 266                               real *distance,
 267                               real *distance1,
 268                               real *distance2)
 269 {
 270   t_pbc pbc;
 271
 272   /* distance is distance between centers, distance 1 between center of plane
 273      and atom one of vector, distance 2 same for atom two
 274   */
 275
 276   rvec  normal1,normal2;        /* normals on planes of interest */
 277   rvec  center1,center2;
 278     /* center of triangle of pts to define plane,
 279      * or center of vector if a vector is given
 280      */
 281   rvec  h1,h2,h3,h4,h5;         /* temp. vectors */
 282
 283   if (box)
 284     set_pbc(&pbc,ePBC,box);
 285
 286   switch(gnx1) {
 287   case 3:           /* group 1 defines plane */
 288     calculate_normal(index1,xzero,normal1,center1);
 289     break;
 290   case 2:           /* group 1 defines vector */
 291     rvec_sub(xzero[index1[0]],xzero[index1[1]],normal1);
 292     rvec_add(xzero[index1[0]],xzero[index1[1]],h1);
 293     svmul(0.5,h1,center1);  /* center is geometric mean */
 294     break;
 295   default:          /* group 1 does none of the above */
 296     gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
 297   }
 298
 299   switch(gnx2) {
 300   case 3:          /* group 2 defines plane */
 301     calculate_normal(index2,x,normal2,center2);
 302     break;
 303   case 2:          /* group 2 defines vector */
 304     rvec_sub(x[index2[0]],x[index2[1]],normal2);
 305     rvec_add(x[index2[0]],x[index2[1]],h2);
 306     svmul(0.5,h2,center2);  /* center is geometric mean */
 307     break;
 308   default:         /* group 2 does none of the above */
 309     gmx_fatal(FARGS,"Something wrong with contents of index file.\n");
 310   }
 311
 312   *angle = cos_angle(normal1,normal2);
 313
 314   if (box)
 315     pbc_dx(&pbc,center1,center2,h3);
 316   else
 317     rvec_sub(center1,center2,h3);
 318   *distance = norm(h3);
 319
 320   if (gnx1 == 3 && gnx2 == 2) {
 321     if (box) {
 322       pbc_dx(&pbc,center1,x[index2[0]],h4);
 323       pbc_dx(&pbc,center1,x[index2[1]],h5);
 324     } else {
 325       rvec_sub(center1,x[index2[0]],h4);
 326       rvec_sub(center1,x[index2[1]],h5);
 327     }
 328     *distance1 = norm(h4);
 329     *distance2 = norm(h5);
 330   } else if (gnx1 == 2 && gnx2 ==3) {
 331     rvec_sub(center1,xzero[index1[0]],h4);
 332     rvec_sub(center1,xzero[index1[1]],h5);
 333     *distance1 = norm(h4);
 334     *distance2 = norm(h5);
 335   } else {
 336     *distance1 = 0; *distance2 = 0;
 337   }
 338 }
 339
 340
 341 void sgangle_plot_single(const char *fn,const char *afile,const char *dfile,
 342                          const char *d1file, const char *d2file,
 343                          atom_id index1[], int gnx1, char *grpn1,
 344                          atom_id index2[], int gnx2, char *grpn2,
 345                          t_topology *top,int ePBC, const output_env_t oenv)
 346 {
 347   FILE
 348     *sg_angle,           /* xvgr file with angles */
 349     *sg_distance = NULL, /* xvgr file with distances */
 350     *sg_distance1 = NULL,/* xvgr file with distance between plane and atom */
 351     *sg_distance2 = NULL;/* xvgr file with distance between plane and atom2 */
 352   real
 353     t,                   /* time */
 354     angle,               /* cosine of angle between two groups */
 355     distance,            /* distance between two groups. */
 356     distance1,           /* distance between plane and one of two atoms */
 357     distance2;           /* same for second of two atoms */
 358   int        status,natoms,teller=0;
 359   int        i;
 360   rvec       *x0;   /* coordinates, and coordinates corrected for pb */
 361   rvec       *xzero;
 362   matrix     box;
 363   char       buf[256];   /* for xvgr title */
 364
 365   if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
 366     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 367
 368   sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
 369   sg_angle = xvgropen(afile,buf,"Time (ps)","Cos(angle) ",oenv);
 370
 371   if (dfile) {
 372     sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
 373     sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
 374   }
 375
 376   if (d1file) {
 377     sprintf(buf,"Distance between plane and first atom of vector");
 378     sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
 379   }
 380
 381   if (d2file) {
 382     sprintf(buf,"Distance between plane and second atom of vector");
 383     sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
 384   }
 385
 386   snew(xzero,natoms);
 387
 388   do {
 389     teller++;
 390
 391     rm_pbc(&(top->idef),ePBC,natoms,box,x0,x0);
 392     if (teller==1) {
 393       for(i=0;i<natoms;i++)
 394         copy_rvec(x0[i],xzero[i]);
 395     }
 396
 397
 398     calc_angle_single(ePBC,box,
 399                       xzero,x0,index1,index2,gnx1,gnx2,&angle,
 400                       &distance,&distance1,&distance2);
 401
 402     fprintf(sg_angle,"%12g  %12g  %12g\n",t,angle,acos(angle)*180.0/M_PI);
 403     if (dfile)
 404       fprintf(sg_distance,"%12g  %12g\n",t,distance);
 405     if (d1file)
 406       fprintf(sg_distance1,"%12g  %12g\n",t,distance1);
 407     if (d2file)
 408       fprintf(sg_distance2,"%12g  %12g\n",t,distance1);
 409
 410   } while (read_next_x(oenv,status,&t,natoms,x0,box));
 411
 412   fprintf(stderr,"\n");
 413   close_trj(status);
 414   ffclose(sg_angle);
 415   if (dfile)
 416     ffclose(sg_distance);
 417   if (d1file)
 418     ffclose(sg_distance1);
 419   if (d2file)
 420     ffclose(sg_distance2);
 421
 422   sfree(x0);
 423 }
 424
 425
 426
 427 int gmx_sgangle(int argc,char *argv[])
 428 {
 429   const char *desc[] = {
 430     "Compute the angle and distance between two groups. ",
 431     "The groups are defined by a number of atoms given in an index file and",
 432     "may be two or three atoms in size.",
 433     "If -one is set, only one group should be specified in the index",
 434     "file and the angle between this group at time 0 and t will be computed.",
 435     "The angles calculated depend on the order in which the atoms are ",
 436     "given. Giving for instance 5 6 will rotate the vector 5-6 with ",
 437     "180 degrees compared to giving 6 5. [PAR]If three atoms are given, ",
 438     "the normal on the plane spanned by those three atoms will be",
 439     "calculated, using the formula  P1P2 x P1P3.",
 440     "The cos of the angle is calculated, using the inproduct of the two",
 441     "normalized vectors.[PAR]",
 442     "Here is what some of the file options do:[BR]",
 443     "-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]",
 444     "-od: Distance between two groups. Distance is taken from the center of one group to the center of the other group.[BR]",
 445     "-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]",
 446     "-od2: For two planes this option has no meaning."
 447   };
 448
 449   output_env_t oenv;
 450   const char *fna, *fnd, *fnd1, *fnd2;
 451   char *      grpname[2];                       /* name of the two groups */
 452   int       gnx[2];                             /* size of the two groups */
 453   t_topology *top;                              /* topology             */
 454   int       ePBC;
 455   atom_id   *index[2];
 456   static bool bOne = FALSE, bZ=FALSE;
 457   t_pargs pa[] = {
 458     { "-one", FALSE, etBOOL, {&bOne},
 459       "Only one group compute angle between vector at time zero and time t"},
 460     { "-z", FALSE, etBOOL, {&bZ},
 461         "Use the Z-axis as reference" }
 462   };
 463 #define NPA asize(pa)
 464
 465   t_filenm  fnm[] = {                           /* files for g_sgangle  */
 466     { efTRX, "-f", NULL,  ffREAD },             /* trajectory file      */
 467     { efNDX, NULL, NULL,  ffREAD },             /* index file           */
 468     { efTPX, NULL, NULL,  ffREAD },             /* topology file        */
 469     { efXVG,"-oa","sg_angle",ffWRITE },         /* xvgr output file     */
 470     { efXVG, "-od","sg_dist",ffOPTWR },         /* xvgr output file     */
 471     { efXVG, "-od1", "sg_dist1",ffOPTWR },      /* xvgr output file     */
 472     { efXVG, "-od2", "sg_dist2",ffOPTWR }       /* xvgr output file     */
 473   };
 474
 475 #define NFILE asize(fnm)
 476
 477   CopyRight(stderr,argv[0]);
 478   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 479                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 480
 481
 482   top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);     /* read topology file */
 483
 484   fna = opt2fn("-oa",NFILE,fnm);
 485   fnd = opt2fn_null("-od",NFILE,fnm);
 486   fnd1 = opt2fn_null("-od1",NFILE,fnm);
 487   fnd2 = opt2fn_null("-od2",NFILE,fnm);
 488
 489   /* read index file. */
 490   if(bOne) {
 491     rd_index(ftp2fn(efNDX,NFILE,fnm),1,gnx,index,grpname);
 492     print_types(index[0],gnx[0],grpname[0],
 493                 index[0],gnx[0],grpname[0],top);
 494
 495     sgangle_plot_single(ftp2fn(efTRX,NFILE,fnm), fna, fnd, fnd1, fnd2,
 496                         index[0],gnx[0],grpname[0],
 497                         index[0],gnx[0],grpname[0],
 498                         top,ePBC,oenv);
 499   }  else {
 500     rd_index(ftp2fn(efNDX,NFILE,fnm),bZ ? 1 : 2,gnx,index,grpname);
 501     if (!bZ)
 502       print_types(index[0],gnx[0],grpname[0],
 503                 index[1],gnx[1],grpname[1],top);
 504     else {
 505       gnx[1] = 0;
 506       grpname[1] = strdup("Z-axis");
 507     }
 508     sgangle_plot(ftp2fn(efTRX,NFILE,fnm), fna, fnd, fnd1, fnd2,
 509                  index[0],gnx[0],grpname[0],
 510                  index[1],gnx[1],grpname[1],
 511                  top,ePBC,oenv);
 512   }
 513
 514   do_view(oenv,fna,"-nxy");     /* view xvgr file */
 515   do_view(oenv,fnd,"-nxy");     /* view xvgr file */
 516   do_view(oenv,fnd1,"-nxy");     /* view xvgr file */
 517   do_view(oenv,fnd2,"-nxy");     /* view xvgr file */
 518
 519   thanx(stderr);
 520   return 0;
 521 }