src/tools/gmx_dih.c

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