2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, 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.
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.
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.
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.
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.
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.
42 #include "gromacs/fileio/confio.h"
44 #include "gromacs/utility/futil.h"
47 #include "gromacs/legacyheaders/macros.h"
48 #include "gromacs/math/utilities.h"
49 #include "gromacs/fileio/tpxio.h"
50 #include "gromacs/fileio/trxio.h"
51 #include "gromacs/utility/smalloc.h"
52 #include "gromacs/legacyheaders/txtdump.h"
53 #include "gromacs/legacyheaders/typedefs.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/legacyheaders/viewit.h"
58 #include "gromacs/commandline/pargs.h"
59 #include "gromacs/fileio/xvgr.h"
60 #include "gromacs/pbcutil/rmpbc.h"
61 #include "gromacs/utility/fatalerror.h"
63 int gmx_helix(int argc, char *argv[])
65 const char *desc[] = {
66 "[THISMODULE] computes all kinds of helix properties. First, the peptide",
67 "is checked to find the longest helical part, as determined by",
68 "hydrogen bonds and [GRK]phi[grk]/[GRK]psi[grk] angles.",
70 "to an ideal helix around the [IT]z[it]-axis and centered around the origin.",
71 "Then the following properties are computed:[PAR]",
72 "[BB]1.[bb] Helix radius (file [TT]radius.xvg[tt]). This is merely the",
73 "RMS deviation in two dimensions for all C[GRK]alpha[grk] atoms.",
74 "it is calculated as [SQRT]([SUM][sum][SUB]i[sub] (x^2(i)+y^2(i)))/N[sqrt] where N is the number",
75 "of backbone atoms. For an ideal helix the radius is 0.23 nm[BR]",
76 "[BB]2.[bb] Twist (file [TT]twist.xvg[tt]). The average helical angle per",
77 "residue is calculated. For an [GRK]alpha[grk]-helix it is 100 degrees,",
78 "for 3-10 helices it will be smaller, and ",
79 "for 5-helices it will be larger.[BR]",
80 "[BB]3.[bb] Rise per residue (file [TT]rise.xvg[tt]). The helical rise per",
81 "residue is plotted as the difference in [IT]z[it]-coordinate between C[GRK]alpha[grk]",
82 "atoms. For an ideal helix, this is 0.15 nm[BR]",
83 "[BB]4.[bb] Total helix length (file [TT]len-ahx.xvg[tt]). The total length",
85 "helix in nm. This is simply the average rise (see above) times the",
86 "number of helical residues (see below).[BR]",
87 "[BB]5.[bb] Helix dipole, backbone only (file [TT]dip-ahx.xvg[tt]).[BR]",
88 "[BB]6.[bb] RMS deviation from ideal helix, calculated for the C[GRK]alpha[grk]",
89 "atoms only (file [TT]rms-ahx.xvg[tt]).[BR]",
90 "[BB]7.[bb] Average C[GRK]alpha[grk] - C[GRK]alpha[grk] dihedral angle (file [TT]phi-ahx.xvg[tt]).[BR]",
91 "[BB]8.[bb] Average [GRK]phi[grk] and [GRK]psi[grk] angles (file [TT]phipsi.xvg[tt]).[BR]",
92 "[BB]9.[bb] Ellipticity at 222 nm according to Hirst and Brooks.",
95 static gmx_bool bCheck = FALSE, bFit = TRUE, bDBG = FALSE, bEV = FALSE;
96 static int rStart = 0, rEnd = 0, r0 = 1;
98 { "-r0", FALSE, etINT, {&r0},
99 "The first residue number in the sequence" },
100 { "-q", FALSE, etBOOL, {&bCheck},
101 "Check at every step which part of the sequence is helical" },
102 { "-F", FALSE, etBOOL, {&bFit},
103 "Toggle fit to a perfect helix" },
104 { "-db", FALSE, etBOOL, {&bDBG},
105 "Print debug info" },
106 { "-ev", FALSE, etBOOL, {&bEV},
107 "Write a new 'trajectory' file for ED" },
108 { "-ahxstart", FALSE, etINT, {&rStart},
109 "First residue in helix" },
110 { "-ahxend", FALSE, etINT, {&rEnd},
111 "Last residue in helix" }
124 t_xvgrfile xf[efhNR] = {
125 { NULL, NULL, TRUE, "radius", "Helix radius", NULL, "r (nm)", 0.0 },
126 { NULL, NULL, TRUE, "twist", "Twist per residue", NULL, "Angle (deg)", 0.0 },
127 { NULL, NULL, TRUE, "rise", "Rise per residue", NULL, "Rise (nm)", 0.0 },
128 { NULL, NULL, FALSE, "len-ahx", "Length of the Helix", NULL, "Length (nm)", 0.0 },
129 { NULL, NULL, FALSE, "dip-ahx", "Helix Backbone Dipole", NULL, "rq (nm e)", 0.0 },
130 { NULL, NULL, TRUE, "rms-ahx", "RMS Deviation from Ideal Helix", NULL, "RMS (nm)", 0.0 },
131 { NULL, NULL, FALSE, "rmsa-ahx", "Average RMSD per Residue", "Residue", "RMS (nm)", 0.0 },
132 { NULL, NULL, FALSE, "cd222", "Ellipticity at 222 nm", NULL, "nm", 0.0 },
133 { NULL, NULL, TRUE, "pprms", "RMS Distance from \\8a\\4-helix", NULL, "deg", 0.0 },
134 { NULL, NULL, TRUE, "caphi", "Average Ca-Ca Dihedral", NULL, "\\8F\\4(deg)", 0.0 },
135 { NULL, NULL, TRUE, "phi", "Average \\8F\\4 angles", NULL, "deg", 0.0 },
136 { NULL, NULL, TRUE, "psi", "Average \\8Y\\4 angles", NULL, "deg", 0.0 },
137 { NULL, NULL, TRUE, "hb3", "Average n-n+3 hbond length", NULL, "nm", 0.0 },
138 { NULL, NULL, TRUE, "hb4", "Average n-n+4 hbond length", NULL, "nm", 0.0 },
139 { NULL, NULL, TRUE, "hb5", "Average n-n+5 hbond length", NULL, "nm", 0.0 },
140 { NULL, NULL, FALSE, "JCaHa", "J-Coupling Values", "Residue", "Hz", 0.0 },
141 { NULL, NULL, FALSE, "helicity", "Helicity per Residue", "Residue", "% of time", 0.0 }
147 int natoms, nre, nres;
149 int i, j, m, nall, nbb, nca, teller, nSel = 0;
150 atom_id *bbindex, *caindex, *allindex;
157 gmx_rmpbc_t gpbc = NULL;
160 { efTPX, NULL, NULL, ffREAD },
161 { efNDX, NULL, NULL, ffREAD },
162 { efTRX, "-f", NULL, ffREAD },
163 { efSTO, "-cz", "zconf", ffWRITE },
165 #define NFILE asize(fnm)
167 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
168 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
173 bRange = (opt2parg_bSet("-ahxstart", asize(pa), pa) &&
174 opt2parg_bSet("-ahxend", asize(pa), pa));
176 top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC);
178 natoms = read_first_x(oenv, &status, opt2fn("-f", NFILE, fnm), &t, &x, box);
180 if (natoms != top->atoms.nr)
182 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)",
183 top->atoms.nr, natoms);
186 bb = mkbbind(ftp2fn(efNDX, NFILE, fnm), &nres, &nbb, r0, &nall, &allindex,
187 top->atoms.atomname, top->atoms.atom, top->atoms.resinfo);
188 snew(bbindex, natoms);
191 fprintf(stderr, "nall=%d\n", nall);
193 /* Open output files, default x-axis is time */
194 for (i = 0; (i < efhNR); i++)
196 sprintf(buf, "%s.xvg", xf[i].filenm);
198 xf[i].fp = xvgropen(buf, xf[i].title,
199 xf[i].xaxis ? xf[i].xaxis : "Time (ps)",
203 sprintf(buf, "%s.out", xf[i].filenm);
205 xf[i].fp2 = gmx_ffopen(buf, "w");
209 /* Read reference frame from tpx file to compute helix length */
210 snew(xref, top->atoms.nr);
211 read_tpx(ftp2fn(efTPX, NFILE, fnm),
212 NULL, NULL, &natoms, xref, NULL, NULL, NULL);
213 calc_hxprops(nres, bb, xref);
214 do_start_end(nres, bb, &nbb, bbindex, &nca, caindex, bRange, rStart, rEnd);
218 fprintf(stderr, "nca=%d, nbb=%d\n", nca, nbb);
219 pr_bb(stdout, nres, bb);
222 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
227 if ((teller++ % 10) == 0)
229 fprintf(stderr, "\rt=%.2f", t);
231 gmx_rmpbc(gpbc, natoms, box, x);
234 calc_hxprops(nres, bb, x);
237 do_start_end(nres, bb, &nbb, bbindex, &nca, caindex, FALSE, 0, 0);
242 rms = fit_ahx(nres, bb, natoms, nall, allindex, x, nca, caindex, bFit);
246 write_sto_conf(opt2fn("-cz", NFILE, fnm), "Helix fitted to Z-Axis",
247 &(top->atoms), x, NULL, ePBC, box);
250 xf[efhRAD].val = radius(xf[efhRAD].fp2, nca, caindex, x);
251 xf[efhTWIST].val = twist(nca, caindex, x);
252 xf[efhRISE].val = rise(nca, caindex, x);
253 xf[efhLEN].val = ahx_len(nca, caindex, x);
254 xf[efhCD222].val = ellipticity(nres, bb);
255 xf[efhDIP].val = dip(nbb, bbindex, x, top->atoms.atom);
256 xf[efhRMS].val = rms;
257 xf[efhCPHI].val = ca_phi(nca, caindex, x);
258 xf[efhPPRMS].val = pprms(xf[efhPPRMS].fp2, nres, bb);
260 for (j = 0; (j <= efhCPHI); j++)
262 fprintf(xf[j].fp, "%10g %10g\n", t, xf[j].val);
265 av_phipsi(xf[efhPHI].fp, xf[efhPSI].fp, xf[efhPHI].fp2, xf[efhPSI].fp2,
267 av_hblen(xf[efhHB3].fp, xf[efhHB3].fp2,
268 xf[efhHB4].fp, xf[efhHB4].fp2,
269 xf[efhHB5].fp, xf[efhHB5].fp2,
273 while (read_next_x(oenv, status, &t, x, box));
274 fprintf(stderr, "\n");
276 gmx_rmpbc_done(gpbc);
280 for (i = 0; (i < nres); i++)
284 fprintf(xf[efhRMSA].fp, "%10d %10g\n", r0+i, bb[i].rmsa/bb[i].nrms);
286 fprintf(xf[efhAHX].fp, "%10d %10g\n", r0+i, (bb[i].nhx*100.0)/(real )teller);
287 fprintf(xf[efhJCA].fp, "%10d %10g\n",
288 r0+i, 140.3+(bb[i].jcaha/(double)teller));
291 for (i = 0; (i < efhNR); i++)
293 gmx_ffclose(xf[i].fp);
296 gmx_ffclose(xf[i].fp2);
298 do_view(oenv, xf[i].filenm, "-nxy");