src/tools/g_sgangle.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_g_sgangle_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "sysstuff.h"
  33 #include "string.h"
  34 #include "typedefs.h"
  35 #include "smalloc.h"
  36 #include "macros.h"
  37 #include "lutab.h"
  38 #include "rmpbc.h"
  39 #include "vec.h"
  40 #include "xvgr.h"
  41 #include "pbc.h"
  42 #include "copyrite.h"
  43 #include "futil.h"
  44 #include "statutil.h"
  45 #include "rdgroup.h"
  46
  47 /* this version only works correctly if one of the entries in the index file
  48    is a plane (three atoms specified) and the other a vector. Distance
  49    is calculated from the center of the plane to both atoms of the vector */
  50
  51 static void print_types(atom_id index1[], int gnx1, char *group1,
  52                        atom_id index2[], int gnx2, char *group2,
  53                        t_topology *top)
  54 {
  55   int i,j;
  56
  57   fprintf(stderr,"\n");
  58   fprintf(stderr,"Group %s contains the following atoms: \n",group1);
  59   for(i=0;i<gnx1;i++)
  60     fprintf(stderr,"Atomname %d: %s\n",i,*(top->atoms.atomname[index1[i]]));
  61   fprintf(stderr,"\n");
  62
  63   fprintf(stderr,"Group %s contains the following atoms: \n",group2);
  64   for(j=0;j<gnx2;j++)
  65     fprintf(stderr,"Atomname %d: %s\n",j,*(top->atoms.atomname[index2[j]]));
  66   fprintf(stderr,"\n");
  67
  68   fprintf(stderr,"Careful: distance only makes sense in some situations.\n\n");
  69 }
  70
  71 static void calculate_normal(atom_id index[],rvec x[],rvec result,rvec center)
  72 {
  73   rvec c1,c2;
  74   int i;
  75
  76   /* calculate centroid of triangle spanned by the three points */
  77   for(i=0;i<3;i++)
  78     center[i] = (x[index[0]][i] + x[index[1]][i] + x[index[2]][i])/3;
  79
  80   /* use P1P2 x P1P3 to calculate normal, given three points P1-P3 */
  81   rvec_sub(x[index[1]],x[index[0]],c1);    /* find two vectors */
  82   rvec_sub(x[index[2]],x[index[0]],c2);
  83
  84   oprod(c1,c2,result);                    /* take crossproduct between these */
  85 }
  86
  87 /* calculate the angle and distance between the two groups */
  88 static void calc_angle(matrix box,rvec x[], atom_id index1[],
  89                        atom_id index2[], int gnx1, int gnx2,
  90                        real *angle,      real *distance,
  91                        real *distance1,  real *distance2)
  92
  93 /* distance is distance between centers, distance 1 between center of plane
  94    and atom one of vector, distance 2 same for atom two
  95 */
  96
  97 {
  98   rvec
  99     normal1,normal2,    /* normals on planes of interest */
 100     center1,center2,    /* center of triangle of points given to define plane,*/
 101                         /* or center of vector if a vector is given */
 102     h1,h2,h3,h4,h5;     /* temp. vectors */
 103
 104   switch(gnx1)
 105     {
 106     case 3:           /* group 1 defines plane */
 107       calculate_normal(index1,x,normal1,center1);
 108       break;
 109     case 2:           /* group 1 defines vector */
 110       rvec_sub(x[index1[0]],x[index1[1]],normal1);
 111       rvec_add(x[index1[0]],x[index1[1]],h1);
 112       svmul(0.5,h1,center1);  /* center is geometric mean */
 113       break;
 114     default:          /* group 1 does none of the above */
 115       fatal_error(0,"Something wrong with contents of index file.\n");
 116     }
 117
 118   switch(gnx2)
 119     {
 120     case 3:          /* group 2 defines plane */
 121       calculate_normal(index2,x,normal2,center2);
 122       break;
 123     case 2:          /* group 2 defines vector */
 124       rvec_sub(x[index2[0]],x[index2[1]],normal2);
 125       rvec_add(x[index2[0]],x[index2[1]],h2);
 126       svmul(0.5,h2,center2);  /* center is geometric mean */
 127       break;
 128     default:         /* group 2 does none of the above */
 129       fatal_error(0,"Something wrong with contents of index file.\n");
 130     }
 131
 132   *angle = cos_angle(normal1,normal2);
 133
 134   rvec_sub(center1,center2,h3);
 135   *distance = norm(h3);
 136
 137   if (gnx1 == 3 && gnx2 == 2) {
 138     rvec_sub(center1,x[index2[0]],h4);
 139     rvec_sub(center1,x[index2[1]],h5);
 140     *distance1 = norm(h4);
 141     *distance2 = norm(h5);
 142   }
 143   else if (gnx1 == 2 && gnx2 ==3) {
 144     rvec_sub(center1,x[index1[0]],h4);
 145     rvec_sub(center1,x[index1[1]],h5);
 146     *distance1 = norm(h4);
 147     *distance2 = norm(h5);
 148   }
 149   else {
 150     *distance1 = 0; *distance2 = 0;
 151   }
 152
 153 }
 154
 155 void sgangle_plot(char *fn,char *afile,char *bfile,
 156                   char *cfile, char *dfile,
 157                   atom_id index1[], int gnx1, char *grpn1,
 158                   atom_id index2[], int gnx2, char *grpn2,
 159                   t_topology *top)
 160 {
 161   FILE
 162     *sg_angle,           /* xvgr file with angles */
 163     *sg_distance,        /* xvgr file with distances */
 164     *sg_distance1,       /* xvgr file with distance between plane and atom */
 165     *sg_distance2;       /* xvgr file with distance between plane and atom2 */
 166   real
 167     t,                   /* time */
 168     angle,               /* cosine of angle between two groups */
 169     distance,            /* distance between two groups. */
 170     distance1,           /* distance between plane and one of two atoms */
 171     distance2;           /* same for second of two atoms */
 172   int        status,natoms,teller=0;
 173   rvec       *x0;   /* coordinates, and coordinates corrected for pb */
 174   matrix     box;
 175   char       buf[256];   /* for xvgr title */
 176
 177   if ((natoms = read_first_x(&status,fn,&t,&x0,box)) == 0)
 178     fatal_error(0,"Could not read coordinates from statusfile\n");
 179
 180   sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
 181   sg_angle = xvgropen(afile,buf,"Time (ps)","Cos(angle) ");
 182
 183   sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
 184   sg_distance = xvgropen(bfile,buf,"Time (ps)","Distance (nm)");
 185
 186   sprintf(buf,"Distance between plane and first atom of vector");
 187   sg_distance1 = xvgropen(cfile,buf,"Time (ps)","Distance (nm)");
 188
 189   sprintf(buf,"Distance between plane and second atom of vector");
 190   sg_distance2 = xvgropen(dfile,buf,"Time (ps","Distance (nm)");
 191
 192   do
 193     {
 194       teller++;
 195
 196       rm_pbc(&(top->idef),top->atoms.nr,box,x0,x0);
 197
 198       calc_angle(box,x0,index1,index2,gnx1,gnx2,&angle,
 199                  &distance,&distance1,&distance2);
 200
 201       fprintf(sg_angle,"%12g  %12g  %12g\n",t,angle,acos(angle)*180.0/M_PI);
 202       fprintf(sg_distance,"%12g  %12g\n",t,distance);
 203       fprintf(sg_distance1,"%12g  %12g\n",t,distance1);
 204       fprintf(sg_distance2,"%12g  %12g\n",t,distance1);
 205
 206     } while (read_next_x(status,&t,natoms,x0,box));
 207
 208   fprintf(stderr,"\n");
 209   close_trj(status);
 210   fclose(sg_angle);
 211   fclose(sg_distance);
 212   fclose(sg_distance1);
 213   fclose(sg_distance2);
 214
 215   sfree(x0);
 216 }
 217
 218 int main(int argc,char *argv[])
 219 {
 220   static char *desc[] = {
 221     "Compute the angle and distance between two groups. ",
 222     "The groups are defined by a number of atoms given in an index file and",
 223     "may be two or three atoms in size.",
 224     "The angles calculated depend on the order in which the atoms are ",
 225     "given. Giving for instance 5 6 will rotate the vector 5-6 with ",
 226     "180 degrees compared to giving 6 5. [PAR]If three atoms are given, ",
 227     "the normal on the plane spanned by those three atoms will be",
 228     "calculated, using the formula  P1P2 x P1P3.",
 229     "The cos of the angle is calculated, using the inproduct of the two",
 230     "normalized vectors.[PAR]",
 231     "Here is what some of the file options do:[BR]",
 232     "-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]",
 233     "-od: Distance between two groups. Distance is taken from the center of one group to the center of the other group.[BR]",
 234     "-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]",
 235     "-od2: For two planes this option has no meaning."
 236   };
 237
 238   char      *grpname[2];                        /* name of the two groups */
 239   int       gnx[2];                             /* size of the two groups */
 240   t_topology *top;                              /* topology             */
 241   atom_id   *index[2];                          /* atom_id's of the atoms inthe groups */
 242   t_filenm  fnm[] = {                           /* files for g_sgangle  */
 243     { efTRX, "-f", NULL,  ffREAD },             /* trajectory file      */
 244     { efNDX, NULL, NULL,  ffREAD },             /* index file           */
 245     { efTPX, NULL, NULL,  ffREAD },             /* topology file        */
 246     { efXVG,"-oa","sg_angle",ffWRITE },         /* xvgr output file     */
 247     { efXVG, "-od","sg_dist",ffWRITE },         /* xvgr output file     */
 248     { efXVG, "-od1", "sg_dist1",ffWRITE },      /* xvgr output file     */
 249     { efXVG, "-od2", "sg_dist2",ffWRITE }       /* xvgr output file     */
 250   };
 251
 252 #define NFILE asize(fnm)
 253
 254   CopyRight(stderr,argv[0]);
 255   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
 256                     NFILE,fnm,0,NULL,asize(desc),desc,0,NULL);
 257
 258   init_lookup_table(stdout);
 259
 260   top = read_top(ftp2fn(efTPX,NFILE,fnm));     /* read topology file */
 261
 262   /* read index file. */
 263   rd_index(ftp2fn(efNDX,NFILE,fnm),2,gnx,index,grpname);
 264
 265   print_types(index[0],gnx[0],grpname[0],      /* show atomtypes, to check */
 266              index[1],gnx[1],grpname[1],top);  /* if index file is correct */
 267
 268
 269   sgangle_plot(ftp2fn(efTRX,NFILE,fnm),
 270                opt2fn("-oa",NFILE,fnm), /* make plot */
 271                opt2fn("-od",NFILE,fnm),
 272                opt2fn("-od1",NFILE,fnm),
 273                opt2fn("-od2",NFILE,fnm),
 274                index[0],gnx[0],grpname[0],
 275                index[1],gnx[1],grpname[1],
 276                top);
 277
 278   xvgr_file(opt2fn("-oa",NFILE,fnm),NULL);     /* view xvgr file */
 279   xvgr_file(opt2fn("-od",NFILE,fnm),NULL);     /* view xvgr file */
 280   xvgr_file(opt2fn("-od1",NFILE,fnm),NULL);
 281   xvgr_file(opt2fn("-od2",NFILE,fnm),NULL);
 282
 283   thanx(stdout);
 284   return 0;
 285 }
 286
 287
 288
 289
 290
 291
 292
 293
 294