src/tools/gmx_angle.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 #include <string.h>
  40
  41 #include "sysstuff.h"
  42 #include "physics.h"
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "futil.h"
  46 #include "statutil.h"
  47 #include "copyrite.h"
  48 #include "vec.h"
  49 #include "index.h"
  50 #include "macros.h"
  51 #include "gmx_fatal.h"
  52 #include "xvgr.h"
  53 #include "gstat.h"
  54 #include "trnio.h"
  55 #include "gmx_ana.h"
  56
  57
  58 static void dump_dih_trn(int nframes,int nangles,real **dih,const char *fn,
  59                          real dt)
  60 {
  61   int    i,j,k,l,m,na;
  62   t_fileio *trn;
  63   rvec   *x;
  64   matrix box = {{2,0,0},{0,2,0},{0,0,2}};
  65
  66   na = (nangles*2);
  67   if ((na % 3) != 0)
  68     na = 1+na/3;
  69   else
  70     na = na/3;
  71   printf("There are %d dihedrals. Will fill %d atom positions with cos/sin\n",
  72          nangles,na);
  73   snew(x,na);
  74   trn = open_trn(fn,"w");
  75   for(i=0; (i<nframes); i++) {
  76     k = l = 0;
  77     for(j=0; (j<nangles); j++) {
  78       for(m=0; (m<2); m++) {
  79         x[k][l] = (m == 0) ? cos(dih[j][i]) : sin(dih[j][i]);
  80         l++;
  81         if (l == DIM) {
  82           l = 0;
  83           k++;
  84         }
  85       }
  86     }
  87     fwrite_trn(trn,i,(real)i*dt,0,box,na,x,NULL,NULL);
  88   }
  89   close_trn(trn);
  90   sfree(x);
  91 }
  92
  93 int gmx_g_angle(int argc,char *argv[])
  94 {
  95   static const char *desc[] = {
  96     "g_angle computes the angle distribution for a number of angles",
  97     "or dihedrals. This way you can check whether your simulation",
  98     "is correct. With option -ov you can plot the average angle of",
  99     "a group of angles as a function of time. With the -all option",
 100     "the first graph is the average, the rest are the individual angles.[PAR]",
 101     "With the -of option g_angle also calculates the fraction of trans",
 102     "dihedrals (only for dihedrals) as function of time, but this is",
 103     "probably only fun for a selected few.[PAR]",
 104     "With option -oc a dihedral correlation function is calculated.[PAR]",
 105     "It should be noted that the indexfile should contain",
 106     "atom-triples for angles or atom-quadruplets for dihedrals.",
 107     "If this is not the case, the program will crash.[PAR]",
 108     "With option [TT]-or[tt] a trajectory file is dumped containing cos and",
 109     "sin of selected dihedral angles which subsequently can be used as",
 110     "input for a PCA analysis using [TT]g_covar[tt]."
 111   };
 112   static const char *opt[] = { NULL, "angle", "dihedral", "improper", "ryckaert-bellemans", NULL };
 113   static gmx_bool bALL=FALSE,bChandler=FALSE,bAverCorr=FALSE,bPBC=TRUE;
 114   static real binwidth=1;
 115   t_pargs pa[] = {
 116     { "-type", FALSE, etENUM, {opt},
 117       "Type of angle to analyse" },
 118     { "-all",    FALSE,  etBOOL, {&bALL},
 119       "Plot all angles separately in the averages file, in the order of appearance in the index file." },
 120     { "-binwidth", FALSE, etREAL, {&binwidth},
 121       "binwidth (degrees) for calculating the distribution" },
 122     { "-periodic", FALSE, etBOOL, {&bPBC},
 123       "Print dihedral angles modulo 360 degrees" },
 124     { "-chandler", FALSE,  etBOOL, {&bChandler},
 125       "Use Chandler correlation function (N[trans] = 1, N[gauche] = 0) rather than cosine correlation function. Trans is defined as phi < -60 || phi > 60." },
 126     { "-avercorr", FALSE,  etBOOL, {&bAverCorr},
 127       "Average the correlation functions for the individual angles/dihedrals" }
 128   };
 129   static const char *bugs[] = {
 130     "Counting transitions only works for dihedrals with multiplicity 3"
 131   };
 132
 133   FILE       *out;
 134   real       tmp,dt;
 135   int        status,isize;
 136   atom_id    *index;
 137   char       *grpname;
 138   real       maxang,Jc,S2,norm_fac,maxstat;
 139   unsigned long mode;
 140   int        nframes,maxangstat,mult,*angstat;
 141   int        i,j,total,nangles,natoms,nat2,first,last,angind;
 142   gmx_bool       bAver,bRb,bPeriodic,
 143     bFrac,          /* calculate fraction too?  */
 144     bTrans,         /* worry about transtions too? */
 145     bCorr;          /* correlation function ? */
 146   real       t,aa,aver,aver2,aversig,fraction;       /* fraction trans dihedrals */
 147   double     tfrac=0;
 148   char       title[256];
 149   real       **dih=NULL;          /* mega array with all dih. angles at all times*/
 150   char       buf[80];
 151   real       *time,*trans_frac,*aver_angle;
 152   t_filenm   fnm[] = {
 153     { efTRX, "-f", NULL,  ffREAD  },
 154     { efNDX, NULL, "angle",  ffREAD  },
 155     { efXVG, "-od", "angdist",  ffWRITE },
 156     { efXVG, "-ov", "angaver",  ffOPTWR },
 157     { efXVG, "-of", "dihfrac",  ffOPTWR },
 158     { efXVG, "-ot", "dihtrans", ffOPTWR },
 159     { efXVG, "-oh", "trhisto",  ffOPTWR },
 160     { efXVG, "-oc", "dihcorr",  ffOPTWR },
 161     { efTRR, "-or", NULL,       ffOPTWR }
 162   };
 163 #define NFILE asize(fnm)
 164   int     npargs;
 165   t_pargs *ppa;
 166   output_env_t oenv;
 167
 168   CopyRight(stderr,argv[0]);
 169   npargs = asize(pa);
 170   ppa    = add_acf_pargs(&npargs,pa);
 171   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 172                     NFILE,fnm,npargs,ppa,asize(desc),desc,asize(bugs),bugs,
 173                     &oenv);
 174
 175   mult   = 4;
 176   maxang = 360.0;
 177   bRb    = FALSE;
 178   switch(opt[0][0]) {
 179   case 'a':
 180     mult   = 3;
 181     maxang = 180.0;
 182     break;
 183   case 'd':
 184     break;
 185   case 'i':
 186     break;
 187   case 'r':
 188     bRb = TRUE;
 189     break;
 190   }
 191
 192   if (opt2bSet("-or",NFILE,fnm)) {
 193     if (mult != 4)
 194       gmx_fatal(FARGS,"Can not combine angles with trn dump");
 195     else
 196       please_cite(stdout,"Mu2005a");
 197   }
 198
 199   /* Calculate bin size */
 200   maxangstat=(int)(maxang/binwidth+0.5);
 201   binwidth=maxang/maxangstat;
 202
 203   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 204   nangles=isize/mult;
 205   if ((isize % mult) != 0)
 206     gmx_fatal(FARGS,"number of index elements not multiple of %d, "
 207                 "these can not be %s\n",
 208                 mult,(mult==3) ? "angle triplets" : "dihedral quadruplets");
 209
 210
 211   /* Check whether specific analysis has to be performed */
 212   bCorr=opt2bSet("-oc",NFILE,fnm);
 213   bAver=opt2bSet("-ov",NFILE,fnm);
 214   bTrans=opt2bSet("-ot",NFILE,fnm);
 215   bFrac=opt2bSet("-of",NFILE,fnm);
 216
 217   if (bChandler && !bCorr)
 218     bCorr=TRUE;
 219
 220   if (bFrac && !bRb) {
 221     fprintf(stderr,"Warning:"
 222             " calculating fractions as defined in this program\n"
 223             "makes sense for Ryckaert Bellemans dihs. only. Ignoring -of\n\n");
 224     bFrac = FALSE;
 225   }
 226
 227   if ( (bTrans || bFrac || bCorr) && mult==3)
 228     gmx_fatal(FARGS,"Can only do transition, fraction or correlation\n"
 229                 "on dihedrals. Select -d\n");
 230
 231   /*
 232    * We need to know the nr of frames so we can allocate memory for an array
 233    * with all dihedral angles at all timesteps. Works for me.
 234    */
 235   if (bTrans || bCorr  || bALL || opt2bSet("-or",NFILE,fnm))
 236     snew(dih,nangles);
 237
 238   snew(angstat,maxangstat);
 239
 240   read_ang_dih(ftp2fn(efTRX,NFILE,fnm),(mult == 3),
 241                bALL || bCorr || bTrans || opt2bSet("-or",NFILE,fnm),
 242                bRb,bPBC,maxangstat,angstat,
 243                &nframes,&time,isize,index,&trans_frac,&aver_angle,dih,
 244                oenv);
 245
 246   dt=(time[nframes-1]-time[0])/(nframes-1);
 247
 248   if (bAver) {
 249     sprintf(title,"Average Angle: %s",grpname);
 250     out=xvgropen(opt2fn("-ov",NFILE,fnm),
 251                  title,"Time (ps)","Angle (degrees)",oenv);
 252     for(i=0; (i<nframes); i++) {
 253       fprintf(out,"%10.5f  %8.3f",time[i],aver_angle[i]*RAD2DEG);
 254       if (bALL) {
 255         for(j=0; (j<nangles); j++)
 256           if (bPBC) {
 257             real dd = dih[j][i];
 258             fprintf(out,"  %8.3f",atan2(sin(dd),cos(dd))*RAD2DEG);
 259           }
 260           else
 261             fprintf(out,"  %8.3f",dih[j][i]*RAD2DEG);
 262       }
 263       fprintf(out,"\n");
 264     }
 265     ffclose(out);
 266   }
 267   if (opt2bSet("-or",NFILE,fnm))
 268     dump_dih_trn(nframes,nangles,dih,opt2fn("-or",NFILE,fnm),dt);
 269
 270   if (bFrac) {
 271     sprintf(title,"Trans fraction: %s",grpname);
 272     out=xvgropen(opt2fn("-of",NFILE,fnm),
 273                   title,"Time (ps)","Fraction",oenv);
 274     tfrac = 0.0;
 275     for(i=0; (i<nframes); i++) {
 276       fprintf(out,"%10.5f  %10.3f\n",time[i],trans_frac[i]);
 277       tfrac += trans_frac[i];
 278     }
 279     ffclose(out);
 280
 281     tfrac/=nframes;
 282     fprintf(stderr,"Average trans fraction: %g\n",tfrac);
 283   }
 284   sfree(trans_frac);
 285
 286   if (bTrans)
 287     ana_dih_trans(opt2fn("-ot",NFILE,fnm),opt2fn("-oh",NFILE,fnm),
 288                   dih,nframes,nangles,grpname,time[0],dt,bRb,oenv);
 289
 290   if (bCorr) {
 291     /* Autocorrelation function */
 292     if (nframes < 2)
 293       fprintf(stderr,"Not enough frames for correlation function\n");
 294     else {
 295
 296       if (bChandler) {
 297         real dval,sixty=DEG2RAD*60;
 298         gmx_bool bTest;
 299
 300         for(i=0; (i<nangles); i++)
 301           for(j=0; (j<nframes); j++) {
 302             dval = dih[i][j];
 303             if (bRb)
 304               bTest=(dval > -sixty) && (dval < sixty);
 305             else
 306               bTest=(dval < -sixty) || (dval > sixty);
 307             if (bTest)
 308               dih[i][j] = dval-tfrac;
 309             else
 310               dih[i][j] = -tfrac;
 311           }
 312       }
 313       if (bChandler)
 314         mode = eacNormal;
 315       else
 316         mode = eacCos;
 317       do_autocorr(opt2fn("-oc",NFILE,fnm), oenv,
 318                   "Dihedral Autocorrelation Function",
 319                   nframes,nangles,dih,dt,mode,bAverCorr);
 320     }
 321   }
 322
 323
 324   /* Determine the non-zero part of the distribution */
 325   for(first=0; (first < maxangstat-1) && (angstat[first+1] == 0); first++)
 326     ;
 327   for(last=maxangstat-1; (last > 0) && (angstat[last-1] == 0) ; last--)
 328     ;
 329
 330   aver=aver2=0;
 331   for(i=0; (i<nframes); i++) {
 332     aver  += RAD2DEG*aver_angle[i];
 333     aver2 += sqr(RAD2DEG*aver_angle[i]);
 334   }
 335   aver   /= (real) nframes;
 336   aver2  /= (real) nframes;
 337   aversig = sqrt(aver2-sqr(aver));
 338   printf("Found points in the range from %d to %d (max %d)\n",
 339          first,last,maxangstat);
 340   printf(" < angle >  = %g\n",aver);
 341   printf("< angle^2 > = %g\n",aver2);
 342   printf("Std. Dev.   = %g\n",aversig);
 343
 344   if (mult == 3)
 345     sprintf(title,"Angle Distribution: %s",grpname);
 346   else {
 347     sprintf(title,"Dihedral Distribution: %s",grpname);
 348
 349     calc_distribution_props(maxangstat,angstat,-180.0,0,NULL,&S2);
 350     fprintf(stderr,"Order parameter S^2 = %g\n",S2);
 351   }
 352
 353   bPeriodic=(mult==4) && (first==0) && (last==maxangstat-1);
 354
 355   out=xvgropen(opt2fn("-od",NFILE,fnm),title,"Degrees","",oenv);
 356   if (output_env_get_print_xvgr_codes(oenv))
 357     fprintf(out,"@    subtitle \"average angle: %g\\So\\N\"\n",aver);
 358   norm_fac=1.0/(nangles*nframes*binwidth);
 359   if (bPeriodic) {
 360     maxstat=0;
 361     for(i=first; (i<=last); i++)
 362       maxstat=max(maxstat,angstat[i]*norm_fac);
 363     fprintf(out,"@with g0\n");
 364     fprintf(out,"@    world xmin -180\n");
 365     fprintf(out,"@    world xmax  180\n");
 366     fprintf(out,"@    world ymin 0\n");
 367     fprintf(out,"@    world ymax %g\n",maxstat*1.05);
 368     fprintf(out,"@    xaxis  tick major 60\n");
 369     fprintf(out,"@    xaxis  tick minor 30\n");
 370     fprintf(out,"@    yaxis  tick major 0.005\n");
 371     fprintf(out,"@    yaxis  tick minor 0.0025\n");
 372   }
 373   for(i=first; (i<=last); i++)
 374     fprintf(out,"%10g  %10f\n",i*binwidth+180.0-maxang,angstat[i]*norm_fac);
 375   if ( bPeriodic )
 376     /* print first bin again as last one */
 377     fprintf(out,"%10g  %10f\n",180.0,angstat[0]*norm_fac);
 378
 379   ffclose(out);
 380
 381   do_view(oenv,opt2fn("-od",NFILE,fnm),"-nxy");
 382   if (bAver)
 383     do_view(oenv,opt2fn("-ov",NFILE,fnm),"-nxy");
 384
 385   thanx(stderr);
 386
 387   return 0;
 388 }