src/gromacs/gmxana/gmx_helix.cpp

   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  * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cmath>
  40 #include <cstring>
  41
  42 #include "gromacs/commandline/pargs.h"
  43 #include "gromacs/commandline/viewit.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/fileio/tpxio.h"
  46 #include "gromacs/fileio/trxio.h"
  47 #include "gromacs/fileio/xvgr.h"
  48 #include "gromacs/gmxana/fitahx.h"
  49 #include "gromacs/gmxana/gmx_ana.h"
  50 #include "gromacs/gmxana/gstat.h"
  51 #include "gromacs/gmxana/hxprops.h"
  52 #include "gromacs/math/utilities.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/pbcutil/rmpbc.h"
  55 #include "gromacs/topology/topology.h"
  56 #include "gromacs/utility/arraysize.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/futil.h"
  59 #include "gromacs/utility/smalloc.h"
  60
  61 int gmx_helix(int argc, char* argv[])
  62 {
  63     const char* desc[] = {
  64         "[THISMODULE] computes all kinds of helix properties. First, the peptide",
  65         "is checked to find the longest helical part, as determined by",
  66         "hydrogen bonds and [GRK]phi[grk]/[GRK]psi[grk] angles.",
  67         "That bit is fitted",
  68         "to an ideal helix around the [IT]z[it]-axis and centered around the origin.",
  69         "Then the following properties are computed:",
  70         "",
  71         " * Helix radius (file [TT]radius.xvg[tt]). This is merely the",
  72         "   RMS deviation in two dimensions for all C[GRK]alpha[grk] atoms.",
  73         "   it is calculated as [SQRT]([SUM][sum][SUB]i[sub] (x^2(i)+y^2(i)))/N[sqrt] where N is ",
  74         "   the number of backbone atoms. For an ideal helix the radius is 0.23 nm.",
  75         " * Twist (file [TT]twist.xvg[tt]). The average helical angle per",
  76         "   residue is calculated. For an [GRK]alpha[grk]-helix it is 100 degrees,",
  77         "   for 3-10 helices it will be smaller, and ",
  78         "   for 5-helices it will be larger.",
  79         " * Rise per residue (file [TT]rise.xvg[tt]). The helical rise per",
  80         "   residue is plotted as the difference in [IT]z[it]-coordinate between C[GRK]alpha[grk]",
  81         "   atoms. For an ideal helix, this is 0.15 nm.",
  82         " * Total helix length (file [TT]len-ahx.xvg[tt]). The total length",
  83         "   of the",
  84         "   helix in nm. This is simply the average rise (see above) times the",
  85         "   number of helical residues (see below).",
  86         " * Helix dipole, backbone only (file [TT]dip-ahx.xvg[tt]).",
  87         " * RMS deviation from ideal helix, calculated for the C[GRK]alpha[grk]",
  88         "   atoms only (file [TT]rms-ahx.xvg[tt]).",
  89         " * Average C[GRK]alpha[grk] - C[GRK]alpha[grk] dihedral angle (file [TT]phi-ahx.xvg[tt]).",
  90         " * Average [GRK]phi[grk] and [GRK]psi[grk] angles (file [TT]phipsi.xvg[tt]).",
  91         " * Ellipticity at 222 nm according to Hirst and Brooks."
  92     };
  93     static gmx_bool bCheck = FALSE, bFit = TRUE, bDBG = FALSE, bEV = FALSE;
  94     static int      rStart = 0, rEnd = 0, r0 = 1;
  95     t_pargs pa[] = { { "-r0", FALSE, etINT, { &r0 }, "The first residue number in the sequence" },
  96                      { "-q",
  97                        FALSE,
  98                        etBOOL,
  99                        { &bCheck },
 100                        "Check at every step which part of the sequence is helical" },
 101                      { "-F", FALSE, etBOOL, { &bFit }, "Toggle fit to a perfect helix" },
 102                      { "-db", FALSE, etBOOL, { &bDBG }, "Print debug info" },
 103                      { "-ev", FALSE, etBOOL, { &bEV }, "Write a new 'trajectory' file for ED" },
 104                      { "-ahxstart", FALSE, etINT, { &rStart }, "First residue in helix" },
 105                      { "-ahxend", FALSE, etINT, { &rEnd }, "Last residue in helix" } };
 106
 107     typedef struct
 108     { //NOLINT(clang-analyzer-optin.performance.Padding)
 109         FILE *      fp, *fp2;
 110         gmx_bool    bfp2;
 111         const char* filenm;
 112         const char* title;
 113         const char* xaxis;
 114         const char* yaxis;
 115         real        val;
 116     } t_xvgrfile;
 117
 118     t_xvgrfile xf[efhNR] = {
 119         { nullptr, nullptr, TRUE, "radius", "Helix radius", nullptr, "r (nm)", 0.0 },
 120         { nullptr, nullptr, TRUE, "twist", "Twist per residue", nullptr, "Angle (deg)", 0.0 },
 121         { nullptr, nullptr, TRUE, "rise", "Rise per residue", nullptr, "Rise (nm)", 0.0 },
 122         { nullptr, nullptr, FALSE, "len-ahx", "Length of the Helix", nullptr, "Length (nm)", 0.0 },
 123         { nullptr, nullptr, FALSE, "dip-ahx", "Helix Backbone Dipole", nullptr, "rq (nm e)", 0.0 },
 124         { nullptr, nullptr, TRUE, "rms-ahx", "RMS Deviation from Ideal Helix", nullptr, "RMS (nm)", 0.0 },
 125         { nullptr, nullptr, FALSE, "rmsa-ahx", "Average RMSD per Residue", "Residue", "RMS (nm)", 0.0 },
 126         { nullptr, nullptr, FALSE, "cd222", "Ellipticity at 222 nm", nullptr, "nm", 0.0 },
 127         { nullptr, nullptr, TRUE, "pprms", "RMS Distance from \\8a\\4-helix", nullptr, "deg", 0.0 },
 128         { nullptr, nullptr, TRUE, "caphi", "Average Ca-Ca Dihedral", nullptr, "\\8F\\4(deg)", 0.0 },
 129         { nullptr, nullptr, TRUE, "phi", "Average \\8F\\4 angles", nullptr, "deg", 0.0 },
 130         { nullptr, nullptr, TRUE, "psi", "Average \\8Y\\4 angles", nullptr, "deg", 0.0 },
 131         { nullptr, nullptr, TRUE, "hb3", "Average n-n+3 hbond length", nullptr, "nm", 0.0 },
 132         { nullptr, nullptr, TRUE, "hb4", "Average n-n+4 hbond length", nullptr, "nm", 0.0 },
 133         { nullptr, nullptr, TRUE, "hb5", "Average n-n+5 hbond length", nullptr, "nm", 0.0 },
 134         { nullptr, nullptr, FALSE, "JCaHa", "J-Coupling Values", "Residue", "Hz", 0.0 },
 135         { nullptr, nullptr, FALSE, "helicity", "Helicity per Residue", "Residue", "% of time", 0.0 }
 136     };
 137
 138     gmx_output_env_t* oenv;
 139     char              buf[54];
 140     t_trxstatus*      status;
 141     int               natoms, nres;
 142     t_bb*             bb;
 143     int               i, j, nall, nbb, nca, teller;
 144     int *             bbindex, *caindex, *allindex;
 145     t_topology*       top;
 146     int               ePBC;
 147     rvec *            x, *xref;
 148     real              t;
 149     real              rms;
 150     matrix            box;
 151     gmx_rmpbc_t       gpbc = nullptr;
 152     gmx_bool          bRange;
 153     t_filenm          fnm[] = {
 154         { efTPR, nullptr, nullptr, ffREAD },
 155         { efNDX, nullptr, nullptr, ffREAD },
 156         { efTRX, "-f", nullptr, ffREAD },
 157         { efSTO, "-cz", "zconf", ffWRITE },
 158     };
 159 #define NFILE asize(fnm)
 160
 161     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa,
 162                            asize(desc), desc, 0, nullptr, &oenv))
 163     {
 164         return 0;
 165     }
 166
 167     bRange = (opt2parg_bSet("-ahxstart", asize(pa), pa) && opt2parg_bSet("-ahxend", asize(pa), pa));
 168
 169     top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC);
 170
 171     natoms = read_first_x(oenv, &status, opt2fn("-f", NFILE, fnm), &t, &x, box);
 172
 173     if (natoms != top->atoms.nr)
 174     {
 175         gmx_fatal(FARGS,
 176                   "Sorry can only run when the number of atoms in the run input file (%d) is equal "
 177                   "to the number in the trajectory (%d)",
 178                   top->atoms.nr, natoms);
 179     }
 180
 181     bb = mkbbind(ftp2fn(efNDX, NFILE, fnm), &nres, &nbb, r0, &nall, &allindex, top->atoms.atomname,
 182                  top->atoms.atom, top->atoms.resinfo);
 183     snew(bbindex, natoms);
 184     snew(caindex, nres);
 185
 186     fprintf(stderr, "nall=%d\n", nall);
 187
 188     /* Open output files, default x-axis is time */
 189     for (i = 0; (i < efhNR); i++)
 190     {
 191         sprintf(buf, "%s.xvg", xf[i].filenm);
 192         remove(buf);
 193         xf[i].fp = xvgropen(buf, xf[i].title, xf[i].xaxis ? xf[i].xaxis : "Time (ps)", xf[i].yaxis, oenv);
 194         if (xf[i].bfp2)
 195         {
 196             sprintf(buf, "%s.out", xf[i].filenm);
 197             remove(buf);
 198             xf[i].fp2 = gmx_ffopen(buf, "w");
 199         }
 200     }
 201
 202     /* Read reference frame from tpx file to compute helix length */
 203     snew(xref, top->atoms.nr);
 204     read_tpx(ftp2fn(efTPR, NFILE, fnm), nullptr, nullptr, nullptr, xref, nullptr, nullptr);
 205     calc_hxprops(nres, bb, xref);
 206     do_start_end(nres, bb, &nbb, bbindex, &nca, caindex, bRange, rStart, rEnd);
 207     sfree(xref);
 208     if (bDBG)
 209     {
 210         fprintf(stderr, "nca=%d, nbb=%d\n", nca, nbb);
 211         pr_bb(stdout, nres, bb);
 212     }
 213
 214     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
 215
 216     teller = 0;
 217     do
 218     {
 219         if ((teller++ % 10) == 0)
 220         {
 221             fprintf(stderr, "\rt=%.2f", t);
 222             fflush(stderr);
 223         }
 224         gmx_rmpbc(gpbc, natoms, box, x);
 225
 226
 227         calc_hxprops(nres, bb, x);
 228         if (bCheck)
 229         {
 230             do_start_end(nres, bb, &nbb, bbindex, &nca, caindex, FALSE, 0, 0);
 231         }
 232
 233         if (nca >= 5)
 234         {
 235             rms = fit_ahx(nres, bb, natoms, nall, allindex, x, nca, caindex, bFit);
 236
 237             if (teller == 1)
 238             {
 239                 write_sto_conf(opt2fn("-cz", NFILE, fnm), "Helix fitted to Z-Axis", &(top->atoms),
 240                                x, nullptr, ePBC, box);
 241             }
 242
 243             xf[efhRAD].val   = radius(xf[efhRAD].fp2, nca, caindex, x);
 244             xf[efhTWIST].val = twist(nca, caindex, x);
 245             xf[efhRISE].val  = rise(nca, caindex, x);
 246             xf[efhLEN].val   = ahx_len(nca, caindex, x);
 247             xf[efhCD222].val = ellipticity(nres, bb);
 248             xf[efhDIP].val   = dip(nbb, bbindex, x, top->atoms.atom);
 249             xf[efhRMS].val   = rms;
 250             xf[efhCPHI].val  = ca_phi(nca, caindex, x);
 251             xf[efhPPRMS].val = pprms(xf[efhPPRMS].fp2, nres, bb);
 252
 253             for (j = 0; (j <= efhCPHI); j++)
 254             {
 255                 fprintf(xf[j].fp, "%10g  %10g\n", t, xf[j].val);
 256             }
 257
 258             av_phipsi(xf[efhPHI].fp, xf[efhPSI].fp, xf[efhPHI].fp2, xf[efhPSI].fp2, t, nres, bb);
 259             av_hblen(xf[efhHB3].fp, xf[efhHB3].fp2, xf[efhHB4].fp, xf[efhHB4].fp2, xf[efhHB5].fp,
 260                      xf[efhHB5].fp2, t, nres, bb);
 261         }
 262     } while (read_next_x(oenv, status, &t, x, box));
 263     fprintf(stderr, "\n");
 264
 265     gmx_rmpbc_done(gpbc);
 266
 267     close_trx(status);
 268
 269     for (i = 0; (i < nres); i++)
 270     {
 271         if (bb[i].nrms > 0)
 272         {
 273             fprintf(xf[efhRMSA].fp, "%10d  %10g\n", r0 + i, bb[i].rmsa / bb[i].nrms);
 274         }
 275         fprintf(xf[efhAHX].fp, "%10d  %10g\n", r0 + i, (bb[i].nhx * 100.0) / static_cast<real>(teller));
 276         fprintf(xf[efhJCA].fp, "%10d  %10g\n", r0 + i,
 277                 140.3 + (bb[i].jcaha / static_cast<double>(teller)));
 278     }
 279
 280     for (i = 0; (i < efhNR); i++)
 281     {
 282         xvgrclose(xf[i].fp);
 283         if (xf[i].bfp2)
 284         {
 285             xvgrclose(xf[i].fp2);
 286         }
 287         do_view(oenv, xf[i].filenm, "-nxy");
 288     }
 289
 290     return 0;
 291 }