src/tools/gmx_dih.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
  40 #include "sysstuff.h"
  41 #include "string2.h"
  42 #include "copyrite.h"
  43 #include "futil.h"
  44 #include "smalloc.h"
  45 #include "statutil.h"
  46 #include "nrama.h"
  47 #include "physics.h"
  48 #include "macros.h"
  49 #include "xvgr.h"
  50 #include "vec.h"
  51 #include "gmx_ana.h"
  52
  53
  54 #define NOMIN 'X'
  55
  56 static void ana_dih(FILE *out,char *index,int nframes,real dih[],t_dih *dd)
  57 {
  58   int i;
  59   real mind,maxd,sum,av,var,prev,width;
  60   gmx_bool bTrans;
  61
  62   mind=5400,maxd=-5400,sum=0,av=0,var=0;
  63
  64   prev=dih[0];
  65   for(i=0; (i<nframes); i++) {
  66     if ((dih[i]-prev) > 180) {
  67       /* PBC.. */
  68       dih[i]-=360;
  69     }
  70     else if ((dih[i]-prev) < -180)
  71       dih[i]+=360;
  72     prev=dih[i];
  73
  74     sum+=dih[i];
  75     mind=min(mind,dih[i]);
  76     maxd=max(maxd,dih[i]);
  77   }
  78   av=sum/nframes;
  79   for(i=0; (i<nframes); i++)
  80     var+=sqr(dih[i]-av);
  81   var/=nframes;
  82   width=(360.0/dd->mult);
  83   bTrans=((maxd - mind) > width);
  84
  85   fprintf(out,"%-10s %10.3f %10.3f %10.3f %10.3f %10.3f %-10s%3.0f\n",
  86           index,mind,av,maxd,var,sqrt(var),
  87           bTrans ? "Yep" : "",width);
  88 }
  89
  90 static int find_min(real phi,int ntab,real phitab[])
  91 {
  92   int  i,imin;
  93   real mind,mm;
  94   real width;
  95
  96   /* Set closest minimum to the first one */
  97   width=360.0/ntab;
  98   mind=fabs(phi-phitab[0]);
  99   imin=0;
 100   for(i=1; (i<ntab); i++) {
 101     mm=fabs(phi-phitab[i]);
 102     if (mm < mind) {
 103       imin=i;
 104       mind=mm;
 105     }
 106   }
 107   if (mind < width*0.5 )
 108     return imin;
 109   else
 110     return -1;
 111 }
 112
 113 static int vphi(t_dih *dih,real phi,int mult)
 114 {
 115   static real m2[] = { 90, 270 };
 116   static real m3[] = { 60, 180, 300 };
 117   static real m4[] = { 45, 135, 225, 315 };
 118   static real m6[] = { 30, 90, 150, 210, 270, 330 };
 119
 120   real phiref;
 121   int  vpp=0;
 122
 123   phiref=RAD2DEG*(phi-dih->phi0);
 124   while (phiref < 0)
 125     phiref+=360;
 126   while (phiref > 360)
 127     phiref-=360;
 128
 129   switch(mult) {
 130   case 2:
 131     vpp=find_min(phiref,2,m2);
 132     break;
 133   case 3:
 134     vpp=find_min(phiref,3,m3);
 135     break;
 136   case 4:
 137     vpp=find_min(phiref,4,m4);
 138     break;
 139   case 6:
 140     vpp=find_min(phiref,6,m6);
 141     break;
 142   default:
 143     gmx_fatal(FARGS,"No such multiplicity %d",dih->mult);
 144   }
 145
 146   if (vpp == -1)
 147     return NOMIN;
 148   else
 149     return vpp+'0';
 150 }
 151
 152 typedef struct t_cluster {
 153   int    ndih;
 154   int    freq;
 155   char   *minimum;
 156   struct t_cluster *next;
 157 } t_cluster;
 158
 159 static t_cluster *search_cluster(t_cluster *cl,char *minimum)
 160 {
 161   t_cluster *ccl=cl;
 162
 163   while (ccl != NULL) {
 164     if (strcmp(minimum,ccl->minimum)==0)
 165       return ccl;
 166     ccl=ccl->next;
 167   }
 168   return NULL;
 169 }
 170
 171 static void add_cluster(t_cluster **cl,int ndih,char *minimum)
 172 {
 173   t_cluster *loper;
 174   t_cluster *ccl;
 175
 176   snew(ccl,1);
 177   ccl->ndih=ndih;
 178   ccl->freq=1;
 179   ccl->minimum=strdup(minimum);
 180   ccl->next=NULL;
 181
 182   if (*cl == NULL)
 183     *cl=ccl;
 184   else {
 185     loper=*cl;
 186     while (loper->next != NULL)
 187       loper=loper->next;
 188     loper->next=ccl;
 189   }
 190 }
 191
 192 static void p_cluster(FILE *out,t_cluster *cl)
 193 {
 194   t_cluster *loper;
 195
 196   fprintf(out,"* * * C L U S T E R   A N A L Y S I S * * *\n\n");
 197   fprintf(out," Frequency  Dihedral minima\n");
 198   loper=cl;
 199   while (loper != NULL) {
 200     fprintf(out,"%10d  %s\n",loper->freq,loper->minimum);
 201     loper=loper->next;
 202   }
 203 }
 204
 205 static void ana_cluster(FILE *out, t_xrama *xr,real **dih,real time[],
 206                         t_topology *top,int nframes,int mult)
 207 {
 208   t_cluster *cl=NULL,*scl;
 209   char      *minimum;
 210   int       i,j,nx;
 211
 212   /* Number of dihedrals + terminating NULL
 213    * this allows for using string routines
 214    */
 215   snew(minimum,xr->ndih+1);
 216
 217   for(i=0; (i<nframes); i++) {
 218     nx=0;
 219     for(j=0; (j<xr->ndih); j++) {
 220       minimum[j] = vphi(&xr->dih[j],dih[j][i],
 221                         mult == -1 ? xr->dih[j].mult : mult);
 222       if (minimum[j] == NOMIN)
 223         nx++;
 224     }
 225     if (nx == 0) {
 226       if ((scl=search_cluster(cl,minimum)) == NULL)
 227         add_cluster(&cl,xr->ndih,minimum);
 228       else
 229         scl->freq++;
 230     }
 231   }
 232   p_cluster(out,cl);
 233
 234   sfree(minimum);
 235 }
 236
 237 static void ana_trans(FILE *out, t_xrama *xr,real **dih,real time[],
 238                       t_topology *top,int nframes, const output_env_t oenv)
 239 {
 240   FILE *outd;
 241   real prev_phi,prev_psi;
 242   int  i,j,phi,psi;
 243   char buf[10];
 244
 245   fprintf(out,"\n\t* * * D I H E D R A L    S T A T I S T I C S * * *\n\n");
 246   fprintf(out,"%-10s %10s %10s %10s %10s %10s %10s\n",
 247           "index","minimum","average","maximum","variance","std.dev",
 248           "transition");
 249   for(i=0; (i<xr->ndih); i++) {
 250     sprintf(buf,"dih-%d",i);
 251     ana_dih(out,buf,nframes,dih[i],&(xr->dih[i]));
 252   }
 253   for(i=0; (i<xr->npp); i++) {
 254     sprintf(buf,"%s",xr->pp[i].label);
 255     outd=xvgropen(buf,"Dihedral Angles","Time (ps)","Degrees",oenv);
 256
 257     phi=xr->pp[i].iphi;
 258     psi=xr->pp[i].ipsi;
 259     prev_phi=dih[phi][0];
 260     prev_psi=dih[psi][0];
 261     for(j=0; (j<nframes); j++) {
 262       /* PBC.. */
 263       if ((dih[phi][j]-prev_phi) > 180)
 264         dih[phi][j]-=360;
 265       else if ((dih[phi][j]-prev_phi) < -180)
 266         dih[phi][j]+=360;
 267       prev_phi=dih[phi][j];
 268       if ((dih[psi][j]-prev_psi) > 180)
 269         dih[psi][j]-=360;
 270       else if ((dih[psi][j]-prev_psi) < -180)
 271         dih[psi][j]+=360;
 272       prev_psi=dih[psi][j];
 273       fprintf(outd,"%10g  %10g  %10g\n",time[j],prev_phi,prev_psi);
 274     }
 275     ffclose(outd);
 276   }
 277 }
 278
 279 int gmx_dih(int argc,char *argv[])
 280 {
 281   const char *desc[] = {
 282     "[TT]g_dih[tt] can do two things. The default is to analyze dihedral transitions",
 283     "by merely computing all the dihedral angles defined in your topology",
 284     "for the whole trajectory. When a dihedral flips over to another minimum",
 285     "an angle/time plot is made.[PAR]",
 286     "The opther option is to discretize the dihedral space into a number of",
 287     "bins, and group each conformation in dihedral space in the",
 288     "appropriate bin. The output is then given as a number of dihedral",
 289     "conformations sorted according to occupancy."
 290   };
 291   static int  mult = -1;
 292   static gmx_bool bSA  = FALSE;
 293   t_pargs pa[] = {
 294     { "-sa", FALSE, etBOOL, {&bSA},
 295       "Perform cluster analysis in dihedral space instead of analysing dihedral transitions." },
 296     { "-mult", FALSE, etINT, {&mult},
 297       "mulitiplicity for dihedral angles (by default read from topology)" }
 298   };
 299   FILE       *out;
 300   t_xrama    *xr;
 301   t_topology *top;
 302   real       **dih,*time;
 303   real       dd;
 304   int        i,nframes,maxframes=1000;
 305   output_env_t oenv;
 306   t_filenm   fnm[] = {
 307     { efTRX, "-f", NULL, ffREAD },
 308     { efTPX, NULL, NULL, ffREAD },
 309     { efOUT, NULL, NULL, ffWRITE }
 310   };
 311 #define NFILE asize(fnm)
 312
 313   CopyRight(stderr,argv[0]);
 314   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 315                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 316
 317   if (mult != -1)
 318     fprintf(stderr,"Using %d for dihedral multiplicity rather than topology values\n",mult);
 319
 320   snew(xr,1);
 321   init_rama(oenv,ftp2fn(efTRX,NFILE,fnm),
 322             ftp2fn(efTPX,NFILE,fnm),xr,3);
 323   top=read_top(ftp2fn(efTPX,NFILE,fnm),NULL);
 324
 325   /* Brute force malloc, may be too big... */
 326   snew(dih,xr->ndih);
 327   for(i=0; (i<xr->ndih); i++)
 328     snew(dih[i],maxframes);
 329   snew(time,maxframes);
 330
 331   fprintf(stderr,"\n");
 332   nframes = 0;
 333   while (new_data(xr)) {
 334     for(i=0; (i<xr->ndih); i++) {
 335       dd=xr->dih[i].ang*RAD2DEG;
 336       while (dd < 0)
 337         dd+=360;
 338       while (dd > 360)
 339         dd-=360;
 340       dih[i][nframes]=dd;
 341     }
 342     time[nframes]=xr->t;
 343     nframes++;
 344     if (nframes > maxframes) {
 345       maxframes += 1000;
 346       for(i=0; (i<xr->ndih); i++)
 347         srenew(dih[i],maxframes);
 348       srenew(time,maxframes);
 349     }
 350   }
 351
 352   fprintf(stderr,"\nCalculated all dihedrals, now analysing...\n");
 353
 354   out=ftp2FILE(efOUT,NFILE,fnm,"w");
 355
 356   if (bSA) {
 357     /* Cluster and structure analysis */
 358     ana_cluster(out,xr,dih,time,top,nframes,mult);
 359   }
 360   else {
 361     /* Analyse transitions... */
 362     ana_trans(out,xr,dih,time,top,nframes,oenv);
 363   }
 364   ffclose(out);
 365
 366   thanx(stderr);
 367
 368   return 0;
 369 }