Move gmx_fatal.* to utility/fatalerror.*
[alexxy/gromacs.git] / src / gromacs / gmxana / gmx_bundle.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 #ifdef HAVE_CONFIG_H
38 #include <config.h>
39 #endif
40 #include <math.h>
41 #include <string.h>
42
43 #include "gromacs/commandline/pargs.h"
44 #include "sysstuff.h"
45 #include "typedefs.h"
46 #include "gromacs/utility/smalloc.h"
47 #include "macros.h"
48 #include "vec.h"
49 #include "gromacs/fileio/futil.h"
50 #include "index.h"
51 #include "xvgr.h"
52 #include "rmpbc.h"
53 #include "gromacs/fileio/tpxio.h"
54 #include "gromacs/fileio/trxio.h"
55 #include "physics.h"
56 #include "gmx_ana.h"
57
58 #include "gromacs/utility/fatalerror.h"
59
60 #define MAX_ENDS 3
61
62 typedef struct {
63     int   n;
64     int   nend;
65     rvec *end[MAX_ENDS];
66     rvec *mid;
67     rvec *dir;
68     real *len;
69 } t_bundle;
70
71 static void rotate_ends(t_bundle *bun, rvec axis, int c0, int c1)
72 {
73     int  end, i;
74     rvec ax, tmp;
75
76     unitv(axis, ax);
77     for (end = 0; end < bun->nend; end++)
78     {
79         for (i = 0; i < bun->n; i++)
80         {
81             copy_rvec(bun->end[end][i], tmp);
82             bun->end[end][i][c0] = ax[c1]*tmp[c0] - ax[c0]*tmp[c1];
83             bun->end[end][i][c1] = ax[c0]*tmp[c0] + ax[c1]*tmp[c1];
84         }
85     }
86     copy_rvec(axis, tmp);
87     axis[c0] = ax[c1]*tmp[c0] - ax[c0]*tmp[c1];
88     axis[c1] = ax[c0]*tmp[c0] + ax[c1]*tmp[c1];
89 }
90
91 static void calc_axes(rvec x[], t_atom atom[], int gnx[], atom_id *index[],
92                       gmx_bool bRot, t_bundle *bun)
93 {
94     int   end, i, div, d;
95     real *mtot, m;
96     rvec  axis[MAX_ENDS], cent;
97
98     snew(mtot, bun->n);
99
100     for (end = 0; end < bun->nend; end++)
101     {
102         for (i = 0; i < bun->n; i++)
103         {
104             clear_rvec(bun->end[end][i]);
105             mtot[i] = 0;
106         }
107         div = gnx[end]/bun->n;
108         for (i = 0; i < gnx[end]; i++)
109         {
110             m = atom[index[end][i]].m;
111             for (d = 0; d < DIM; d++)
112             {
113                 bun->end[end][i/div][d] += m*x[index[end][i]][d];
114             }
115             mtot[i/div] += m;
116         }
117         clear_rvec(axis[end]);
118         for (i = 0; i < bun->n; i++)
119         {
120             svmul(1.0/mtot[i], bun->end[end][i], bun->end[end][i]);
121             rvec_inc(axis[end], bun->end[end][i]);
122         }
123         svmul(1.0/bun->n, axis[end], axis[end]);
124     }
125     sfree(mtot);
126
127     rvec_add(axis[0], axis[1], cent);
128     svmul(0.5, cent, cent);
129     /* center the bundle on the origin */
130     for (end = 0; end < bun->nend; end++)
131     {
132         rvec_dec(axis[end], cent);
133         for (i = 0; i < bun->n; i++)
134         {
135             rvec_dec(bun->end[end][i], cent);
136         }
137     }
138     if (bRot)
139     {
140         /* rotate the axis parallel to the z-axis */
141         rotate_ends(bun, axis[0], YY, ZZ);
142         rotate_ends(bun, axis[0], XX, ZZ);
143     }
144     for (i = 0; i < bun->n; i++)
145     {
146         rvec_add(bun->end[0][i], bun->end[1][i], bun->mid[i]);
147         svmul(0.5, bun->mid[i], bun->mid[i]);
148         rvec_sub(bun->end[0][i], bun->end[1][i], bun->dir[i]);
149         bun->len[i] = norm(bun->dir[i]);
150         unitv(bun->dir[i], bun->dir[i]);
151     }
152 }
153
154 static void dump_axes(t_trxstatus *status, t_trxframe *fr, t_atoms *outat,
155                       t_bundle *bun)
156 {
157     t_trxframe   frout;
158     static rvec *xout = NULL;
159     int          i, end;
160
161     if (xout == NULL)
162     {
163         snew(xout, outat->nr);
164     }
165
166     for (i = 0; i < bun->n; i++)
167     {
168         copy_rvec(bun->end[0][i], xout[3*i]);
169         if (bun->nend >= 3)
170         {
171             copy_rvec(bun->end[2][i], xout[3*i+1]);
172         }
173         else
174         {
175             copy_rvec(bun->mid[i], xout[3*i+1]);
176         }
177         copy_rvec(bun->end[1][i], xout[3*i+2]);
178     }
179     frout        = *fr;
180     frout.bV     = FALSE;
181     frout.bF     = FALSE;
182     frout.bBox   = FALSE;
183     frout.bAtoms = TRUE;
184     frout.natoms = outat->nr;
185     frout.atoms  = outat;
186     frout.x      = xout;
187     write_trxframe(status, &frout, NULL);
188 }
189
190 int gmx_bundle(int argc, char *argv[])
191 {
192     const char     *desc[] = {
193         "[THISMODULE] analyzes bundles of axes. The axes can be for instance",
194         "helix axes. The program reads two index groups and divides both",
195         "of them in [TT]-na[tt] parts. The centers of mass of these parts",
196         "define the tops and bottoms of the axes.",
197         "Several quantities are written to file:",
198         "the axis length, the distance and the z-shift of the axis mid-points",
199         "with respect to the average center of all axes, the total tilt,",
200         "the radial tilt and the lateral tilt with respect to the average axis.",
201         "[PAR]",
202         "With options [TT]-ok[tt], [TT]-okr[tt] and [TT]-okl[tt] the total,",
203         "radial and lateral kinks of the axes are plotted. An extra index",
204         "group of kink atoms is required, which is also divided into [TT]-na[tt]",
205         "parts. The kink angle is defined as the angle between the kink-top and",
206         "the bottom-kink vectors.",
207         "[PAR]",
208         "With option [TT]-oa[tt] the top, mid (or kink when [TT]-ok[tt] is set)",
209         "and bottom points of each axis",
210         "are written to a [TT].pdb[tt] file each frame. The residue numbers correspond",
211         "to the axis numbers. When viewing this file with Rasmol, use the",
212         "command line option [TT]-nmrpdb[tt], and type [TT]set axis true[tt] to",
213         "display the reference axis."
214     };
215     static int      n    = 0;
216     static gmx_bool bZ   = FALSE;
217     t_pargs         pa[] = {
218         { "-na", FALSE, etINT, {&n},
219           "Number of axes" },
220         { "-z", FALSE, etBOOL, {&bZ},
221           "Use the [IT]z[it]-axis as reference instead of the average axis" }
222     };
223     FILE           *out, *flen, *fdist, *fz, *ftilt, *ftiltr, *ftiltl;
224     FILE           *fkink = NULL, *fkinkr = NULL, *fkinkl = NULL;
225     t_trxstatus    *status;
226     t_trxstatus    *fpdb;
227     t_topology      top;
228     int             ePBC;
229     rvec           *xtop;
230     matrix          box;
231     t_trxframe      fr;
232     t_atoms         outatoms;
233     real            t, comp;
234     int             natoms;
235     char           *grpname[MAX_ENDS], title[256];
236     /* FIXME: The constness should not be cast away */
237     char           *anm = (char *)"CA", *rnm = (char *)"GLY";
238     int             i, j, gnx[MAX_ENDS];
239     atom_id        *index[MAX_ENDS];
240     t_bundle        bun;
241     gmx_bool        bKink;
242     rvec            va, vb, vc, vr, vl;
243     output_env_t    oenv;
244     gmx_rmpbc_t     gpbc = NULL;
245
246 #define NLEG asize(leg)
247     t_filenm fnm[] = {
248         { efTRX, "-f", NULL, ffREAD },
249         { efTPS, NULL, NULL, ffREAD },
250         { efNDX, NULL, NULL, ffOPTRD },
251         { efXVG, "-ol", "bun_len", ffWRITE },
252         { efXVG, "-od", "bun_dist", ffWRITE },
253         { efXVG, "-oz", "bun_z", ffWRITE },
254         { efXVG, "-ot", "bun_tilt", ffWRITE },
255         { efXVG, "-otr", "bun_tiltr", ffWRITE },
256         { efXVG, "-otl", "bun_tiltl", ffWRITE },
257         { efXVG, "-ok", "bun_kink", ffOPTWR },
258         { efXVG, "-okr", "bun_kinkr", ffOPTWR },
259         { efXVG, "-okl", "bun_kinkl", ffOPTWR },
260         { efPDB, "-oa", "axes", ffOPTWR }
261     };
262 #define NFILE asize(fnm)
263
264     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
265                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
266     {
267         return 0;
268     }
269
270     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, NULL, box, TRUE);
271
272     bKink = opt2bSet("-ok", NFILE, fnm) || opt2bSet("-okr", NFILE, fnm)
273         || opt2bSet("-okl", NFILE, fnm);
274     if (bKink)
275     {
276         bun.nend = 3;
277     }
278     else
279     {
280         bun.nend = 2;
281     }
282
283     fprintf(stderr, "Select a group of top and a group of bottom ");
284     if (bKink)
285     {
286         fprintf(stderr, "and a group of kink ");
287     }
288     fprintf(stderr, "atoms\n");
289     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), bun.nend,
290               gnx, index, grpname);
291
292     if (n <= 0 || gnx[0] % n || gnx[1] % n || (bKink && gnx[2] % n))
293     {
294         gmx_fatal(FARGS,
295                   "The size of one of your index groups is not a multiple of n");
296     }
297     bun.n = n;
298     snew(bun.end[0], n);
299     snew(bun.end[1], n);
300     if (bKink)
301     {
302         snew(bun.end[2], n);
303     }
304     snew(bun.mid, n);
305     snew(bun.dir, n);
306     snew(bun.len, n);
307
308     flen   = xvgropen(opt2fn("-ol", NFILE, fnm), "Axis lengths",
309                       output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
310     fdist  = xvgropen(opt2fn("-od", NFILE, fnm), "Distance of axis centers",
311                       output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
312     fz     = xvgropen(opt2fn("-oz", NFILE, fnm), "Z-shift of axis centers",
313                       output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
314     ftilt  = xvgropen(opt2fn("-ot", NFILE, fnm), "Axis tilts",
315                       output_env_get_xvgr_tlabel(oenv), "(degrees)", oenv);
316     ftiltr = xvgropen(opt2fn("-otr", NFILE, fnm), "Radial axis tilts",
317                       output_env_get_xvgr_tlabel(oenv), "(degrees)", oenv);
318     ftiltl = xvgropen(opt2fn("-otl", NFILE, fnm), "Lateral axis tilts",
319                       output_env_get_xvgr_tlabel(oenv), "(degrees)", oenv);
320
321     if (bKink)
322     {
323         fkink  = xvgropen(opt2fn("-ok", NFILE, fnm), "Kink angles",
324                           output_env_get_xvgr_tlabel(oenv), "(degrees)", oenv);
325         fkinkr = xvgropen(opt2fn("-okr", NFILE, fnm), "Radial kink angles",
326                           output_env_get_xvgr_tlabel(oenv), "(degrees)", oenv);
327         if (output_env_get_print_xvgr_codes(oenv))
328         {
329             fprintf(fkinkr, "@ subtitle \"+ = ) (   - = ( )\"\n");
330         }
331         fkinkl = xvgropen(opt2fn("-okl", NFILE, fnm), "Lateral kink angles",
332                           output_env_get_xvgr_tlabel(oenv), "(degrees)", oenv);
333     }
334
335     if (opt2bSet("-oa", NFILE, fnm))
336     {
337         init_t_atoms(&outatoms, 3*n, FALSE);
338         outatoms.nr = 3*n;
339         for (i = 0; i < 3*n; i++)
340         {
341             outatoms.atomname[i]       = &anm;
342             outatoms.atom[i].resind    = i/3;
343             outatoms.resinfo[i/3].name = &rnm;
344             outatoms.resinfo[i/3].nr   = i/3 + 1;
345             outatoms.resinfo[i/3].ic   = ' ';
346         }
347         fpdb = open_trx(opt2fn("-oa", NFILE, fnm), "w");
348     }
349     else
350     {
351         fpdb = NULL;
352     }
353
354     read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, TRX_NEED_X);
355     gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);
356
357     do
358     {
359         gmx_rmpbc_trxfr(gpbc, &fr);
360         calc_axes(fr.x, top.atoms.atom, gnx, index, !bZ, &bun);
361         t = output_env_conv_time(oenv, fr.time);
362         fprintf(flen, " %10g", t);
363         fprintf(fdist, " %10g", t);
364         fprintf(fz, " %10g", t);
365         fprintf(ftilt, " %10g", t);
366         fprintf(ftiltr, " %10g", t);
367         fprintf(ftiltl, " %10g", t);
368         if (bKink)
369         {
370             fprintf(fkink, " %10g", t);
371             fprintf(fkinkr, " %10g", t);
372             fprintf(fkinkl, " %10g", t);
373         }
374
375         for (i = 0; i < bun.n; i++)
376         {
377             fprintf(flen, " %6g", bun.len[i]);
378             fprintf(fdist, " %6g", norm(bun.mid[i]));
379             fprintf(fz, " %6g", bun.mid[i][ZZ]);
380             fprintf(ftilt, " %6g", RAD2DEG*acos(bun.dir[i][ZZ]));
381             comp = bun.mid[i][XX]*bun.dir[i][XX]+bun.mid[i][YY]*bun.dir[i][YY];
382             fprintf(ftiltr, " %6g", RAD2DEG*
383                     asin(comp/sqrt(sqr(comp)+sqr(bun.dir[i][ZZ]))));
384             comp = bun.mid[i][YY]*bun.dir[i][XX]-bun.mid[i][XX]*bun.dir[i][YY];
385             fprintf(ftiltl, " %6g", RAD2DEG*
386                     asin(comp/sqrt(sqr(comp)+sqr(bun.dir[i][ZZ]))));
387             if (bKink)
388             {
389                 rvec_sub(bun.end[0][i], bun.end[2][i], va);
390                 rvec_sub(bun.end[2][i], bun.end[1][i], vb);
391                 unitv_no_table(va, va);
392                 unitv_no_table(vb, vb);
393                 fprintf(fkink, " %6g", RAD2DEG*acos(iprod(va, vb)));
394                 cprod(va, vb, vc);
395                 copy_rvec(bun.mid[i], vr);
396                 vr[ZZ] = 0;
397                 unitv_no_table(vr, vr);
398                 fprintf(fkinkr, " %6g", RAD2DEG*asin(iprod(vc, vr)));
399                 vl[XX] = vr[YY];
400                 vl[YY] = -vr[XX];
401                 vl[ZZ] = 0;
402                 fprintf(fkinkl, " %6g", RAD2DEG*asin(iprod(vc, vl)));
403             }
404         }
405         fprintf(flen, "\n");
406         fprintf(fdist, "\n");
407         fprintf(fz, "\n");
408         fprintf(ftilt, "\n");
409         fprintf(ftiltr, "\n");
410         fprintf(ftiltl, "\n");
411         if (bKink)
412         {
413             fprintf(fkink, "\n");
414             fprintf(fkinkr, "\n");
415             fprintf(fkinkl, "\n");
416         }
417         if (fpdb)
418         {
419             dump_axes(fpdb, &fr, &outatoms, &bun);
420         }
421     }
422     while (read_next_frame(oenv, status, &fr));
423     gmx_rmpbc_done(gpbc);
424
425     close_trx(status);
426
427     if (fpdb)
428     {
429         close_trx(fpdb);
430     }
431     gmx_ffclose(flen);
432     gmx_ffclose(fdist);
433     gmx_ffclose(fz);
434     gmx_ffclose(ftilt);
435     gmx_ffclose(ftiltr);
436     gmx_ffclose(ftiltl);
437     if (bKink)
438     {
439         gmx_ffclose(fkink);
440         gmx_ffclose(fkinkr);
441         gmx_ffclose(fkinkl);
442     }
443
444     return 0;
445 }