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