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