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,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/commandline/viewit.h"
47 #include "gromacs/correlationfunctions/autocorr.h"
48 #include "gromacs/fileio/trrio.h"
49 #include "gromacs/fileio/xvgr.h"
50 #include "gromacs/gmxana/angle_correction.h"
51 #include "gromacs/gmxana/gmx_ana.h"
52 #include "gromacs/gmxana/gstat.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/topology/index.h"
57 #include "gromacs/utility/arraysize.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/futil.h"
60 #include "gromacs/utility/gmxassert.h"
61 #include "gromacs/utility/pleasecite.h"
62 #include "gromacs/utility/smalloc.h"
64 static void dump_dih_trr(int nframes, int nangles, real** dih, const char* fn, real* time)
66 int i, j, k, l, m, na;
69 matrix box = { { 2, 0, 0 }, { 0, 2, 0 }, { 0, 0, 2 } };
80 printf("There are %d dihedrals. Will fill %d atom positions with cos/sin\n", nangles, na);
82 fio = gmx_trr_open(fn, "w");
83 for (i = 0; (i < nframes); i++)
86 for (j = 0; (j < nangles); j++)
88 for (m = 0; (m < 2); m++)
90 // This is just because the compler and static-analyzer cannot
91 // know that dih[j][i] is always valid. Since it occurs in the innermost
92 // loop over angles and will only trigger on coding errors, we
93 // only enable it for debug builds.
94 GMX_ASSERT(dih != nullptr && dih[j] != nullptr, "Incorrect dihedral array data");
95 x[k][l] = (m == 0) ? std::cos(dih[j][i]) : std::sin(dih[j][i]);
104 gmx_trr_write_frame(fio, i, time[i], 0, box, na, x, nullptr, nullptr);
110 int gmx_g_angle(int argc, char* argv[])
112 static const char* desc[] = {
113 "[THISMODULE] computes the angle distribution for a number of angles",
114 "or dihedrals.[PAR]",
115 "With option [TT]-ov[tt], you can plot the average angle of",
116 "a group of angles as a function of time. With the [TT]-all[tt] option,",
117 "the first graph is the average and the rest are the individual angles.[PAR]",
118 "With the [TT]-of[tt] option, [THISMODULE] also calculates the fraction of trans",
119 "dihedrals (only for dihedrals) as function of time, but this is",
120 "probably only fun for a select few.[PAR]",
121 "With option [TT]-oc[tt], a dihedral correlation function is calculated.[PAR]",
122 "It should be noted that the index file must contain",
123 "atom triplets for angles or atom quadruplets for dihedrals.",
124 "If this is not the case, the program will crash.[PAR]",
125 "With option [TT]-or[tt], a trajectory file is dumped containing cos and",
126 "sin of selected dihedral angles, which subsequently can be used as",
127 "input for a principal components analysis using [gmx-covar].[PAR]",
128 "Option [TT]-ot[tt] plots when transitions occur between",
129 "dihedral rotamers of multiplicity 3 and [TT]-oh[tt]",
130 "records a histogram of the times between such transitions,",
131 "assuming the input trajectory frames are equally spaced in time."
133 static const char* opt[] = { nullptr, "angle", "dihedral", "improper", "ryckaert-bellemans",
135 static gmx_bool bALL = FALSE, bChandler = FALSE, bAverCorr = FALSE, bPBC = TRUE;
136 static real binwidth = 1;
138 { "-type", FALSE, etENUM, { opt }, "Type of angle to analyse" },
143 "Plot all angles separately in the averages file, in the order of appearance in the "
149 "binwidth (degrees) for calculating the distribution" },
150 { "-periodic", FALSE, etBOOL, { &bPBC }, "Print dihedral angles modulo 360 degrees" },
155 "Use Chandler correlation function (N[trans] = 1, N[gauche] = 0) rather than cosine "
156 "correlation function. Trans is defined as phi < -60 or phi > 60." },
161 "Average the correlation functions for the individual angles/dihedrals" }
163 static const char* bugs[] = {
164 "Counting transitions only works for dihedrals with multiplicity 3"
172 real maxang, S2, norm_fac, maxstat;
174 int nframes, maxangstat, mult, *angstat;
175 int i, j, nangles, first, last;
176 gmx_bool bAver, bRb, bPeriodic, bFrac, /* calculate fraction too? */
177 bTrans, /* worry about transtions too? */
178 bCorr; /* correlation function ? */
181 real** dih = nullptr; /* mega array with all dih. angles at all times*/
182 real * time, *trans_frac, *aver_angle;
183 t_filenm fnm[] = { { efTRX, "-f", nullptr, ffREAD }, { efNDX, nullptr, "angle", ffREAD },
184 { efXVG, "-od", "angdist", ffWRITE }, { efXVG, "-ov", "angaver", ffOPTWR },
185 { efXVG, "-of", "dihfrac", ffOPTWR }, { efXVG, "-ot", "dihtrans", ffOPTWR },
186 { efXVG, "-oh", "trhisto", ffOPTWR }, { efXVG, "-oc", "dihcorr", ffOPTWR },
187 { efTRR, "-or", nullptr, ffOPTWR } };
188 #define NFILE asize(fnm)
191 gmx_output_env_t* oenv;
194 ppa = add_acf_pargs(&npargs, pa);
195 if (!parse_common_args(
196 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
206 GMX_RELEASE_ASSERT(opt[0] != nullptr,
207 "Internal option inconsistency; opt[0]==NULL after processing");
215 case 'd': // Intended fall through
217 case 'r': bRb = TRUE; break;
220 if (opt2bSet("-or", NFILE, fnm))
224 gmx_fatal(FARGS, "Can not combine angles with trr dump");
228 please_cite(stdout, "Mu2005a");
232 /* Calculate bin size */
233 maxangstat = gmx::roundToInt(maxang / binwidth);
234 binwidth = maxang / maxangstat;
236 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
237 nangles = isize / mult;
238 if ((isize % mult) != 0)
241 "number of index elements not multiple of %d, "
242 "these can not be %s\n",
244 (mult == 3) ? "angle triplets" : "dihedral quadruplets");
248 /* Check whether specific analysis has to be performed */
249 bCorr = opt2bSet("-oc", NFILE, fnm);
250 bAver = opt2bSet("-ov", NFILE, fnm);
251 bTrans = opt2bSet("-ot", NFILE, fnm);
252 bFrac = opt2bSet("-of", NFILE, fnm);
253 if (bTrans && opt[0][0] != 'd')
255 fprintf(stderr, "Option -ot should only accompany -type dihedral. Disabling -ot.\n");
259 if (bChandler && !bCorr)
268 " calculating fractions as defined in this program\n"
269 "makes sense for Ryckaert Bellemans dihs. only. Ignoring -of\n\n");
273 if ((bTrans || bFrac || bCorr) && mult == 3)
276 "Can only do transition, fraction or correlation\n"
277 "on dihedrals. Select -d\n");
281 * We need to know the nr of frames so we can allocate memory for an array
282 * with all dihedral angles at all timesteps. Works for me.
284 if (bTrans || bCorr || bALL || opt2bSet("-or", NFILE, fnm))
289 snew(angstat, maxangstat);
291 read_ang_dih(ftp2fn(efTRX, NFILE, fnm),
293 bALL || bCorr || bTrans || opt2bSet("-or", NFILE, fnm),
307 dt = (time[nframes - 1] - time[0]) / (nframes - 1);
311 sprintf(title, "Average Angle: %s", grpname);
312 out = xvgropen(opt2fn("-ov", NFILE, fnm), title, "Time (ps)", "Angle (degrees)", oenv);
313 for (i = 0; (i < nframes); i++)
315 fprintf(out, "%10.5f %8.3f", time[i], aver_angle[i] * RAD2DEG);
318 for (j = 0; (j < nangles); j++)
323 fprintf(out, " %8.3f", std::atan2(std::sin(dd), std::cos(dd)) * RAD2DEG);
327 fprintf(out, " %8.3f", dih[j][i] * RAD2DEG);
335 if (opt2bSet("-or", NFILE, fnm))
337 dump_dih_trr(nframes, nangles, dih, opt2fn("-or", NFILE, fnm), time);
342 sprintf(title, "Trans fraction: %s", grpname);
343 out = xvgropen(opt2fn("-of", NFILE, fnm), title, "Time (ps)", "Fraction", oenv);
345 for (i = 0; (i < nframes); i++)
347 fprintf(out, "%10.5f %10.3f\n", time[i], trans_frac[i]);
348 tfrac += trans_frac[i];
353 fprintf(stderr, "Average trans fraction: %g\n", tfrac);
360 opt2fn("-ot", NFILE, fnm), opt2fn("-oh", NFILE, fnm), dih, nframes, nangles, grpname, time, bRb, oenv);
365 /* Autocorrelation function */
368 fprintf(stderr, "Not enough frames for correlation function\n");
375 real dval, sixty = DEG2RAD * 60;
378 for (i = 0; (i < nangles); i++)
380 for (j = 0; (j < nframes); j++)
385 bTest = (dval > -sixty) && (dval < sixty);
389 bTest = (dval < -sixty) || (dval > sixty);
393 dih[i][j] = dval - tfrac;
410 do_autocorr(opt2fn("-oc", NFILE, fnm),
412 "Dihedral Autocorrelation Function",
423 /* Determine the non-zero part of the distribution */
424 for (first = 0; (first < maxangstat - 1) && (angstat[first + 1] == 0); first++) {}
425 for (last = maxangstat - 1; (last > 0) && (angstat[last - 1] == 0); last--) {}
428 printf("Found points in the range from %d to %d (max %d)\n", first, last, maxangstat);
429 if (bTrans || bCorr || bALL || opt2bSet("-or", NFILE, fnm))
430 { /* It's better to re-calculate Std. Dev per sample */
431 real b_aver = aver_angle[0];
434 for (int i = 0; (i < nframes); i++)
436 delta = correctRadianAngleRange(aver_angle[i] - b_aver);
439 for (int j = 0; (j < nangles); j++)
441 delta = correctRadianAngleRange(dih[j][i] - b);
447 { /* Incorrect for Std. Dev. */
448 real delta, b_aver = aver_angle[0];
449 for (i = 0; (i < nframes); i++)
451 delta = correctRadianAngleRange(aver_angle[i] - b_aver);
457 double aversig = correctRadianAngleRange(aver);
460 printf(" < angle > = %g\n", aversig);
464 sprintf(title, "Angle Distribution: %s", grpname);
468 sprintf(title, "Dihedral Distribution: %s", grpname);
470 calc_distribution_props(maxangstat, angstat, -180.0, 0, nullptr, &S2);
471 fprintf(stderr, "Order parameter S^2 = %g\n", S2);
474 bPeriodic = (mult == 4) && (first == 0) && (last == maxangstat - 1);
476 out = xvgropen(opt2fn("-od", NFILE, fnm), title, "Degrees", "", oenv);
477 if (output_env_get_print_xvgr_codes(oenv))
479 fprintf(out, "@ subtitle \"average angle: %g\\So\\N\"\n", aver);
481 norm_fac = 1.0 / (nangles * nframes * binwidth);
485 for (i = first; (i <= last); i++)
487 maxstat = std::max(maxstat, angstat[i] * norm_fac);
489 if (output_env_get_print_xvgr_codes(oenv))
491 fprintf(out, "@with g0\n");
492 fprintf(out, "@ world xmin -180\n");
493 fprintf(out, "@ world xmax 180\n");
494 fprintf(out, "@ world ymin 0\n");
495 fprintf(out, "@ world ymax %g\n", maxstat * 1.05);
496 fprintf(out, "@ xaxis tick major 60\n");
497 fprintf(out, "@ xaxis tick minor 30\n");
498 fprintf(out, "@ yaxis tick major 0.005\n");
499 fprintf(out, "@ yaxis tick minor 0.0025\n");
502 for (i = first; (i <= last); i++)
504 fprintf(out, "%10g %10f\n", i * binwidth + 180.0 - maxang, angstat[i] * norm_fac);
508 /* print first bin again as last one */
509 fprintf(out, "%10g %10f\n", 180.0, angstat[0] * norm_fac);
514 do_view(oenv, opt2fn("-od", NFILE, fnm), "-nxy");
517 do_view(oenv, opt2fn("-ov", NFILE, fnm), "-nxy");