src/tools/g_helix.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_helix_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "confio.h"
  33 #include "copyrite.h"
  34 #include "fatal.h"
  35 #include "fitahx.h"
  36 #include "futil.h"
  37 #include "gstat.h"
  38 #include "wgms.h"
  39 #include "hxprops.h"
  40 #include "macros.h"
  41 #include "maths.h"
  42 #include "pbc.h"
  43 #include "tpxio.h"
  44 #include "physics.h"
  45 #include "rdgroup.h"
  46 #include "smalloc.h"
  47 #include "statutil.h"
  48 #include "string.h"
  49 #include "sysstuff.h"
  50 #include "txtdump.h"
  51 #include "typedefs.h"
  52 #include "vec.h"
  53 #include "xvgr.h"
  54
  55 void dump_ahx(int nres,
  56               t_bb bb[],rvec x[],matrix box,int teller)
  57 {
  58   FILE *fp;
  59   char buf[256];
  60   int  i;
  61
  62   sprintf(buf,"dump%d.gro",teller);
  63   fp=ffopen(buf,"w");
  64   fprintf(fp,"Dumping fitted helix frame %d\n",teller);
  65   fprintf(fp,"%5d\n",nres*5);
  66   for(i=0; (i<nres); i++) {
  67 #define PR(AA) fprintf(fp,"%5d%5s%5s%5d%8.3f%8.3f%8.3f\n",i+1,"GLY",#AA,bb[i].AA,x[bb[i].AA][XX],x[bb[i].AA][YY],x[bb[i].AA][ZZ]); fflush(fp)
  68     if (bb[i].bHelix) {
  69       PR(N);
  70       PR(H);
  71       PR(CA);
  72       PR(C);
  73       PR(O);
  74     }
  75   }
  76   for(i=0; (i<DIM); i++)
  77     fprintf(fp,"%10.5f",box[i][i]);
  78   fprintf(fp,"\n");
  79   fclose(fp);
  80 }
  81
  82 int get_prop(char prop[])
  83 {
  84   static char *ppp[efhNR] = {
  85     "RAD", "TWIST", "RISE", "LEN", "NHX", "DIP", "RMS", "CPHI",
  86     "RMSA", "PHI", "PSI", "HB3", "HB4", "HB5", "CD222"
  87   };
  88   int i,n;
  89
  90   fprintf(stderr,"Select something\n");
  91   for(i=0; (i<efhNR); ) {
  92     fprintf(stderr,"%2d: %10s",i,ppp[i]); i++;
  93     if (i<efhNR) {
  94       fprintf(stderr,"      %2d: %10s\n",i,ppp[i]);
  95       i++;
  96     }
  97     else
  98       fprintf(stderr,"\n");
  99   }
 100   do {
 101     scanf("%d",&n);
 102   } while ((n < 0) || (n >= efhNR));
 103
 104   strcpy(prop,ppp[n]);
 105
 106   return n;
 107 }
 108
 109
 110 void dump_otrj(FILE *otrj,int natoms,atom_id all_index[],rvec x[],
 111                real fac,rvec xav[])
 112 {
 113   static FILE *fp=NULL;
 114   static real fac0=1.0;
 115
 116   int  i,ai,m;
 117   real xa,xm;
 118
 119   if (fp==NULL) {
 120     fp=ffopen("WEDGAMMA10.DAT","w");
 121     fac0=fac;
 122   }
 123   fac/=fac0;
 124
 125   fprintf(fp,"%10g\n",fac);
 126
 127   for(i=0; (i<natoms); i++) {
 128     ai=all_index[i];
 129     for(m=0; (m<DIM); m++) {
 130       xa = xav[ai][m];
 131       xm = x[ai][m]*fac;
 132       xav[ai][m] = xa+xm;
 133       x[ai][m]   = xm;
 134     }
 135   }
 136   write_gms_ndx(otrj,natoms,all_index,x,NULL);
 137 }
 138
 139 int main(int argc,char *argv[])
 140 {
 141   static char *desc[] = {
 142     "g_helix computes all kind of helix properties. First, the peptide",
 143     "is checked to find the longest helical part. This is determined by",
 144     "Hydrogen bonds and Phi/Psi angles.",
 145     "That bit is fitted",
 146     "to an ideal helix around the Z-axis and centered around the origin.",
 147     "Then the following properties are computed:[PAR]",
 148     "[BB]1.[bb] Helix radius (file radius.xvg). This is merely the",
 149     "RMS deviation in two dimensions for all Calpha atoms.",
 150     "it is calced as sqrt((SUM i(x^2(i)+y^2(i)))/N), where N is the number",
 151     "of backbone atoms. For an ideal helix the radius is 0.23 nm[BR]",
 152     "[BB]2.[bb] Twist (file twist.xvg). The average helical angle per",
 153     "residue is calculated. For alpha helix it is 100 degrees,",
 154     "for 3-10 helices it will be smaller,",
 155     "for 5-helices it will be larger.[BR]",
 156     "[BB]3.[bb] Rise per residue (file rise.xvg). The helical rise per",
 157     "residue is plotted as the difference in Z-coordinate between Ca",
 158     "atoms. For an ideal helix this is 0.15 nm[BR]",
 159     "[BB]4.[bb] Total helix length (file len-ahx.xvg). The total length",
 160     "of the",
 161     "helix in nm. This is simply the average rise (see above) times the",
 162     "number of helical residues (see below).[BR]",
 163     "[BB]5.[bb] Number of helical residues (file n-ahx.xvg). The title says",
 164     "it all.[BR]",
 165     "[BB]6.[bb] Helix Dipole, backbone only (file dip-ahx.xvg).[BR]",
 166     "[BB]7.[bb] RMS deviation from ideal helix, calculated for the Calpha",
 167     "atoms only (file rms-ahx.xvg).[BR]",
 168     "[BB]8.[bb] Average Calpha-Calpha dihedral angle (file phi-ahx.xvg).[BR]",
 169     "[BB]9.[bb] Average Phi and Psi angles (file phipsi.xvg).[BR]",
 170     "[BB]10.[bb] Ellipticity at 222 nm according to [IT]Hirst and Brooks[it]",
 171     "[PAR]"
 172   };
 173   static bool bCheck=FALSE,bFit=TRUE,bDBG=FALSE,bEV=FALSE;
 174   static int  rStart=0,rEnd=0,r0=1;
 175   t_pargs pa [] = {
 176     { "-r0", FALSE, etINT, {&r0},
 177       "The first residue number in the sequence" },
 178     { "-q",  FALSE, etBOOL,{&bCheck},
 179       "Check at every step which part of the sequence is helical" },
 180     { "-F",  FALSE, etBOOL,{&bFit},
 181       "Toggle fit to a perfect helix" },
 182     { "-db", FALSE, etBOOL,{&bDBG},
 183       "Print debug info" },
 184     { "-ev", FALSE, etBOOL,{&bEV},
 185       "Write a new 'trajectory' file for ED" },
 186     { "-ahxstart", FALSE, etINT, {&rStart},
 187       "First residue in helix" },
 188     { "-ahxend", FALSE, etINT, {&rEnd},
 189       "Last residue in helix" }
 190   };
 191
 192   typedef struct {
 193     FILE *fp,*fp2;
 194     bool bfp2;
 195     char *filenm;
 196     char *title;
 197     char *xaxis;
 198     char *yaxis;
 199     real val;
 200   } t_xvgrfile;
 201
 202   t_xvgrfile xf[efhNR] = {
 203     { NULL, NULL, TRUE,  "radius",  "Helix radius",               NULL, "r (nm)" , 0.0 },
 204     { NULL, NULL, TRUE,  "twist",   "Twist per residue",          NULL, "Angle (deg)", 0.0 },
 205     { NULL, NULL, TRUE,  "rise",    "Rise per residue",           NULL, "Rise (nm)", 0.0 },
 206     { NULL, NULL, FALSE, "len-ahx", "Length of the Helix",        NULL, "Length (nm)", 0.0 },
 207     { NULL, NULL, FALSE, "dip-ahx", "Helix Backbone Dipole",      NULL, "rq (nm e)", 0.0 },
 208     { NULL, NULL, TRUE,  "rms-ahx", "RMS Deviation from Ideal Helix", NULL, "RMS (nm)", 0.0 },
 209     { NULL, NULL, FALSE, "rmsa-ahx","Average RMSD per Residue",   "Residue", "RMS (nm)", 0.0 },
 210     { NULL, NULL,FALSE,  "cd222",   "Ellipticity at 222 nm", NULL, "nm", 0.0 },
 211     { NULL, NULL, TRUE,  "pprms",   "RMS Distance from \\8a\\4-helix", NULL, "deg" , 0.0 },
 212     { NULL, NULL, TRUE,  "caphi",   "Average Ca-Ca Dihedral",     NULL, "\\8F\\4(deg)", 0.0 },
 213     { NULL, NULL, TRUE,  "phi",     "Average \\8F\\4 angles", NULL, "deg" , 0.0 },
 214     { NULL, NULL, TRUE,  "psi",     "Average \\8Y\\4 angles", NULL, "deg" , 0.0 },
 215     { NULL, NULL, TRUE,  "hb3",     "Average n-n+3 hbond length", NULL, "nm" , 0.0 },
 216     { NULL, NULL, TRUE,  "hb4",     "Average n-n+4 hbond length", NULL, "nm" , 0.0 },
 217     { NULL, NULL, TRUE,  "hb5",     "Average n-n+5 hbond length", NULL, "nm" , 0.0 },
 218     { NULL, NULL,FALSE,  "JCaHa",   "J-Coupling Values",        "Residue", "Hz" , 0.0 },
 219     { NULL, NULL,FALSE,  "helicity","Helicity per Residue",     "Residue", "% of time" , 0.0 }
 220   };
 221
 222   FILE       *otrj;
 223   char       buf[54],prop[256];
 224   int        status;
 225   int        natoms,nre,nres,step;
 226   t_bb       *bb;
 227   int        i,j,ai,m,nall,nbb,nca,teller,nSel=0;
 228   atom_id    *bbindex,*caindex,*allindex;
 229   t_topology *top;
 230   rvec       *x,*xref,*xav;
 231   real       t,lambda;
 232   real       rms,fac;
 233   matrix     box;
 234   bool       bRange;
 235   t_filenm  fnm[] = {
 236     { efTPX, NULL,  NULL,   ffREAD  },
 237     { efNDX, NULL,  NULL,   ffREAD  },
 238     { efTRX, "-f",  NULL,   ffREAD  },
 239     { efG87, "-to", NULL,   ffOPTWR },
 240     { efSTO, "-cz", "zconf",ffWRITE },
 241     { efSTO, "-co", "waver",ffWRITE }
 242   };
 243 #define NFILE asize(fnm)
 244
 245   CopyRight(stderr,argv[0]);
 246   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
 247                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 248   bRange=(opt2parg_bSet("-ahxstart",asize(pa),pa) &&
 249           opt2parg_bSet("-ahxend",asize(pa),pa));
 250
 251   top=read_top(ftp2fn(efTPX,NFILE,fnm));
 252
 253   natoms=read_first_x(&status,opt2fn("-f",NFILE,fnm),&t,&x,box);
 254
 255   if (opt2bSet("-to",NFILE,fnm)) {
 256     otrj=opt2FILE("-to",NFILE,fnm,"w");
 257     nSel=get_prop(prop);
 258     fprintf(otrj,"%s Weighted Trajectory: %d atoms, NO box\n",prop,natoms);
 259   }
 260   else
 261     otrj=NULL;
 262
 263   if (natoms != top->atoms.nr)
 264     fatal_error(0,"Sorry can only run when the number of atoms in the run input file (%d) is equal to the number in the trajectory (%d)",
 265             top->atoms.nr,natoms);
 266
 267   bb=mkbbind(ftp2fn(efNDX,NFILE,fnm),&nres,&nbb,r0,&nall,&allindex,
 268              top->atoms.atomname,top->atoms.atom,top->atoms.resname);
 269   snew(bbindex,natoms);
 270   snew(caindex,nres);
 271
 272   fprintf(stderr,"nall=%d\n",nall);
 273
 274   /* Open output files, default x-axis is time */
 275   for(i=0; (i<efhNR); i++) {
 276     sprintf(buf,"%s.xvg",xf[i].filenm);
 277     remove(buf);
 278     xf[i].fp=xvgropen(buf,xf[i].title,xf[i].xaxis ? xf[i].xaxis : "Time (ps)",
 279                       xf[i].yaxis);
 280     if (xf[i].bfp2) {
 281       sprintf(buf,"%s.out",xf[i].filenm);
 282       remove(buf);
 283       xf[i].fp2=ffopen(buf,"w");
 284     }
 285   }
 286
 287   /* Read reference frame from tpx file to compute helix length */
 288   snew(xref,top->atoms.nr);
 289   read_tpx(ftp2fn(efTPX,NFILE,fnm),
 290            &step,&t,&lambda,NULL,NULL,
 291            &natoms,xref,NULL,NULL,NULL);
 292   calc_hxprops(nres,bb,xref,box);
 293   do_start_end(nres,bb,xref,&nbb,bbindex,&nca,caindex,bRange,rStart,rEnd);
 294   sfree(xref);
 295   if (bDBG) {
 296     fprintf(stderr,"nca=%d, nbb=%d\n",nca,nbb);
 297     pr_bb(stdout,nres,bb);
 298   }
 299
 300   snew(xav,natoms);
 301   teller=0;
 302   do {
 303     if ((teller++ % 10) == 0)
 304       fprintf(stderr,"\rt=%.2f",t);
 305     rm_pbc(&(top->idef),top->atoms.nr,box,x,x);
 306
 307     calc_hxprops(nres,bb,x,box);
 308     if (bCheck)
 309       do_start_end(nres,bb,x,&nbb,bbindex,&nca,caindex,FALSE,0,0);
 310
 311     if (nca >= 5) {
 312       rms=fit_ahx(nres,bb,natoms,nall,allindex,x,nca,caindex,box,bFit);
 313
 314       if (teller == 1) {
 315         write_sto_conf(opt2fn("-cz",NFILE,fnm),"Helix fitted to Z-Axis",
 316                        &(top->atoms),x,NULL,box);
 317       }
 318
 319       xf[efhRAD].val   = radius(xf[efhRAD].fp2,nca,caindex,x);
 320       xf[efhTWIST].val = twist(xf[efhTWIST].fp2,nca,caindex,x);
 321       xf[efhRISE].val  = rise(nca,caindex,x);
 322       xf[efhLEN].val   = ahx_len(nca,caindex,x,box);
 323       xf[efhCD222].val = ellipticity(nres,bb);
 324       xf[efhDIP].val   = dip(nbb,bbindex,x,top->atoms.atom);
 325       xf[efhRMS].val   = rms;
 326       xf[efhCPHI].val  = ca_phi(nca,caindex,x,box);
 327       xf[efhPPRMS].val = pprms(xf[efhPPRMS].fp2,nres,bb);
 328
 329       for(j=0; (j<=efhCPHI); j++)
 330         fprintf(xf[j].fp,   "%10g  %10g\n",t,xf[j].val);
 331
 332       av_phipsi(xf[efhPHI].fp,xf[efhPSI].fp,xf[efhPHI].fp2,xf[efhPSI].fp2,
 333                 t,nres,bb);
 334       av_hblen(xf[efhHB3].fp,xf[efhHB3].fp2,
 335                xf[efhHB4].fp,xf[efhHB4].fp2,
 336                xf[efhHB5].fp,xf[efhHB5].fp2,
 337                t,nres,bb);
 338
 339       if (otrj)
 340         dump_otrj(otrj,nall,allindex,x,xf[nSel].val,xav);
 341     }
 342   } while (read_next_x(status,&t,natoms,x,box));
 343   fprintf(stderr,"\n");
 344
 345   close_trj(status);
 346
 347   if (otrj) {
 348     fclose(otrj);
 349     fac=1.0/teller;
 350     for(i=0; (i<nall); i++) {
 351       ai=allindex[i];
 352       for(m=0; (m<DIM); m++)
 353         xav[ai][m]*=fac;
 354     }
 355     write_sto_conf_indexed(opt2fn("-co",NFILE,fnm),
 356                            "Weighted and Averaged conformation",
 357                            &(top->atoms),xav,NULL,box,nall,allindex);
 358   }
 359
 360   for(i=0; (i<nres); i++) {
 361     if (bb[i].nrms > 0) {
 362       fprintf(xf[efhRMSA].fp,"%10d  %10g\n",r0+i,bb[i].rmsa/bb[i].nrms);
 363     }
 364     fprintf(xf[efhAHX].fp,"%10d  %10g\n",r0+i,(bb[i].nhx*100.0)/(real )teller);
 365     fprintf(xf[efhJCA].fp,"%10d  %10g\n",
 366             r0+i,140.3+(bb[i].jcaha/(double)teller));
 367   }
 368
 369   for(i=0; (i<efhNR); i++) {
 370     fclose(xf[i].fp);
 371     if (xf[i].bfp2)
 372       fclose(xf[i].fp2);
 373     xvgr_file(xf[i].filenm,NULL);
 374   }
 375
 376   thanx(stdout);
 377
 378   return 0;
 379 }
 380