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,2015, 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/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"
61 int gmx_helix(int argc
, char *argv
[])
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.",
68 "to an ideal helix around the [IT]z[it]-axis and centered around the origin.",
69 "Then the following properties are computed:",
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 the number",
74 " 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",
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."
93 static gmx_bool bCheck
= FALSE
, bFit
= TRUE
, bDBG
= FALSE
, bEV
= FALSE
;
94 static int rStart
= 0, rEnd
= 0, r0
= 1;
96 { "-r0", FALSE
, etINT
, {&r0
},
97 "The first residue number in the sequence" },
98 { "-q", FALSE
, etBOOL
, {&bCheck
},
99 "Check at every step which part of the sequence is helical" },
100 { "-F", FALSE
, etBOOL
, {&bFit
},
101 "Toggle fit to a perfect helix" },
102 { "-db", FALSE
, etBOOL
, {&bDBG
},
103 "Print debug info" },
104 { "-ev", FALSE
, etBOOL
, {&bEV
},
105 "Write a new 'trajectory' file for ED" },
106 { "-ahxstart", FALSE
, etINT
, {&rStart
},
107 "First residue in helix" },
108 { "-ahxend", FALSE
, etINT
, {&rEnd
},
109 "Last residue in helix" }
122 t_xvgrfile xf
[efhNR
] = {
123 { NULL
, NULL
, TRUE
, "radius", "Helix radius", NULL
, "r (nm)", 0.0 },
124 { NULL
, NULL
, TRUE
, "twist", "Twist per residue", NULL
, "Angle (deg)", 0.0 },
125 { NULL
, NULL
, TRUE
, "rise", "Rise per residue", NULL
, "Rise (nm)", 0.0 },
126 { NULL
, NULL
, FALSE
, "len-ahx", "Length of the Helix", NULL
, "Length (nm)", 0.0 },
127 { NULL
, NULL
, FALSE
, "dip-ahx", "Helix Backbone Dipole", NULL
, "rq (nm e)", 0.0 },
128 { NULL
, NULL
, TRUE
, "rms-ahx", "RMS Deviation from Ideal Helix", NULL
, "RMS (nm)", 0.0 },
129 { NULL
, NULL
, FALSE
, "rmsa-ahx", "Average RMSD per Residue", "Residue", "RMS (nm)", 0.0 },
130 { NULL
, NULL
, FALSE
, "cd222", "Ellipticity at 222 nm", NULL
, "nm", 0.0 },
131 { NULL
, NULL
, TRUE
, "pprms", "RMS Distance from \\8a\\4-helix", NULL
, "deg", 0.0 },
132 { NULL
, NULL
, TRUE
, "caphi", "Average Ca-Ca Dihedral", NULL
, "\\8F\\4(deg)", 0.0 },
133 { NULL
, NULL
, TRUE
, "phi", "Average \\8F\\4 angles", NULL
, "deg", 0.0 },
134 { NULL
, NULL
, TRUE
, "psi", "Average \\8Y\\4 angles", NULL
, "deg", 0.0 },
135 { NULL
, NULL
, TRUE
, "hb3", "Average n-n+3 hbond length", NULL
, "nm", 0.0 },
136 { NULL
, NULL
, TRUE
, "hb4", "Average n-n+4 hbond length", NULL
, "nm", 0.0 },
137 { NULL
, NULL
, TRUE
, "hb5", "Average n-n+5 hbond length", NULL
, "nm", 0.0 },
138 { NULL
, NULL
, FALSE
, "JCaHa", "J-Coupling Values", "Residue", "Hz", 0.0 },
139 { NULL
, NULL
, FALSE
, "helicity", "Helicity per Residue", "Residue", "% of time", 0.0 }
142 gmx_output_env_t
*oenv
;
147 int i
, j
, nall
, nbb
, nca
, teller
;
148 int *bbindex
, *caindex
, *allindex
;
155 gmx_rmpbc_t gpbc
= NULL
;
158 { efTPR
, NULL
, NULL
, ffREAD
},
159 { efNDX
, NULL
, NULL
, ffREAD
},
160 { efTRX
, "-f", NULL
, ffREAD
},
161 { efSTO
, "-cz", "zconf", ffWRITE
},
163 #define NFILE asize(fnm)
165 if (!parse_common_args(&argc
, argv
, PCA_CAN_VIEW
| PCA_CAN_TIME
,
166 NFILE
, fnm
, asize(pa
), pa
, asize(desc
), desc
, 0, NULL
, &oenv
))
171 bRange
= (opt2parg_bSet("-ahxstart", asize(pa
), pa
) &&
172 opt2parg_bSet("-ahxend", asize(pa
), pa
));
174 top
= read_top(ftp2fn(efTPR
, NFILE
, fnm
), &ePBC
);
176 natoms
= read_first_x(oenv
, &status
, opt2fn("-f", NFILE
, fnm
), &t
, &x
, box
);
178 if (natoms
!= top
->atoms
.nr
)
180 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)",
181 top
->atoms
.nr
, natoms
);
184 bb
= mkbbind(ftp2fn(efNDX
, NFILE
, fnm
), &nres
, &nbb
, r0
, &nall
, &allindex
,
185 top
->atoms
.atomname
, top
->atoms
.atom
, top
->atoms
.resinfo
);
186 snew(bbindex
, natoms
);
189 fprintf(stderr
, "nall=%d\n", nall
);
191 /* Open output files, default x-axis is time */
192 for (i
= 0; (i
< efhNR
); i
++)
194 sprintf(buf
, "%s.xvg", xf
[i
].filenm
);
196 xf
[i
].fp
= xvgropen(buf
, xf
[i
].title
,
197 xf
[i
].xaxis
? xf
[i
].xaxis
: "Time (ps)",
201 sprintf(buf
, "%s.out", xf
[i
].filenm
);
203 xf
[i
].fp2
= gmx_ffopen(buf
, "w");
207 /* Read reference frame from tpx file to compute helix length */
208 snew(xref
, top
->atoms
.nr
);
209 read_tpx(ftp2fn(efTPR
, NFILE
, fnm
),
210 NULL
, NULL
, &natoms
, xref
, NULL
, NULL
);
211 calc_hxprops(nres
, bb
, xref
);
212 do_start_end(nres
, bb
, &nbb
, bbindex
, &nca
, caindex
, bRange
, rStart
, rEnd
);
216 fprintf(stderr
, "nca=%d, nbb=%d\n", nca
, nbb
);
217 pr_bb(stdout
, nres
, bb
);
220 gpbc
= gmx_rmpbc_init(&top
->idef
, ePBC
, natoms
);
225 if ((teller
++ % 10) == 0)
227 fprintf(stderr
, "\rt=%.2f", t
);
229 gmx_rmpbc(gpbc
, natoms
, box
, x
);
232 calc_hxprops(nres
, bb
, x
);
235 do_start_end(nres
, bb
, &nbb
, bbindex
, &nca
, caindex
, FALSE
, 0, 0);
240 rms
= fit_ahx(nres
, bb
, natoms
, nall
, allindex
, x
, nca
, caindex
, bFit
);
244 write_sto_conf(opt2fn("-cz", NFILE
, fnm
), "Helix fitted to Z-Axis",
245 &(top
->atoms
), x
, NULL
, ePBC
, box
);
248 xf
[efhRAD
].val
= radius(xf
[efhRAD
].fp2
, nca
, caindex
, x
);
249 xf
[efhTWIST
].val
= twist(nca
, caindex
, x
);
250 xf
[efhRISE
].val
= rise(nca
, caindex
, x
);
251 xf
[efhLEN
].val
= ahx_len(nca
, caindex
, x
);
252 xf
[efhCD222
].val
= ellipticity(nres
, bb
);
253 xf
[efhDIP
].val
= dip(nbb
, bbindex
, x
, top
->atoms
.atom
);
254 xf
[efhRMS
].val
= rms
;
255 xf
[efhCPHI
].val
= ca_phi(nca
, caindex
, x
);
256 xf
[efhPPRMS
].val
= pprms(xf
[efhPPRMS
].fp2
, nres
, bb
);
258 for (j
= 0; (j
<= efhCPHI
); j
++)
260 fprintf(xf
[j
].fp
, "%10g %10g\n", t
, xf
[j
].val
);
263 av_phipsi(xf
[efhPHI
].fp
, xf
[efhPSI
].fp
, xf
[efhPHI
].fp2
, xf
[efhPSI
].fp2
,
265 av_hblen(xf
[efhHB3
].fp
, xf
[efhHB3
].fp2
,
266 xf
[efhHB4
].fp
, xf
[efhHB4
].fp2
,
267 xf
[efhHB5
].fp
, xf
[efhHB5
].fp2
,
271 while (read_next_x(oenv
, status
, &t
, x
, box
));
272 fprintf(stderr
, "\n");
274 gmx_rmpbc_done(gpbc
);
278 for (i
= 0; (i
< nres
); i
++)
282 fprintf(xf
[efhRMSA
].fp
, "%10d %10g\n", r0
+i
, bb
[i
].rmsa
/bb
[i
].nrms
);
284 fprintf(xf
[efhAHX
].fp
, "%10d %10g\n", r0
+i
, (bb
[i
].nhx
*100.0)/static_cast<real
>(teller
));
285 fprintf(xf
[efhJCA
].fp
, "%10d %10g\n",
286 r0
+i
, 140.3+(bb
[i
].jcaha
/static_cast<double>(teller
)));
289 for (i
= 0; (i
< efhNR
); i
++)
294 xvgrclose(xf
[i
].fp2
);
296 do_view(oenv
, xf
[i
].filenm
, "-nxy");