Sort all includes in src/gromacs
[alexxy/gromacs.git] / src / gromacs / gmxana / 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  * 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.
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 <math.h>
40 #include <string.h>
41
42 #include "gromacs/commandline/pargs.h"
43 #include "gromacs/fileio/confio.h"
44 #include "gromacs/fileio/tpxio.h"
45 #include "gromacs/fileio/trxio.h"
46 #include "gromacs/fileio/xvgr.h"
47 #include "gromacs/gmxana/fitahx.h"
48 #include "gromacs/gmxana/gmx_ana.h"
49 #include "gromacs/gmxana/gstat.h"
50 #include "gromacs/gmxana/hxprops.h"
51 #include "gromacs/legacyheaders/macros.h"
52 #include "gromacs/legacyheaders/txtdump.h"
53 #include "gromacs/legacyheaders/typedefs.h"
54 #include "gromacs/legacyheaders/viewit.h"
55 #include "gromacs/math/utilities.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/pbcutil/rmpbc.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/futil.h"
60 #include "gromacs/utility/smalloc.h"
61
62 int gmx_helix(int argc, char *argv[])
63 {
64     const char        *desc[] = {
65         "[THISMODULE] computes all kinds of helix properties. First, the peptide",
66         "is checked to find the longest helical part, as determined by",
67         "hydrogen bonds and [GRK]phi[grk]/[GRK]psi[grk] angles.",
68         "That bit is fitted",
69         "to an ideal helix around the [IT]z[it]-axis and centered around the origin.",
70         "Then the following properties are computed:[PAR]",
71         "[BB]1.[bb] 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[BR]",
75         "[BB]2.[bb] 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.[BR]",
79         "[BB]3.[bb] 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[BR]",
82         "[BB]4.[bb] 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).[BR]",
86         "[BB]5.[bb] Helix dipole, backbone only (file [TT]dip-ahx.xvg[tt]).[BR]",
87         "[BB]6.[bb] RMS deviation from ideal helix, calculated for the C[GRK]alpha[grk]",
88         "atoms only (file [TT]rms-ahx.xvg[tt]).[BR]",
89         "[BB]7.[bb] Average C[GRK]alpha[grk] - C[GRK]alpha[grk] dihedral angle (file [TT]phi-ahx.xvg[tt]).[BR]",
90         "[BB]8.[bb] Average [GRK]phi[grk] and [GRK]psi[grk] angles (file [TT]phipsi.xvg[tt]).[BR]",
91         "[BB]9.[bb] Ellipticity at 222 nm according to Hirst and Brooks.",
92         "[PAR]"
93     };
94     static gmx_bool    bCheck = FALSE, bFit = TRUE, bDBG = FALSE, bEV = FALSE;
95     static int         rStart = 0, rEnd = 0, r0 = 1;
96     t_pargs            pa []  = {
97         { "-r0", FALSE, etINT, {&r0},
98           "The first residue number in the sequence" },
99         { "-q",  FALSE, etBOOL, {&bCheck},
100           "Check at every step which part of the sequence is helical" },
101         { "-F",  FALSE, etBOOL, {&bFit},
102           "Toggle fit to a perfect helix" },
103         { "-db", FALSE, etBOOL, {&bDBG},
104           "Print debug info" },
105         { "-ev", FALSE, etBOOL, {&bEV},
106           "Write a new 'trajectory' file for ED" },
107         { "-ahxstart", FALSE, etINT, {&rStart},
108           "First residue in helix" },
109         { "-ahxend", FALSE, etINT, {&rEnd},
110           "Last residue in helix" }
111     };
112
113     typedef struct {
114         FILE       *fp, *fp2;
115         gmx_bool    bfp2;
116         const char *filenm;
117         const char *title;
118         const char *xaxis;
119         const char *yaxis;
120         real        val;
121     } t_xvgrfile;
122
123     t_xvgrfile     xf[efhNR] = {
124         { NULL, NULL, TRUE,  "radius",  "Helix radius",               NULL, "r (nm)", 0.0 },
125         { NULL, NULL, TRUE,  "twist",   "Twist per residue",          NULL, "Angle (deg)", 0.0 },
126         { NULL, NULL, TRUE,  "rise",    "Rise per residue",           NULL, "Rise (nm)", 0.0 },
127         { NULL, NULL, FALSE, "len-ahx", "Length of the Helix",        NULL, "Length (nm)", 0.0 },
128         { NULL, NULL, FALSE, "dip-ahx", "Helix Backbone Dipole",      NULL, "rq (nm e)", 0.0 },
129         { NULL, NULL, TRUE,  "rms-ahx", "RMS Deviation from Ideal Helix", NULL, "RMS (nm)", 0.0 },
130         { NULL, NULL, FALSE, "rmsa-ahx", "Average RMSD per Residue",   "Residue", "RMS (nm)", 0.0 },
131         { NULL, NULL, FALSE,  "cd222",   "Ellipticity at 222 nm", NULL, "nm", 0.0 },
132         { NULL, NULL, TRUE,  "pprms",   "RMS Distance from \\8a\\4-helix", NULL, "deg", 0.0 },
133         { NULL, NULL, TRUE,  "caphi",   "Average Ca-Ca Dihedral",     NULL, "\\8F\\4(deg)", 0.0 },
134         { NULL, NULL, TRUE,  "phi",     "Average \\8F\\4 angles", NULL, "deg", 0.0 },
135         { NULL, NULL, TRUE,  "psi",     "Average \\8Y\\4 angles", NULL, "deg", 0.0 },
136         { NULL, NULL, TRUE,  "hb3",     "Average n-n+3 hbond length", NULL, "nm", 0.0 },
137         { NULL, NULL, TRUE,  "hb4",     "Average n-n+4 hbond length", NULL, "nm", 0.0 },
138         { NULL, NULL, TRUE,  "hb5",     "Average n-n+5 hbond length", NULL, "nm", 0.0 },
139         { NULL, NULL, FALSE,  "JCaHa",   "J-Coupling Values",        "Residue", "Hz", 0.0 },
140         { NULL, NULL, FALSE,  "helicity", "Helicity per Residue",     "Residue", "% of time", 0.0 }
141     };
142
143     output_env_t   oenv;
144     char           buf[54];
145     t_trxstatus   *status;
146     int            natoms, nre, nres;
147     t_bb          *bb;
148     int            i, j, m, nall, nbb, nca, teller, nSel = 0;
149     atom_id       *bbindex, *caindex, *allindex;
150     t_topology    *top;
151     int            ePBC;
152     rvec          *x, *xref;
153     real           t;
154     real           rms;
155     matrix         box;
156     gmx_rmpbc_t    gpbc = NULL;
157     gmx_bool       bRange;
158     t_filenm       fnm[] = {
159         { efTPR, NULL,  NULL,   ffREAD  },
160         { efNDX, NULL,  NULL,   ffREAD  },
161         { efTRX, "-f",  NULL,   ffREAD  },
162         { efSTO, "-cz", "zconf", ffWRITE },
163     };
164 #define NFILE asize(fnm)
165
166     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
167                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
168     {
169         return 0;
170     }
171
172     bRange = (opt2parg_bSet("-ahxstart", asize(pa), pa) &&
173               opt2parg_bSet("-ahxend", asize(pa), pa));
174
175     top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC);
176
177     natoms = read_first_x(oenv, &status, opt2fn("-f", NFILE, fnm), &t, &x, box);
178
179     if (natoms != top->atoms.nr)
180     {
181         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)",
182                   top->atoms.nr, natoms);
183     }
184
185     bb = mkbbind(ftp2fn(efNDX, NFILE, fnm), &nres, &nbb, r0, &nall, &allindex,
186                  top->atoms.atomname, top->atoms.atom, top->atoms.resinfo);
187     snew(bbindex, natoms);
188     snew(caindex, nres);
189
190     fprintf(stderr, "nall=%d\n", nall);
191
192     /* Open output files, default x-axis is time */
193     for (i = 0; (i < efhNR); i++)
194     {
195         sprintf(buf, "%s.xvg", xf[i].filenm);
196         remove(buf);
197         xf[i].fp = xvgropen(buf, xf[i].title,
198                             xf[i].xaxis ? xf[i].xaxis : "Time (ps)",
199                             xf[i].yaxis, oenv);
200         if (xf[i].bfp2)
201         {
202             sprintf(buf, "%s.out", xf[i].filenm);
203             remove(buf);
204             xf[i].fp2 = gmx_ffopen(buf, "w");
205         }
206     }
207
208     /* Read reference frame from tpx file to compute helix length */
209     snew(xref, top->atoms.nr);
210     read_tpx(ftp2fn(efTPR, NFILE, fnm),
211              NULL, NULL, &natoms, xref, NULL, NULL, NULL);
212     calc_hxprops(nres, bb, xref);
213     do_start_end(nres, bb, &nbb, bbindex, &nca, caindex, bRange, rStart, rEnd);
214     sfree(xref);
215     if (bDBG)
216     {
217         fprintf(stderr, "nca=%d, nbb=%d\n", nca, nbb);
218         pr_bb(stdout, nres, bb);
219     }
220
221     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
222
223     teller = 0;
224     do
225     {
226         if ((teller++ % 10) == 0)
227         {
228             fprintf(stderr, "\rt=%.2f", t);
229         }
230         gmx_rmpbc(gpbc, natoms, box, x);
231
232
233         calc_hxprops(nres, bb, x);
234         if (bCheck)
235         {
236             do_start_end(nres, bb, &nbb, bbindex, &nca, caindex, FALSE, 0, 0);
237         }
238
239         if (nca >= 5)
240         {
241             rms = fit_ahx(nres, bb, natoms, nall, allindex, x, nca, caindex, bFit);
242
243             if (teller == 1)
244             {
245                 write_sto_conf(opt2fn("-cz", NFILE, fnm), "Helix fitted to Z-Axis",
246                                &(top->atoms), x, NULL, ePBC, box);
247             }
248
249             xf[efhRAD].val   = radius(xf[efhRAD].fp2, nca, caindex, x);
250             xf[efhTWIST].val = twist(nca, caindex, x);
251             xf[efhRISE].val  = rise(nca, caindex, x);
252             xf[efhLEN].val   = ahx_len(nca, caindex, x);
253             xf[efhCD222].val = ellipticity(nres, bb);
254             xf[efhDIP].val   = dip(nbb, bbindex, x, top->atoms.atom);
255             xf[efhRMS].val   = rms;
256             xf[efhCPHI].val  = ca_phi(nca, caindex, x);
257             xf[efhPPRMS].val = pprms(xf[efhPPRMS].fp2, nres, bb);
258
259             for (j = 0; (j <= efhCPHI); j++)
260             {
261                 fprintf(xf[j].fp,   "%10g  %10g\n", t, xf[j].val);
262             }
263
264             av_phipsi(xf[efhPHI].fp, xf[efhPSI].fp, xf[efhPHI].fp2, xf[efhPSI].fp2,
265                       t, nres, bb);
266             av_hblen(xf[efhHB3].fp, xf[efhHB3].fp2,
267                      xf[efhHB4].fp, xf[efhHB4].fp2,
268                      xf[efhHB5].fp, xf[efhHB5].fp2,
269                      t, nres, bb);
270         }
271     }
272     while (read_next_x(oenv, status, &t, x, box));
273     fprintf(stderr, "\n");
274
275     gmx_rmpbc_done(gpbc);
276
277     close_trj(status);
278
279     for (i = 0; (i < nres); i++)
280     {
281         if (bb[i].nrms > 0)
282         {
283             fprintf(xf[efhRMSA].fp, "%10d  %10g\n", r0+i, bb[i].rmsa/bb[i].nrms);
284         }
285         fprintf(xf[efhAHX].fp, "%10d  %10g\n", r0+i, (bb[i].nhx*100.0)/(real )teller);
286         fprintf(xf[efhJCA].fp, "%10d  %10g\n",
287                 r0+i, 140.3+(bb[i].jcaha/(double)teller));
288     }
289
290     for (i = 0; (i < efhNR); i++)
291     {
292         gmx_ffclose(xf[i].fp);
293         if (xf[i].bfp2)
294         {
295             gmx_ffclose(xf[i].fp2);
296         }
297         do_view(oenv, xf[i].filenm, "-nxy");
298     }
299
300     return 0;
301 }