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