src/tools/g_angle.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 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gromacs Runs One Microsecond At Cannonball Speeds
  31  */
  32 static char *SRCID_g_angle_c = "$Id$";
  33 #include <math.h>
  34 #include <string.h>
  35 #include "sysstuff.h"
  36 #include "physics.h"
  37 #include "typedefs.h"
  38 #include "smalloc.h"
  39 #include "futil.h"
  40 #include "statutil.h"
  41 #include "copyrite.h"
  42 #include "vec.h"
  43 #include "rdgroup.h"
  44 #include "macros.h"
  45 #include "fatal.h"
  46 #include "xvgr.h"
  47 #include "gstat.h"
  48
  49 int main(int argc,char *argv[])
  50 {
  51   static char *desc[] = {
  52     "g_angle computes the angle distribution for a number of angles",
  53     "or dihedrals. This way you can check whether your simulation",
  54     "is correct. With option -ov you can plot the average angle of",
  55     "a group of angles as a function of time. With the -all option",
  56     "the first graph is the average, the rest are the individual angles.[PAR]",
  57     "With the -of option g_angle also calculates the fraction of trans",
  58     "dihedrals (only for dihedrals) as function of time, but this is",
  59     "probably only fun for a selected few.[PAR]",
  60     "With option -oc a dihedral correlation function is calculated.[PAR]",
  61     "It should be noted that the indexfile should contain",
  62     "atom-triples for angles or atom-quadruplets for dihedrals.",
  63     "If this is not the case, the program will crash."
  64   };
  65   static char *opt[] = { NULL, "angle", "dihedral", "improper", "ryckaert-bellemans", NULL };
  66   static bool bALL=FALSE,bChandler=FALSE,bAverCorr=FALSE;
  67   static real binwidth=1;
  68   t_pargs pa[] = {
  69     { "-type", FALSE, etENUM, {opt},
  70       "Type of angle to analyse" },
  71     { "-all",    FALSE,  etBOOL, {&bALL},
  72       "Plot all angles separately in the averages file, in the order of appearance in the index file." },
  73     { "-binwidth", FALSE, etREAL, {&binwidth},
  74       "binwidth (degrees) for calculating the distribution" },
  75     { "-chandler", FALSE,  etBOOL, {&bChandler},
  76       "Use Chandler correlation function (N[trans] = 1, N[gauche] = 0) rather than cosine correlation function. Trans is defined as phi < -60 || phi > 60." },
  77     { "-avercorr", FALSE,  etBOOL, {&bAverCorr},
  78       "Average the correlation functions for the individual angles/dihedrals" }
  79   };
  80   static char *bugs[] = {
  81     "Counting transitions only works for dihedrals with multiplicity 3"
  82   };
  83
  84   FILE       *out;
  85   real       tmp,dt;
  86   int        status,isize;
  87   atom_id    *index;
  88   char       *grpname;
  89   real       maxang,Jc,S2,norm_fac,maxstat;
  90   unsigned long mode;
  91   int        nframes,maxangstat,mult,*angstat;
  92   int        i,j,total,nangles,natoms,nat2,first,last,angind;
  93   bool       bAver,bRb,bPeriodic,
  94     bFrac,          /* calculate fraction too?  */
  95     bTrans,         /* worry about transtions too? */
  96     bCorr;          /* correlation function ? */
  97   real       t,aa,aver,aver2,aversig,fraction;       /* fraction trans dihedrals */
  98   double     tfrac=0;
  99   char       title[256];
 100   real       **dih=NULL;          /* mega array with all dih. angles at all times*/
 101   char       buf[80];
 102   real       *time,*trans_frac,*aver_angle;
 103   t_filenm   fnm[] = {
 104     { efTRX, "-f", NULL,  ffREAD  },
 105     { efTPX, NULL, NULL,  ffREAD  },
 106     { efNDX, NULL, "angle",  ffREAD  },
 107     { efXVG, "-od", "angdist",  ffWRITE },
 108     { efXVG, "-ov", "angaver",  ffOPTWR },
 109     { efXVG, "-of", "dihfrac",  ffOPTWR },
 110     { efXVG, "-ot", "dihtrans", ffOPTWR },
 111     { efXVG, "-oh", "trhisto",  ffOPTWR },
 112     { efXVG, "-oc", "dihcorr",  ffOPTWR }
 113   };
 114 #define NFILE asize(fnm)
 115   int     npargs;
 116   t_pargs *ppa;
 117
 118   CopyRight(stderr,argv[0]);
 119   npargs = asize(pa);
 120   ppa    = add_acf_pargs(&npargs,pa);
 121   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 122                     NFILE,fnm,npargs,ppa,asize(desc),desc,asize(bugs),bugs);
 123
 124   mult   = 4;
 125   maxang = 360.0;
 126   bRb    = FALSE;
 127   switch(opt[0][0]) {
 128   case 'a':
 129     mult   = 3;
 130     maxang = 180.0;
 131     break;
 132   case 'd':
 133     break;
 134   case 'i':
 135     break;
 136   case 'r':
 137     bRb = TRUE;
 138     break;
 139   }
 140
 141   /* Calculate bin size */
 142   maxangstat=(int)(maxang/binwidth+0.5);
 143   binwidth=maxang/maxangstat;
 144
 145   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 146   nangles=isize/mult;
 147   if ((isize % mult) != 0)
 148     fatal_error(0,"number of index elements not multiple of %d, "
 149                 "these can not be %s\n",
 150                 mult,(mult==3) ? "angle triplets" : "dihedral quadruplets");
 151
 152
 153   /* Check whether specific analysis has to be performed */
 154   bCorr=opt2bSet("-oc",NFILE,fnm);
 155   bAver=opt2bSet("-ov",NFILE,fnm);
 156   bTrans=opt2bSet("-ot",NFILE,fnm);
 157   bFrac=opt2bSet("-of",NFILE,fnm);
 158
 159   if (bChandler && !bCorr)
 160     bCorr=TRUE;
 161
 162   if (bFrac && !bRb) {
 163     fprintf(stderr,"Warning:"
 164             " calculating fractions as defined in this program\n"
 165             "makes sense for Ryckaert Bellemans dihs. only. Ignoring -of\n\n");
 166     bFrac = FALSE;
 167   }
 168
 169   if ( (bTrans || bFrac || bCorr) && mult==3)
 170     fatal_error(0,"Can only do transition, fraction or correlation\n"
 171                 "on dihedrals. Select -d\n");
 172
 173   /*
 174    * We need to know the nr of frames so we can allocate memory for an array
 175    * with all dihedral angles at all timesteps. Works for me.
 176    */
 177   if (bTrans || bCorr  || bALL)
 178     snew(dih,nangles);
 179
 180   snew(angstat,maxangstat);
 181
 182   read_ang_dih(ftp2fn(efTRX,NFILE,fnm),ftp2fn(efTPX,NFILE,fnm),(mult == 3),
 183                bALL || bCorr || bTrans,bRb,maxangstat,angstat,
 184                &nframes,&time,isize,index,&trans_frac,&aver_angle,dih);
 185
 186   dt=(time[nframes-1]-time[0])/(nframes-1);
 187
 188   if (bAver) {
 189     sprintf(title,"Average Angle: %s",grpname);
 190     out=xvgropen(opt2fn("-ov",NFILE,fnm),
 191                  title,"Time (ps)","Angle (degrees)");
 192     for(i=0; (i<nframes); i++) {
 193       fprintf(out,"%10.5f  %8.3f",time[i],aver_angle[i]*RAD2DEG);
 194       if (bALL)
 195         for(j=0; (j<nangles); j++)
 196           fprintf(out,"  %8.3f",dih[j][i]*RAD2DEG);
 197       fprintf(out,"\n");
 198     }
 199     fclose(out);
 200   }
 201
 202   if (bFrac) {
 203     sprintf(title,"Trans fraction: %s",grpname);
 204     out=xvgropen(opt2fn("-of",NFILE,fnm),
 205                   title,"Time (ps)","Fraction");
 206     tfrac = 0.0;
 207     for(i=0; (i<nframes); i++) {
 208       fprintf(out,"%10.5f  %10.3f\n",time[i],trans_frac[i]);
 209       tfrac += trans_frac[i];
 210     }
 211     fclose(out);
 212
 213     tfrac/=nframes;
 214     fprintf(stderr,"Average trans fraction: %g\n",tfrac);
 215   }
 216   sfree(trans_frac);
 217
 218   if (bTrans)
 219     ana_dih_trans(opt2fn("-ot",NFILE,fnm),opt2fn("-oh",NFILE,fnm),
 220                   dih,nframes,nangles,grpname,time[0],dt,bRb);
 221
 222   if (bCorr) {
 223     /* Autocorrelation function */
 224     if (nframes < 2)
 225       fprintf(stderr,"Not enough frames for correlation function\n");
 226     else {
 227
 228       if (bChandler) {
 229         real dval,sixty=DEG2RAD*60;
 230         bool bTest;
 231
 232         for(i=0; (i<nangles); i++)
 233           for(j=0; (j<nframes); j++) {
 234             dval = dih[i][j];
 235             if (bRb)
 236               bTest=(dval > -sixty) && (dval < sixty);
 237             else
 238               bTest=(dval < -sixty) || (dval > sixty);
 239             if (bTest)
 240               dih[i][j] = dval-tfrac;
 241             else
 242               dih[i][j] = -tfrac;
 243           }
 244       }
 245       if (bChandler)
 246         mode = eacNormal;
 247       else
 248         mode = eacCos;
 249       do_autocorr(opt2fn("-oc",NFILE,fnm),"Dihedral Autocorrelation Function",
 250                   nframes,nangles,dih,dt,mode,bAverCorr);
 251     }
 252   }
 253
 254
 255   /* Determine the non-zero part of the distribution */
 256   for(first=0; (first < maxangstat-1) && (angstat[first+1] == 0); first++)
 257     ;
 258   for(last=maxangstat-1; (last > 0) && (angstat[last-1] == 0) ; last--)
 259     ;
 260
 261   aver=aver2=0;
 262   for(i=0; (i<nframes); i++) {
 263     aver  += RAD2DEG*aver_angle[i];
 264     aver2 += sqr(RAD2DEG*aver_angle[i]);
 265   }
 266   aver   /= (real) nframes;
 267   aver2  /= (real) nframes;
 268   aversig = sqrt(aver2-sqr(aver));
 269   printf("Found points in the range from %d to %d (max %d)\n",
 270          first,last,maxangstat);
 271   printf(" < angle >  = %g\n",aver);
 272   printf("< angle^2 > = %g\n",aver2);
 273   printf("Std. Dev.   = %g\n",aversig);
 274
 275   if (mult == 3)
 276     sprintf(title,"Angle Distribution: %s",grpname);
 277   else {
 278     sprintf(title,"Dihedral Distribution: %s",grpname);
 279
 280     calc_distribution_props(maxangstat,angstat,-180.0,0,NULL,&S2);
 281     fprintf(stderr,"Order parameter S^2 = %g\n",S2);
 282   }
 283
 284   bPeriodic=(mult==4) && (first==0) && (last==maxangstat-1);
 285
 286   out=xvgropen(opt2fn("-od",NFILE,fnm),title,"Degrees","");
 287   fprintf(out,"@    subtitle \"average angle: %g\\So\\N\"\n",aver*RAD2DEG);
 288   norm_fac=1.0/(nangles*nframes*binwidth);
 289   if (bPeriodic) {
 290     maxstat=0;
 291     for(i=first; (i<=last); i++)
 292       maxstat=max(maxstat,angstat[i]*norm_fac);
 293     fprintf(out,"@with g0\n");
 294     fprintf(out,"@    world xmin -180\n");
 295     fprintf(out,"@    world xmax  180\n");
 296     fprintf(out,"@    world ymin 0\n");
 297     fprintf(out,"@    world ymax %g\n",maxstat*1.05);
 298     fprintf(out,"@    xaxis  tick major 60\n");
 299     fprintf(out,"@    xaxis  tick minor 30\n");
 300     fprintf(out,"@    yaxis  tick major 0.005\n");
 301     fprintf(out,"@    yaxis  tick minor 0.0025\n");
 302   }
 303   for(i=first; (i<=last); i++)
 304     fprintf(out,"%10g  %10f\n",i*binwidth+180.0-maxang,angstat[i]*norm_fac);
 305   if ( bPeriodic )
 306     /* print first bin again as last one */
 307     fprintf(out,"%10g  %10f\n",180.0,angstat[0]*norm_fac);
 308
 309   fclose(out);
 310
 311   do_view(opt2fn("-od",NFILE,fnm),NULL);
 312   if (bAver)
 313     do_view(opt2fn("-ov",NFILE,fnm),"-nxy");
 314
 315   thanx(stderr);
 316
 317   return 0;
 318 }