3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
51 #include "gmx_fatal.h"
58 static void dump_dih_trn(int nframes, int nangles, real **dih, const char *fn,
61 int i, j, k, l, m, na;
64 matrix box = {{2, 0, 0}, {0, 2, 0}, {0, 0, 2}};
75 printf("There are %d dihedrals. Will fill %d atom positions with cos/sin\n",
78 trn = open_trn(fn, "w");
79 for (i = 0; (i < nframes); i++)
82 for (j = 0; (j < nangles); j++)
84 for (m = 0; (m < 2); m++)
86 x[k][l] = (m == 0) ? cos(dih[j][i]) : sin(dih[j][i]);
95 fwrite_trn(trn, i, time[i], 0, box, na, x, NULL, NULL);
101 int gmx_g_angle(int argc, char *argv[])
103 static const char *desc[] = {
104 "[TT]g_angle[tt] computes the angle distribution for a number of angles",
105 "or dihedrals.[PAR]",
106 "With option [TT]-ov[tt], you can plot the average angle of",
107 "a group of angles as a function of time. With the [TT]-all[tt] option,",
108 "the first graph is the average and the rest are the individual angles.[PAR]",
109 "With the [TT]-of[tt] option, [TT]g_angle[tt] also calculates the fraction of trans",
110 "dihedrals (only for dihedrals) as function of time, but this is",
111 "probably only fun for a select few.[PAR]",
112 "With option [TT]-oc[tt], a dihedral correlation function is calculated.[PAR]",
113 "It should be noted that the index file must contain",
114 "atom triplets for angles or atom quadruplets for dihedrals.",
115 "If this is not the case, the program will crash.[PAR]",
116 "With option [TT]-or[tt], a trajectory file is dumped containing cos and",
117 "sin of selected dihedral angles, which subsequently can be used as",
118 "input for a principal components analysis using [TT]g_covar[tt].[PAR]",
119 "Option [TT]-ot[tt] plots when transitions occur between",
120 "dihedral rotamers of multiplicity 3 and [TT]-oh[tt]",
121 "records a histogram of the times between such transitions,",
122 "assuming the input trajectory frames are equally spaced in time."
124 static const char *opt[] = { NULL, "angle", "dihedral", "improper", "ryckaert-bellemans", NULL };
125 static gmx_bool bALL = FALSE, bChandler = FALSE, bAverCorr = FALSE, bPBC = TRUE;
126 static real binwidth = 1;
128 { "-type", FALSE, etENUM, {opt},
129 "Type of angle to analyse" },
130 { "-all", FALSE, etBOOL, {&bALL},
131 "Plot all angles separately in the averages file, in the order of appearance in the index file." },
132 { "-binwidth", FALSE, etREAL, {&binwidth},
133 "binwidth (degrees) for calculating the distribution" },
134 { "-periodic", FALSE, etBOOL, {&bPBC},
135 "Print dihedral angles modulo 360 degrees" },
136 { "-chandler", FALSE, etBOOL, {&bChandler},
137 "Use Chandler correlation function (N[trans] = 1, N[gauche] = 0) rather than cosine correlation function. Trans is defined as phi < -60 or phi > 60." },
138 { "-avercorr", FALSE, etBOOL, {&bAverCorr},
139 "Average the correlation functions for the individual angles/dihedrals" }
141 static const char *bugs[] = {
142 "Counting transitions only works for dihedrals with multiplicity 3"
150 real maxang, Jc, S2, norm_fac, maxstat;
152 int nframes, maxangstat, mult, *angstat;
153 int i, j, total, nangles, natoms, nat2, first, last, angind;
154 gmx_bool bAver, bRb, bPeriodic,
155 bFrac, /* calculate fraction too? */
156 bTrans, /* worry about transtions too? */
157 bCorr; /* correlation function ? */
158 real t, aa, aver, aver2, aversig, fraction; /* fraction trans dihedrals */
161 real **dih = NULL; /* mega array with all dih. angles at all times*/
163 real *time, *trans_frac, *aver_angle;
165 { efTRX, "-f", NULL, ffREAD },
166 { efNDX, NULL, "angle", ffREAD },
167 { efXVG, "-od", "angdist", ffWRITE },
168 { efXVG, "-ov", "angaver", ffOPTWR },
169 { efXVG, "-of", "dihfrac", ffOPTWR },
170 { efXVG, "-ot", "dihtrans", ffOPTWR },
171 { efXVG, "-oh", "trhisto", ffOPTWR },
172 { efXVG, "-oc", "dihcorr", ffOPTWR },
173 { efTRR, "-or", NULL, ffOPTWR }
175 #define NFILE asize(fnm)
181 ppa = add_acf_pargs(&npargs, pa);
182 parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
183 NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
204 if (opt2bSet("-or", NFILE, fnm))
208 gmx_fatal(FARGS, "Can not combine angles with trn dump");
212 please_cite(stdout, "Mu2005a");
216 /* Calculate bin size */
217 maxangstat = (int)(maxang/binwidth+0.5);
218 binwidth = maxang/maxangstat;
220 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
221 nangles = isize/mult;
222 if ((isize % mult) != 0)
224 gmx_fatal(FARGS, "number of index elements not multiple of %d, "
225 "these can not be %s\n",
226 mult, (mult == 3) ? "angle triplets" : "dihedral quadruplets");
230 /* Check whether specific analysis has to be performed */
231 bCorr = opt2bSet("-oc", NFILE, fnm);
232 bAver = opt2bSet("-ov", NFILE, fnm);
233 bTrans = opt2bSet("-ot", NFILE, fnm);
234 bFrac = opt2bSet("-of", NFILE, fnm);
235 if (bTrans && opt[0][0] != 'd')
237 fprintf(stderr, "Option -ot should only accompany -type dihedral. Disabling -ot.\n");
241 if (bChandler && !bCorr)
248 fprintf(stderr, "Warning:"
249 " calculating fractions as defined in this program\n"
250 "makes sense for Ryckaert Bellemans dihs. only. Ignoring -of\n\n");
254 if ( (bTrans || bFrac || bCorr) && mult == 3)
256 gmx_fatal(FARGS, "Can only do transition, fraction or correlation\n"
257 "on dihedrals. Select -d\n");
261 * We need to know the nr of frames so we can allocate memory for an array
262 * with all dihedral angles at all timesteps. Works for me.
264 if (bTrans || bCorr || bALL || opt2bSet("-or", NFILE, fnm))
269 snew(angstat, maxangstat);
271 read_ang_dih(ftp2fn(efTRX, NFILE, fnm), (mult == 3),
272 bALL || bCorr || bTrans || opt2bSet("-or", NFILE, fnm),
273 bRb, bPBC, maxangstat, angstat,
274 &nframes, &time, isize, index, &trans_frac, &aver_angle, dih,
277 dt = (time[nframes-1]-time[0])/(nframes-1);
281 sprintf(title, "Average Angle: %s", grpname);
282 out = xvgropen(opt2fn("-ov", NFILE, fnm),
283 title, "Time (ps)", "Angle (degrees)", oenv);
284 for (i = 0; (i < nframes); i++)
286 fprintf(out, "%10.5f %8.3f", time[i], aver_angle[i]*RAD2DEG);
289 for (j = 0; (j < nangles); j++)
294 fprintf(out, " %8.3f", atan2(sin(dd), cos(dd))*RAD2DEG);
298 fprintf(out, " %8.3f", dih[j][i]*RAD2DEG);
306 if (opt2bSet("-or", NFILE, fnm))
308 dump_dih_trn(nframes, nangles, dih, opt2fn("-or", NFILE, fnm), time);
313 sprintf(title, "Trans fraction: %s", grpname);
314 out = xvgropen(opt2fn("-of", NFILE, fnm),
315 title, "Time (ps)", "Fraction", oenv);
317 for (i = 0; (i < nframes); i++)
319 fprintf(out, "%10.5f %10.3f\n", time[i], trans_frac[i]);
320 tfrac += trans_frac[i];
325 fprintf(stderr, "Average trans fraction: %g\n", tfrac);
331 ana_dih_trans(opt2fn("-ot", NFILE, fnm), opt2fn("-oh", NFILE, fnm),
332 dih, nframes, nangles, grpname, time, bRb, oenv);
337 /* Autocorrelation function */
340 fprintf(stderr, "Not enough frames for correlation function\n");
347 real dval, sixty = DEG2RAD*60;
350 for (i = 0; (i < nangles); i++)
352 for (j = 0; (j < nframes); j++)
357 bTest = (dval > -sixty) && (dval < sixty);
361 bTest = (dval < -sixty) || (dval > sixty);
365 dih[i][j] = dval-tfrac;
382 do_autocorr(opt2fn("-oc", NFILE, fnm), oenv,
383 "Dihedral Autocorrelation Function",
384 nframes, nangles, dih, dt, mode, bAverCorr);
389 /* Determine the non-zero part of the distribution */
390 for (first = 0; (first < maxangstat-1) && (angstat[first+1] == 0); first++)
394 for (last = maxangstat-1; (last > 0) && (angstat[last-1] == 0); last--)
400 for (i = 0; (i < nframes); i++)
402 aver += RAD2DEG*aver_angle[i];
403 aver2 += sqr(RAD2DEG*aver_angle[i]);
405 aver /= (real) nframes;
406 aver2 /= (real) nframes;
407 aversig = sqrt(aver2-sqr(aver));
408 printf("Found points in the range from %d to %d (max %d)\n",
409 first, last, maxangstat);
410 printf(" < angle > = %g\n", aver);
411 printf("< angle^2 > = %g\n", aver2);
412 printf("Std. Dev. = %g\n", aversig);
416 sprintf(title, "Angle Distribution: %s", grpname);
420 sprintf(title, "Dihedral Distribution: %s", grpname);
422 calc_distribution_props(maxangstat, angstat, -180.0, 0, NULL, &S2);
423 fprintf(stderr, "Order parameter S^2 = %g\n", S2);
426 bPeriodic = (mult == 4) && (first == 0) && (last == maxangstat-1);
428 out = xvgropen(opt2fn("-od", NFILE, fnm), title, "Degrees", "", oenv);
429 if (output_env_get_print_xvgr_codes(oenv))
431 fprintf(out, "@ subtitle \"average angle: %g\\So\\N\"\n", aver);
433 norm_fac = 1.0/(nangles*nframes*binwidth);
437 for (i = first; (i <= last); i++)
439 maxstat = max(maxstat, angstat[i]*norm_fac);
441 fprintf(out, "@with g0\n");
442 fprintf(out, "@ world xmin -180\n");
443 fprintf(out, "@ world xmax 180\n");
444 fprintf(out, "@ world ymin 0\n");
445 fprintf(out, "@ world ymax %g\n", maxstat*1.05);
446 fprintf(out, "@ xaxis tick major 60\n");
447 fprintf(out, "@ xaxis tick minor 30\n");
448 fprintf(out, "@ yaxis tick major 0.005\n");
449 fprintf(out, "@ yaxis tick minor 0.0025\n");
451 for (i = first; (i <= last); i++)
453 fprintf(out, "%10g %10f\n", i*binwidth+180.0-maxang, angstat[i]*norm_fac);
457 /* print first bin again as last one */
458 fprintf(out, "%10g %10f\n", 180.0, angstat[0]*norm_fac);
463 do_view(oenv, opt2fn("-od", NFILE, fnm), "-nxy");
466 do_view(oenv, opt2fn("-ov", NFILE, fnm), "-nxy");