src/tools/gmx_angle.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41 #include <math.h>
  42 #include <string.h>
  43
  44 #include "sysstuff.h"
  45 #include "physics.h"
  46 #include "typedefs.h"
  47 #include "smalloc.h"
  48 #include "futil.h"
  49 #include "statutil.h"
  50 #include "copyrite.h"
  51 #include "vec.h"
  52 #include "index.h"
  53 #include "macros.h"
  54 #include "gmx_fatal.h"
  55 #include "xvgr.h"
  56 #include "gstat.h"
  57 #include "trnio.h"
  58 #include "gmx_ana.h"
  59
  60
  61 static void dump_dih_trn(int nframes,int nangles,real **dih,const char *fn,
  62                          real *time)
  63 {
  64   int    i,j,k,l,m,na;
  65   t_fileio *trn;
  66   rvec   *x;
  67   matrix box = {{2,0,0},{0,2,0},{0,0,2}};
  68
  69   na = (nangles*2);
  70   if ((na % 3) != 0)
  71     na = 1+na/3;
  72   else
  73     na = na/3;
  74   printf("There are %d dihedrals. Will fill %d atom positions with cos/sin\n",
  75          nangles,na);
  76   snew(x,na);
  77   trn = open_trn(fn,"w");
  78   for(i=0; (i<nframes); i++) {
  79     k = l = 0;
  80     for(j=0; (j<nangles); j++) {
  81       for(m=0; (m<2); m++) {
  82         x[k][l] = (m == 0) ? cos(dih[j][i]) : sin(dih[j][i]);
  83         l++;
  84         if (l == DIM) {
  85           l = 0;
  86           k++;
  87         }
  88       }
  89     }
  90     fwrite_trn(trn,i,time[i],0,box,na,x,NULL,NULL);
  91   }
  92   close_trn(trn);
  93   sfree(x);
  94 }
  95
  96 int gmx_g_angle(int argc,char *argv[])
  97 {
  98   static const char *desc[] = {
  99     "[TT]g_angle[tt] computes the angle distribution for a number of angles",
 100     "or dihedrals.[PAR]",
 101     "With option [TT]-ov[tt], you can plot the average angle of",
 102     "a group of angles as a function of time. With the [TT]-all[tt] option,",
 103     "the first graph is the average and the rest are the individual angles.[PAR]",
 104     "With the [TT]-of[tt] option, [TT]g_angle[tt] also calculates the fraction of trans",
 105     "dihedrals (only for dihedrals) as function of time, but this is",
 106     "probably only fun for a select few.[PAR]",
 107     "With option [TT]-oc[tt], a dihedral correlation function is calculated.[PAR]",
 108     "It should be noted that the index file must contain",
 109     "atom triplets for angles or atom quadruplets for dihedrals.",
 110     "If this is not the case, the program will crash.[PAR]",
 111     "With option [TT]-or[tt], a trajectory file is dumped containing cos and",
 112     "sin of selected dihedral angles, which subsequently can be used as",
 113     "input for a principal components analysis using [TT]g_covar[tt].[PAR]",
 114     "Option [TT]-ot[tt] plots when transitions occur between",
 115     "dihedral rotamers of multiplicity 3 and [TT]-oh[tt]",
 116     "records a histogram of the times between such transitions,",
 117     "assuming the input trajectory frames are equally spaced in time."
 118   };
 119   static const char *opt[] = { NULL, "angle", "dihedral", "improper", "ryckaert-bellemans", NULL };
 120   static gmx_bool bALL=FALSE,bChandler=FALSE,bAverCorr=FALSE,bPBC=TRUE;
 121   static real binwidth=1;
 122   t_pargs pa[] = {
 123     { "-type", FALSE, etENUM, {opt},
 124       "Type of angle to analyse" },
 125     { "-all",    FALSE,  etBOOL, {&bALL},
 126       "Plot all angles separately in the averages file, in the order of appearance in the index file." },
 127     { "-binwidth", FALSE, etREAL, {&binwidth},
 128       "binwidth (degrees) for calculating the distribution" },
 129     { "-periodic", FALSE, etBOOL, {&bPBC},
 130       "Print dihedral angles modulo 360 degrees" },
 131     { "-chandler", FALSE,  etBOOL, {&bChandler},
 132       "Use Chandler correlation function (N[trans] = 1, N[gauche] = 0) rather than cosine correlation function. Trans is defined as phi < -60 or phi > 60." },
 133     { "-avercorr", FALSE,  etBOOL, {&bAverCorr},
 134       "Average the correlation functions for the individual angles/dihedrals" }
 135   };
 136   static const char *bugs[] = {
 137     "Counting transitions only works for dihedrals with multiplicity 3"
 138   };
 139
 140   FILE       *out;
 141   real       tmp,dt;
 142   int        status,isize;
 143   atom_id    *index;
 144   char       *grpname;
 145   real       maxang,Jc,S2,norm_fac,maxstat;
 146   unsigned long mode;
 147   int        nframes,maxangstat,mult,*angstat;
 148   int        i,j,total,nangles,natoms,nat2,first,last,angind;
 149   gmx_bool       bAver,bRb,bPeriodic,
 150     bFrac,          /* calculate fraction too?  */
 151     bTrans,         /* worry about transtions too? */
 152     bCorr;          /* correlation function ? */
 153   real       t,aa,aver,aver2,aversig,fraction;       /* fraction trans dihedrals */
 154   double     tfrac=0;
 155   char       title[256];
 156   real       **dih=NULL;          /* mega array with all dih. angles at all times*/
 157   char       buf[80];
 158   real       *time,*trans_frac,*aver_angle;
 159   t_filenm   fnm[] = {
 160     { efTRX, "-f", NULL,  ffREAD  },
 161     { efNDX, NULL, "angle",  ffREAD  },
 162     { efXVG, "-od", "angdist",  ffWRITE },
 163     { efXVG, "-ov", "angaver",  ffOPTWR },
 164     { efXVG, "-of", "dihfrac",  ffOPTWR },
 165     { efXVG, "-ot", "dihtrans", ffOPTWR },
 166     { efXVG, "-oh", "trhisto",  ffOPTWR },
 167     { efXVG, "-oc", "dihcorr",  ffOPTWR },
 168     { efTRR, "-or", NULL,       ffOPTWR }
 169   };
 170 #define NFILE asize(fnm)
 171   int     npargs;
 172   t_pargs *ppa;
 173   output_env_t oenv;
 174
 175   CopyRight(stderr,argv[0]);
 176   npargs = asize(pa);
 177   ppa    = add_acf_pargs(&npargs,pa);
 178   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 179                     NFILE,fnm,npargs,ppa,asize(desc),desc,asize(bugs),bugs,
 180                     &oenv);
 181
 182   mult   = 4;
 183   maxang = 360.0;
 184   bRb    = FALSE;
 185   switch(opt[0][0]) {
 186   case 'a':
 187     mult   = 3;
 188     maxang = 180.0;
 189     break;
 190   case 'd':
 191     break;
 192   case 'i':
 193     break;
 194   case 'r':
 195     bRb = TRUE;
 196     break;
 197   }
 198
 199   if (opt2bSet("-or",NFILE,fnm)) {
 200     if (mult != 4)
 201       gmx_fatal(FARGS,"Can not combine angles with trn dump");
 202     else
 203       please_cite(stdout,"Mu2005a");
 204   }
 205
 206   /* Calculate bin size */
 207   maxangstat=(int)(maxang/binwidth+0.5);
 208   binwidth=maxang/maxangstat;
 209
 210   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 211   nangles=isize/mult;
 212   if ((isize % mult) != 0)
 213     gmx_fatal(FARGS,"number of index elements not multiple of %d, "
 214                 "these can not be %s\n",
 215                 mult,(mult==3) ? "angle triplets" : "dihedral quadruplets");
 216
 217
 218   /* Check whether specific analysis has to be performed */
 219   bCorr=opt2bSet("-oc",NFILE,fnm);
 220   bAver=opt2bSet("-ov",NFILE,fnm);
 221   bTrans=opt2bSet("-ot",NFILE,fnm);
 222   bFrac=opt2bSet("-of",NFILE,fnm);
 223   if (bTrans && opt[0][0] != 'd') {
 224     fprintf(stderr, "Option -ot should only accompany -type dihedral. Disabling -ot.\n");
 225     bTrans = FALSE;
 226   }
 227
 228   if (bChandler && !bCorr)
 229     bCorr=TRUE;
 230
 231   if (bFrac && !bRb) {
 232     fprintf(stderr,"Warning:"
 233             " calculating fractions as defined in this program\n"
 234             "makes sense for Ryckaert Bellemans dihs. only. Ignoring -of\n\n");
 235     bFrac = FALSE;
 236   }
 237
 238   if ( (bTrans || bFrac || bCorr) && mult==3)
 239     gmx_fatal(FARGS,"Can only do transition, fraction or correlation\n"
 240                 "on dihedrals. Select -d\n");
 241
 242   /*
 243    * We need to know the nr of frames so we can allocate memory for an array
 244    * with all dihedral angles at all timesteps. Works for me.
 245    */
 246   if (bTrans || bCorr  || bALL || opt2bSet("-or",NFILE,fnm))
 247     snew(dih,nangles);
 248
 249   snew(angstat,maxangstat);
 250
 251   read_ang_dih(ftp2fn(efTRX,NFILE,fnm),(mult == 3),
 252                bALL || bCorr || bTrans || opt2bSet("-or",NFILE,fnm),
 253                bRb,bPBC,maxangstat,angstat,
 254                &nframes,&time,isize,index,&trans_frac,&aver_angle,dih,
 255                oenv);
 256
 257   dt=(time[nframes-1]-time[0])/(nframes-1);
 258
 259   if (bAver) {
 260     sprintf(title,"Average Angle: %s",grpname);
 261     out=xvgropen(opt2fn("-ov",NFILE,fnm),
 262                  title,"Time (ps)","Angle (degrees)",oenv);
 263     for(i=0; (i<nframes); i++) {
 264       fprintf(out,"%10.5f  %8.3f",time[i],aver_angle[i]*RAD2DEG);
 265       if (bALL) {
 266         for(j=0; (j<nangles); j++)
 267           if (bPBC) {
 268             real dd = dih[j][i];
 269             fprintf(out,"  %8.3f",atan2(sin(dd),cos(dd))*RAD2DEG);
 270           }
 271           else
 272             fprintf(out,"  %8.3f",dih[j][i]*RAD2DEG);
 273       }
 274       fprintf(out,"\n");
 275     }
 276     ffclose(out);
 277   }
 278   if (opt2bSet("-or",NFILE,fnm))
 279     dump_dih_trn(nframes,nangles,dih,opt2fn("-or",NFILE,fnm),time);
 280
 281   if (bFrac) {
 282     sprintf(title,"Trans fraction: %s",grpname);
 283     out=xvgropen(opt2fn("-of",NFILE,fnm),
 284                   title,"Time (ps)","Fraction",oenv);
 285     tfrac = 0.0;
 286     for(i=0; (i<nframes); i++) {
 287       fprintf(out,"%10.5f  %10.3f\n",time[i],trans_frac[i]);
 288       tfrac += trans_frac[i];
 289     }
 290     ffclose(out);
 291
 292     tfrac/=nframes;
 293     fprintf(stderr,"Average trans fraction: %g\n",tfrac);
 294   }
 295   sfree(trans_frac);
 296
 297   if (bTrans)
 298     ana_dih_trans(opt2fn("-ot",NFILE,fnm),opt2fn("-oh",NFILE,fnm),
 299                   dih,nframes,nangles,grpname,time,bRb,oenv);
 300
 301   if (bCorr) {
 302     /* Autocorrelation function */
 303     if (nframes < 2)
 304       fprintf(stderr,"Not enough frames for correlation function\n");
 305     else {
 306
 307       if (bChandler) {
 308         real dval,sixty=DEG2RAD*60;
 309         gmx_bool bTest;
 310
 311         for(i=0; (i<nangles); i++)
 312           for(j=0; (j<nframes); j++) {
 313             dval = dih[i][j];
 314             if (bRb)
 315               bTest=(dval > -sixty) && (dval < sixty);
 316             else
 317               bTest=(dval < -sixty) || (dval > sixty);
 318             if (bTest)
 319               dih[i][j] = dval-tfrac;
 320             else
 321               dih[i][j] = -tfrac;
 322           }
 323       }
 324       if (bChandler)
 325         mode = eacNormal;
 326       else
 327         mode = eacCos;
 328       do_autocorr(opt2fn("-oc",NFILE,fnm), oenv,
 329                   "Dihedral Autocorrelation Function",
 330                   nframes,nangles,dih,dt,mode,bAverCorr);
 331     }
 332   }
 333
 334
 335   /* Determine the non-zero part of the distribution */
 336   for(first=0; (first < maxangstat-1) && (angstat[first+1] == 0); first++)
 337     ;
 338   for(last=maxangstat-1; (last > 0) && (angstat[last-1] == 0) ; last--)
 339     ;
 340
 341   aver=aver2=0;
 342   for(i=0; (i<nframes); i++) {
 343     aver  += RAD2DEG*aver_angle[i];
 344     aver2 += sqr(RAD2DEG*aver_angle[i]);
 345   }
 346   aver   /= (real) nframes;
 347   aver2  /= (real) nframes;
 348   aversig = sqrt(aver2-sqr(aver));
 349   printf("Found points in the range from %d to %d (max %d)\n",
 350          first,last,maxangstat);
 351   printf(" < angle >  = %g\n",aver);
 352   printf("< angle^2 > = %g\n",aver2);
 353   printf("Std. Dev.   = %g\n",aversig);
 354
 355   if (mult == 3)
 356     sprintf(title,"Angle Distribution: %s",grpname);
 357   else {
 358     sprintf(title,"Dihedral Distribution: %s",grpname);
 359
 360     calc_distribution_props(maxangstat,angstat,-180.0,0,NULL,&S2);
 361     fprintf(stderr,"Order parameter S^2 = %g\n",S2);
 362   }
 363
 364   bPeriodic=(mult==4) && (first==0) && (last==maxangstat-1);
 365
 366   out=xvgropen(opt2fn("-od",NFILE,fnm),title,"Degrees","",oenv);
 367   if (output_env_get_print_xvgr_codes(oenv))
 368     fprintf(out,"@    subtitle \"average angle: %g\\So\\N\"\n",aver);
 369   norm_fac=1.0/(nangles*nframes*binwidth);
 370   if (bPeriodic) {
 371     maxstat=0;
 372     for(i=first; (i<=last); i++)
 373       maxstat=max(maxstat,angstat[i]*norm_fac);
 374     fprintf(out,"@with g0\n");
 375     fprintf(out,"@    world xmin -180\n");
 376     fprintf(out,"@    world xmax  180\n");
 377     fprintf(out,"@    world ymin 0\n");
 378     fprintf(out,"@    world ymax %g\n",maxstat*1.05);
 379     fprintf(out,"@    xaxis  tick major 60\n");
 380     fprintf(out,"@    xaxis  tick minor 30\n");
 381     fprintf(out,"@    yaxis  tick major 0.005\n");
 382     fprintf(out,"@    yaxis  tick minor 0.0025\n");
 383   }
 384   for(i=first; (i<=last); i++)
 385     fprintf(out,"%10g  %10f\n",i*binwidth+180.0-maxang,angstat[i]*norm_fac);
 386   if ( bPeriodic )
 387     /* print first bin again as last one */
 388     fprintf(out,"%10g  %10f\n",180.0,angstat[0]*norm_fac);
 389
 390   ffclose(out);
 391
 392   do_view(oenv,opt2fn("-od",NFILE,fnm),"-nxy");
 393   if (bAver)
 394     do_view(oenv,opt2fn("-ov",NFILE,fnm),"-nxy");
 395
 396   thanx(stderr);
 397
 398   return 0;
 399 }