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/commandline/pargs.h"
43 #include "gromacs/fileio/trnio.h"
44 #include "gromacs/fileio/xvgr.h"
45 #include "gromacs/gmxana/gmx_ana.h"
46 #include "gromacs/gmxana/gstat.h"
47 #include "gromacs/legacyheaders/copyrite.h"
48 #include "gromacs/legacyheaders/macros.h"
49 #include "gromacs/legacyheaders/typedefs.h"
50 #include "gromacs/legacyheaders/viewit.h"
51 #include "gromacs/math/units.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/topology/index.h"
54 #include "gromacs/utility/fatalerror.h"
55 #include "gromacs/utility/futil.h"
56 #include "gromacs/utility/smalloc.h"
59 static void dump_dih_trn(int nframes, int nangles, real **dih, const char *fn,
62 int i, j, k, l, m, na;
65 matrix box = {{2, 0, 0}, {0, 2, 0}, {0, 0, 2}};
76 printf("There are %d dihedrals. Will fill %d atom positions with cos/sin\n",
79 trn = open_trn(fn, "w");
80 for (i = 0; (i < nframes); i++)
83 for (j = 0; (j < nangles); j++)
85 for (m = 0; (m < 2); m++)
87 x[k][l] = (m == 0) ? cos(dih[j][i]) : sin(dih[j][i]);
96 fwrite_trn(trn, i, time[i], 0, box, na, x, NULL, NULL);
102 int gmx_g_angle(int argc, char *argv[])
104 static const char *desc[] = {
105 "[THISMODULE] computes the angle distribution for a number of angles",
106 "or dihedrals.[PAR]",
107 "With option [TT]-ov[tt], you can plot the average angle of",
108 "a group of angles as a function of time. With the [TT]-all[tt] option,",
109 "the first graph is the average and the rest are the individual angles.[PAR]",
110 "With the [TT]-of[tt] option, [THISMODULE] also calculates the fraction of trans",
111 "dihedrals (only for dihedrals) as function of time, but this is",
112 "probably only fun for a select few.[PAR]",
113 "With option [TT]-oc[tt], a dihedral correlation function is calculated.[PAR]",
114 "It should be noted that the index file must contain",
115 "atom triplets for angles or atom quadruplets for dihedrals.",
116 "If this is not the case, the program will crash.[PAR]",
117 "With option [TT]-or[tt], a trajectory file is dumped containing cos and",
118 "sin of selected dihedral angles, which subsequently can be used as",
119 "input for a principal components analysis using [gmx-covar].[PAR]",
120 "Option [TT]-ot[tt] plots when transitions occur between",
121 "dihedral rotamers of multiplicity 3 and [TT]-oh[tt]",
122 "records a histogram of the times between such transitions,",
123 "assuming the input trajectory frames are equally spaced in time."
125 static const char *opt[] = { NULL, "angle", "dihedral", "improper", "ryckaert-bellemans", NULL };
126 static gmx_bool bALL = FALSE, bChandler = FALSE, bAverCorr = FALSE, bPBC = TRUE;
127 static real binwidth = 1;
129 { "-type", FALSE, etENUM, {opt},
130 "Type of angle to analyse" },
131 { "-all", FALSE, etBOOL, {&bALL},
132 "Plot all angles separately in the averages file, in the order of appearance in the index file." },
133 { "-binwidth", FALSE, etREAL, {&binwidth},
134 "binwidth (degrees) for calculating the distribution" },
135 { "-periodic", FALSE, etBOOL, {&bPBC},
136 "Print dihedral angles modulo 360 degrees" },
137 { "-chandler", FALSE, etBOOL, {&bChandler},
138 "Use Chandler correlation function (N[trans] = 1, N[gauche] = 0) rather than cosine correlation function. Trans is defined as phi < -60 or phi > 60." },
139 { "-avercorr", FALSE, etBOOL, {&bAverCorr},
140 "Average the correlation functions for the individual angles/dihedrals" }
142 static const char *bugs[] = {
143 "Counting transitions only works for dihedrals with multiplicity 3"
151 real maxang, Jc, S2, norm_fac, maxstat;
153 int nframes, maxangstat, mult, *angstat;
154 int i, j, total, nangles, natoms, nat2, first, last, angind;
155 gmx_bool bAver, bRb, bPeriodic,
156 bFrac, /* calculate fraction too? */
157 bTrans, /* worry about transtions too? */
158 bCorr; /* correlation function ? */
159 real t, aa, aver, aver2, aversig, fraction; /* fraction trans dihedrals */
162 real **dih = NULL; /* mega array with all dih. angles at all times*/
164 real *time, *trans_frac, *aver_angle;
166 { efTRX, "-f", NULL, ffREAD },
167 { efNDX, NULL, "angle", ffREAD },
168 { efXVG, "-od", "angdist", ffWRITE },
169 { efXVG, "-ov", "angaver", ffOPTWR },
170 { efXVG, "-of", "dihfrac", ffOPTWR },
171 { efXVG, "-ot", "dihtrans", ffOPTWR },
172 { efXVG, "-oh", "trhisto", ffOPTWR },
173 { efXVG, "-oc", "dihcorr", ffOPTWR },
174 { efTRR, "-or", NULL, ffOPTWR }
176 #define NFILE asize(fnm)
182 ppa = add_acf_pargs(&npargs, pa);
183 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
184 NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
208 if (opt2bSet("-or", NFILE, fnm))
212 gmx_fatal(FARGS, "Can not combine angles with trn dump");
216 please_cite(stdout, "Mu2005a");
220 /* Calculate bin size */
221 maxangstat = (int)(maxang/binwidth+0.5);
222 binwidth = maxang/maxangstat;
224 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
225 nangles = isize/mult;
226 if ((isize % mult) != 0)
228 gmx_fatal(FARGS, "number of index elements not multiple of %d, "
229 "these can not be %s\n",
230 mult, (mult == 3) ? "angle triplets" : "dihedral quadruplets");
234 /* Check whether specific analysis has to be performed */
235 bCorr = opt2bSet("-oc", NFILE, fnm);
236 bAver = opt2bSet("-ov", NFILE, fnm);
237 bTrans = opt2bSet("-ot", NFILE, fnm);
238 bFrac = opt2bSet("-of", NFILE, fnm);
239 if (bTrans && opt[0][0] != 'd')
241 fprintf(stderr, "Option -ot should only accompany -type dihedral. Disabling -ot.\n");
245 if (bChandler && !bCorr)
252 fprintf(stderr, "Warning:"
253 " calculating fractions as defined in this program\n"
254 "makes sense for Ryckaert Bellemans dihs. only. Ignoring -of\n\n");
258 if ( (bTrans || bFrac || bCorr) && mult == 3)
260 gmx_fatal(FARGS, "Can only do transition, fraction or correlation\n"
261 "on dihedrals. Select -d\n");
265 * We need to know the nr of frames so we can allocate memory for an array
266 * with all dihedral angles at all timesteps. Works for me.
268 if (bTrans || bCorr || bALL || opt2bSet("-or", NFILE, fnm))
273 snew(angstat, maxangstat);
275 read_ang_dih(ftp2fn(efTRX, NFILE, fnm), (mult == 3),
276 bALL || bCorr || bTrans || opt2bSet("-or", NFILE, fnm),
277 bRb, bPBC, maxangstat, angstat,
278 &nframes, &time, isize, index, &trans_frac, &aver_angle, dih,
281 dt = (time[nframes-1]-time[0])/(nframes-1);
285 sprintf(title, "Average Angle: %s", grpname);
286 out = xvgropen(opt2fn("-ov", NFILE, fnm),
287 title, "Time (ps)", "Angle (degrees)", oenv);
288 for (i = 0; (i < nframes); i++)
290 fprintf(out, "%10.5f %8.3f", time[i], aver_angle[i]*RAD2DEG);
293 for (j = 0; (j < nangles); j++)
298 fprintf(out, " %8.3f", atan2(sin(dd), cos(dd))*RAD2DEG);
302 fprintf(out, " %8.3f", dih[j][i]*RAD2DEG);
310 if (opt2bSet("-or", NFILE, fnm))
312 dump_dih_trn(nframes, nangles, dih, opt2fn("-or", NFILE, fnm), time);
317 sprintf(title, "Trans fraction: %s", grpname);
318 out = xvgropen(opt2fn("-of", NFILE, fnm),
319 title, "Time (ps)", "Fraction", oenv);
321 for (i = 0; (i < nframes); i++)
323 fprintf(out, "%10.5f %10.3f\n", time[i], trans_frac[i]);
324 tfrac += trans_frac[i];
329 fprintf(stderr, "Average trans fraction: %g\n", tfrac);
335 ana_dih_trans(opt2fn("-ot", NFILE, fnm), opt2fn("-oh", NFILE, fnm),
336 dih, nframes, nangles, grpname, time, bRb, oenv);
341 /* Autocorrelation function */
344 fprintf(stderr, "Not enough frames for correlation function\n");
351 real dval, sixty = DEG2RAD*60;
354 for (i = 0; (i < nangles); i++)
356 for (j = 0; (j < nframes); j++)
361 bTest = (dval > -sixty) && (dval < sixty);
365 bTest = (dval < -sixty) || (dval > sixty);
369 dih[i][j] = dval-tfrac;
386 do_autocorr(opt2fn("-oc", NFILE, fnm), oenv,
387 "Dihedral Autocorrelation Function",
388 nframes, nangles, dih, dt, mode, bAverCorr);
393 /* Determine the non-zero part of the distribution */
394 for (first = 0; (first < maxangstat-1) && (angstat[first+1] == 0); first++)
398 for (last = maxangstat-1; (last > 0) && (angstat[last-1] == 0); last--)
404 for (i = 0; (i < nframes); i++)
406 aver += RAD2DEG*aver_angle[i];
407 aver2 += sqr(RAD2DEG*aver_angle[i]);
409 aver /= (real) nframes;
410 aver2 /= (real) nframes;
411 aversig = sqrt(aver2-sqr(aver));
412 printf("Found points in the range from %d to %d (max %d)\n",
413 first, last, maxangstat);
414 printf(" < angle > = %g\n", aver);
415 printf("< angle^2 > = %g\n", aver2);
416 printf("Std. Dev. = %g\n", aversig);
420 sprintf(title, "Angle Distribution: %s", grpname);
424 sprintf(title, "Dihedral Distribution: %s", grpname);
426 calc_distribution_props(maxangstat, angstat, -180.0, 0, NULL, &S2);
427 fprintf(stderr, "Order parameter S^2 = %g\n", S2);
430 bPeriodic = (mult == 4) && (first == 0) && (last == maxangstat-1);
432 out = xvgropen(opt2fn("-od", NFILE, fnm), title, "Degrees", "", oenv);
433 if (output_env_get_print_xvgr_codes(oenv))
435 fprintf(out, "@ subtitle \"average angle: %g\\So\\N\"\n", aver);
437 norm_fac = 1.0/(nangles*nframes*binwidth);
441 for (i = first; (i <= last); i++)
443 maxstat = max(maxstat, angstat[i]*norm_fac);
445 if (output_env_get_print_xvgr_codes(oenv))
447 fprintf(out, "@with g0\n");
448 fprintf(out, "@ world xmin -180\n");
449 fprintf(out, "@ world xmax 180\n");
450 fprintf(out, "@ world ymin 0\n");
451 fprintf(out, "@ world ymax %g\n", maxstat*1.05);
452 fprintf(out, "@ xaxis tick major 60\n");
453 fprintf(out, "@ xaxis tick minor 30\n");
454 fprintf(out, "@ yaxis tick major 0.005\n");
455 fprintf(out, "@ yaxis tick minor 0.0025\n");
458 for (i = first; (i <= last); i++)
460 fprintf(out, "%10g %10f\n", i*binwidth+180.0-maxang, angstat[i]*norm_fac);
464 /* print first bin again as last one */
465 fprintf(out, "%10g %10f\n", 180.0, angstat[0]*norm_fac);
470 do_view(oenv, opt2fn("-od", NFILE, fnm), "-nxy");
473 do_view(oenv, opt2fn("-ov", NFILE, fnm), "-nxy");