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. This way you can check whether your simulation",
 101     "is correct. 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, 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 selected few.[PAR]",
 107     "With option [TT]-oc[tt] a dihedral correlation function is calculated.[PAR]",
 108     "It should be noted that the index file should contain",
 109     "atom-triples 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 PCA 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 }