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