src/gromacs/gmxana/gmx_helix.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40
  41 #include "gromacs/fileio/confio.h"
  42 #include "gmx_fatal.h"
  43 #include "fitahx.h"
  44 #include "gromacs/fileio/futil.h"
  45 #include "gstat.h"
  46 #include "gromacs/fileio/g87io.h"
  47 #include "hxprops.h"
  48 #include "macros.h"
  49 #include "maths.h"
  50 #include "pbc.h"
  51 #include "gromacs/fileio/tpxio.h"
  52 #include "gromacs/fileio/trxio.h"
  53 #include "physics.h"
  54 #include "index.h"
  55 #include "smalloc.h"
  56 #include "statutil.h"
  57 #include <string.h>
  58 #include "sysstuff.h"
  59 #include "txtdump.h"
  60 #include "typedefs.h"
  61 #include "vec.h"
  62 #include "xvgr.h"
  63 #include "gmx_ana.h"
  64
  65 void dump_otrj(FILE *otrj, int natoms, atom_id all_index[], rvec x[],
  66                real fac, rvec xav[])
  67 {
  68     static FILE *fp   = NULL;
  69     static real  fac0 = 1.0;
  70
  71     int          i, ai, m;
  72     real         xa, xm;
  73
  74     if (fp == NULL)
  75     {
  76         fp   = ffopen("WEDGAMMA10.DAT", "w");
  77         fac0 = fac;
  78     }
  79     fac /= fac0;
  80
  81     fprintf(fp, "%10g\n", fac);
  82
  83     for (i = 0; (i < natoms); i++)
  84     {
  85         ai = all_index[i];
  86         for (m = 0; (m < DIM); m++)
  87         {
  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_g87_ndx(otrj, natoms, all_index, x, NULL);
  95 }
  96
  97 int gmx_helix(int argc, char *argv[])
  98 {
  99     const char        *desc[] = {
 100         "[THISMODULE] 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] Helix dipole, backbone only (file [TT]dip-ahx.xvg[tt]).[BR]",
 122         "[BB]6.[bb] RMS deviation from ideal helix, calculated for the C[GRK]alpha[grk]",
 123         "atoms only (file [TT]rms-ahx.xvg[tt]).[BR]",
 124         "[BB]7.[bb] Average C[GRK]alpha[grk] - C[GRK]alpha[grk] dihedral angle (file [TT]phi-ahx.xvg[tt]).[BR]",
 125         "[BB]8.[bb] Average [GRK]phi[grk] and [GRK]psi[grk] angles (file [TT]phipsi.xvg[tt]).[BR]",
 126         "[BB]9.[bb] Ellipticity at 222 nm according to Hirst and Brooks.",
 127         "[PAR]"
 128     };
 129     static const char *ppp[efhNR+2] = {
 130         NULL, "RAD", "TWIST", "RISE", "LEN", "NHX", "DIP", "RMS", "CPHI",
 131         "RMSA", "PHI", "PSI", "HB3", "HB4", "HB5", "CD222", NULL
 132     };
 133     static gmx_bool    bCheck = FALSE, bFit = TRUE, bDBG = FALSE, bEV = FALSE;
 134     static int         rStart = 0, rEnd = 0, r0 = 1;
 135     t_pargs            pa []  = {
 136         { "-r0", FALSE, etINT, {&r0},
 137           "The first residue number in the sequence" },
 138         { "-q",  FALSE, etBOOL, {&bCheck},
 139           "Check at every step which part of the sequence is helical" },
 140         { "-F",  FALSE, etBOOL, {&bFit},
 141           "Toggle fit to a perfect helix" },
 142         { "-db", FALSE, etBOOL, {&bDBG},
 143           "Print debug info" },
 144         { "-prop", FALSE, etENUM, {ppp},
 145           "Select property to weight eigenvectors with. WARNING experimental stuff" },
 146         { "-ev", FALSE, etBOOL, {&bEV},
 147           "Write a new 'trajectory' file for ED" },
 148         { "-ahxstart", FALSE, etINT, {&rStart},
 149           "First residue in helix" },
 150         { "-ahxend", FALSE, etINT, {&rEnd},
 151           "Last residue in helix" }
 152     };
 153
 154     typedef struct {
 155         FILE       *fp, *fp2;
 156         gmx_bool    bfp2;
 157         const char *filenm;
 158         const char *title;
 159         const char *xaxis;
 160         const char *yaxis;
 161         real        val;
 162     } t_xvgrfile;
 163
 164     t_xvgrfile     xf[efhNR] = {
 165         { NULL, NULL, TRUE,  "radius",  "Helix radius",               NULL, "r (nm)", 0.0 },
 166         { NULL, NULL, TRUE,  "twist",   "Twist per residue",          NULL, "Angle (deg)", 0.0 },
 167         { NULL, NULL, TRUE,  "rise",    "Rise per residue",           NULL, "Rise (nm)", 0.0 },
 168         { NULL, NULL, FALSE, "len-ahx", "Length of the Helix",        NULL, "Length (nm)", 0.0 },
 169         { NULL, NULL, FALSE, "dip-ahx", "Helix Backbone Dipole",      NULL, "rq (nm e)", 0.0 },
 170         { NULL, NULL, TRUE,  "rms-ahx", "RMS Deviation from Ideal Helix", NULL, "RMS (nm)", 0.0 },
 171         { NULL, NULL, FALSE, "rmsa-ahx", "Average RMSD per Residue",   "Residue", "RMS (nm)", 0.0 },
 172         { NULL, NULL, FALSE,  "cd222",   "Ellipticity at 222 nm", NULL, "nm", 0.0 },
 173         { NULL, NULL, TRUE,  "pprms",   "RMS Distance from \\8a\\4-helix", NULL, "deg", 0.0 },
 174         { NULL, NULL, TRUE,  "caphi",   "Average Ca-Ca Dihedral",     NULL, "\\8F\\4(deg)", 0.0 },
 175         { NULL, NULL, TRUE,  "phi",     "Average \\8F\\4 angles", NULL, "deg", 0.0 },
 176         { NULL, NULL, TRUE,  "psi",     "Average \\8Y\\4 angles", NULL, "deg", 0.0 },
 177         { NULL, NULL, TRUE,  "hb3",     "Average n-n+3 hbond length", NULL, "nm", 0.0 },
 178         { NULL, NULL, TRUE,  "hb4",     "Average n-n+4 hbond length", NULL, "nm", 0.0 },
 179         { NULL, NULL, TRUE,  "hb5",     "Average n-n+5 hbond length", NULL, "nm", 0.0 },
 180         { NULL, NULL, FALSE,  "JCaHa",   "J-Coupling Values",        "Residue", "Hz", 0.0 },
 181         { NULL, NULL, FALSE,  "helicity", "Helicity per Residue",     "Residue", "% of time", 0.0 }
 182     };
 183
 184     output_env_t   oenv;
 185     FILE          *otrj;
 186     char           buf[54], prop[256];
 187     t_trxstatus   *status;
 188     int            natoms, nre, nres;
 189     t_bb          *bb;
 190     int            i, j, ai, m, nall, nbb, nca, teller, nSel = 0;
 191     atom_id       *bbindex, *caindex, *allindex;
 192     t_topology    *top;
 193     int            ePBC;
 194     rvec          *x, *xref, *xav;
 195     real           t;
 196     real           rms, fac;
 197     matrix         box;
 198     gmx_rmpbc_t    gpbc = NULL;
 199     gmx_bool       bRange;
 200     t_filenm       fnm[] = {
 201         { efTPX, NULL,  NULL,   ffREAD  },
 202         { efNDX, NULL,  NULL,   ffREAD  },
 203         { efTRX, "-f",  NULL,   ffREAD  },
 204         { efG87, "-to", NULL,   ffOPTWR },
 205         { efSTO, "-cz", "zconf", ffWRITE },
 206         { efSTO, "-co", "waver", ffWRITE }
 207     };
 208 #define NFILE asize(fnm)
 209
 210     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 211                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 212     {
 213         return 0;
 214     }
 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     {
 225         otrj = opt2FILE("-to", NFILE, fnm, "w");
 226         strcpy(prop, ppp[0]);
 227         fprintf(otrj, "%s Weighted Trajectory: %d atoms, NO box\n", prop, natoms);
 228     }
 229     else
 230     {
 231         otrj = NULL;
 232     }
 233
 234     if (natoms != top->atoms.nr)
 235     {
 236         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)",
 237                   top->atoms.nr, natoms);
 238     }
 239
 240     bb = mkbbind(ftp2fn(efNDX, NFILE, fnm), &nres, &nbb, r0, &nall, &allindex,
 241                  top->atoms.atomname, top->atoms.atom, top->atoms.resinfo);
 242     snew(bbindex, natoms);
 243     snew(caindex, nres);
 244
 245     fprintf(stderr, "nall=%d\n", nall);
 246
 247     /* Open output files, default x-axis is time */
 248     for (i = 0; (i < efhNR); i++)
 249     {
 250         sprintf(buf, "%s.xvg", xf[i].filenm);
 251         remove(buf);
 252         xf[i].fp = xvgropen(buf, xf[i].title,
 253                             xf[i].xaxis ? xf[i].xaxis : "Time (ps)",
 254                             xf[i].yaxis, oenv);
 255         if (xf[i].bfp2)
 256         {
 257             sprintf(buf, "%s.out", xf[i].filenm);
 258             remove(buf);
 259             xf[i].fp2 = ffopen(buf, "w");
 260         }
 261     }
 262
 263     /* Read reference frame from tpx file to compute helix length */
 264     snew(xref, top->atoms.nr);
 265     read_tpx(ftp2fn(efTPX, NFILE, fnm),
 266              NULL, NULL, &natoms, xref, NULL, NULL, NULL);
 267     calc_hxprops(nres, bb, xref);
 268     do_start_end(nres, bb, &nbb, bbindex, &nca, caindex, bRange, rStart, rEnd);
 269     sfree(xref);
 270     if (bDBG)
 271     {
 272         fprintf(stderr, "nca=%d, nbb=%d\n", nca, nbb);
 273         pr_bb(stdout, nres, bb);
 274     }
 275
 276     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
 277
 278     snew(xav, natoms);
 279     teller = 0;
 280     do
 281     {
 282         if ((teller++ % 10) == 0)
 283         {
 284             fprintf(stderr, "\rt=%.2f", t);
 285         }
 286         gmx_rmpbc(gpbc, natoms, box, x);
 287
 288
 289         calc_hxprops(nres, bb, x);
 290         if (bCheck)
 291         {
 292             do_start_end(nres, bb, &nbb, bbindex, &nca, caindex, FALSE, 0, 0);
 293         }
 294
 295         if (nca >= 5)
 296         {
 297             rms = fit_ahx(nres, bb, natoms, nall, allindex, x, nca, caindex, bFit);
 298
 299             if (teller == 1)
 300             {
 301                 write_sto_conf(opt2fn("-cz", NFILE, fnm), "Helix fitted to Z-Axis",
 302                                &(top->atoms), x, NULL, ePBC, box);
 303             }
 304
 305             xf[efhRAD].val   = radius(xf[efhRAD].fp2, nca, caindex, x);
 306             xf[efhTWIST].val = twist(nca, caindex, x);
 307             xf[efhRISE].val  = rise(nca, caindex, x);
 308             xf[efhLEN].val   = ahx_len(nca, caindex, x);
 309             xf[efhCD222].val = ellipticity(nres, bb);
 310             xf[efhDIP].val   = dip(nbb, bbindex, x, top->atoms.atom);
 311             xf[efhRMS].val   = rms;
 312             xf[efhCPHI].val  = ca_phi(nca, caindex, x);
 313             xf[efhPPRMS].val = pprms(xf[efhPPRMS].fp2, nres, bb);
 314
 315             for (j = 0; (j <= efhCPHI); j++)
 316             {
 317                 fprintf(xf[j].fp,   "%10g  %10g\n", t, xf[j].val);
 318             }
 319
 320             av_phipsi(xf[efhPHI].fp, xf[efhPSI].fp, xf[efhPHI].fp2, xf[efhPSI].fp2,
 321                       t, nres, bb);
 322             av_hblen(xf[efhHB3].fp, xf[efhHB3].fp2,
 323                      xf[efhHB4].fp, xf[efhHB4].fp2,
 324                      xf[efhHB5].fp, xf[efhHB5].fp2,
 325                      t, nres, bb);
 326
 327             if (otrj)
 328             {
 329                 dump_otrj(otrj, nall, allindex, x, xf[nSel].val, xav);
 330             }
 331         }
 332     }
 333     while (read_next_x(oenv, status, &t, x, box));
 334     fprintf(stderr, "\n");
 335
 336     gmx_rmpbc_done(gpbc);
 337
 338     close_trj(status);
 339
 340     if (otrj)
 341     {
 342         ffclose(otrj);
 343         fac = 1.0/teller;
 344         for (i = 0; (i < nall); i++)
 345         {
 346             ai = allindex[i];
 347             for (m = 0; (m < DIM); m++)
 348             {
 349                 xav[ai][m] *= fac;
 350             }
 351         }
 352         write_sto_conf_indexed(opt2fn("-co", NFILE, fnm),
 353                                "Weighted and Averaged conformation",
 354                                &(top->atoms), xav, NULL, ePBC, box, nall, allindex);
 355     }
 356
 357     for (i = 0; (i < nres); i++)
 358     {
 359         if (bb[i].nrms > 0)
 360         {
 361             fprintf(xf[efhRMSA].fp, "%10d  %10g\n", r0+i, bb[i].rmsa/bb[i].nrms);
 362         }
 363         fprintf(xf[efhAHX].fp, "%10d  %10g\n", r0+i, (bb[i].nhx*100.0)/(real )teller);
 364         fprintf(xf[efhJCA].fp, "%10d  %10g\n",
 365                 r0+i, 140.3+(bb[i].jcaha/(double)teller));
 366     }
 367
 368     for (i = 0; (i < efhNR); i++)
 369     {
 370         ffclose(xf[i].fp);
 371         if (xf[i].bfp2)
 372         {
 373             ffclose(xf[i].fp2);
 374         }
 375         do_view(oenv, xf[i].filenm, "-nxy");
 376     }
 377
 378     return 0;
 379 }