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, 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.
44 #include "gromacs/commandline/pargs.h"
45 #include "gromacs/commandline/viewit.h"
46 #include "gromacs/correlationfunctions/autocorr.h"
47 #include "gromacs/fileio/trrio.h"
48 #include "gromacs/fileio/xvgr.h"
49 #include "gromacs/gmxana/gmx_ana.h"
50 #include "gromacs/gmxana/gstat.h"
51 #include "gromacs/math/functions.h"
52 #include "gromacs/math/units.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/topology/index.h"
55 #include "gromacs/utility/arraysize.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/futil.h"
58 #include "gromacs/utility/gmxassert.h"
59 #include "gromacs/utility/pleasecite.h"
60 #include "gromacs/utility/smalloc.h"
62 static void dump_dih_trr(int nframes
, int nangles
, real
**dih
, const char *fn
,
65 int i
, j
, k
, l
, m
, na
;
68 matrix box
= {{2, 0, 0}, {0, 2, 0}, {0, 0, 2}};
79 printf("There are %d dihedrals. Will fill %d atom positions with cos/sin\n",
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
!= NULL
&& dih
[j
] != NULL
, "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
, NULL
, NULL
);
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
[] = { NULL
, "angle", "dihedral", "improper", "ryckaert-bellemans", NULL
};
134 static gmx_bool bALL
= FALSE
, bChandler
= FALSE
, bAverCorr
= FALSE
, bPBC
= TRUE
;
135 static real binwidth
= 1;
137 { "-type", FALSE
, etENUM
, {opt
},
138 "Type of angle to analyse" },
139 { "-all", FALSE
, etBOOL
, {&bALL
},
140 "Plot all angles separately in the averages file, in the order of appearance in the index file." },
141 { "-binwidth", FALSE
, etREAL
, {&binwidth
},
142 "binwidth (degrees) for calculating the distribution" },
143 { "-periodic", FALSE
, etBOOL
, {&bPBC
},
144 "Print dihedral angles modulo 360 degrees" },
145 { "-chandler", FALSE
, etBOOL
, {&bChandler
},
146 "Use Chandler correlation function (N[trans] = 1, N[gauche] = 0) rather than cosine correlation function. Trans is defined as phi < -60 or phi > 60." },
147 { "-avercorr", FALSE
, etBOOL
, {&bAverCorr
},
148 "Average the correlation functions for the individual angles/dihedrals" }
150 static const char *bugs
[] = {
151 "Counting transitions only works for dihedrals with multiplicity 3"
159 real maxang
, S2
, norm_fac
, maxstat
;
161 int nframes
, maxangstat
, mult
, *angstat
;
162 int i
, j
, nangles
, first
, last
;
163 gmx_bool bAver
, bRb
, bPeriodic
,
164 bFrac
, /* calculate fraction too? */
165 bTrans
, /* worry about transtions too? */
166 bCorr
; /* correlation function ? */
167 real aver
, aver2
, aversig
; /* fraction trans dihedrals */
170 real
**dih
= NULL
; /* mega array with all dih. angles at all times*/
171 real
*time
, *trans_frac
, *aver_angle
;
173 { efTRX
, "-f", NULL
, ffREAD
},
174 { efNDX
, NULL
, "angle", ffREAD
},
175 { efXVG
, "-od", "angdist", ffWRITE
},
176 { efXVG
, "-ov", "angaver", ffOPTWR
},
177 { efXVG
, "-of", "dihfrac", ffOPTWR
},
178 { efXVG
, "-ot", "dihtrans", ffOPTWR
},
179 { efXVG
, "-oh", "trhisto", ffOPTWR
},
180 { efXVG
, "-oc", "dihcorr", ffOPTWR
},
181 { efTRR
, "-or", NULL
, ffOPTWR
}
183 #define NFILE asize(fnm)
186 gmx_output_env_t
*oenv
;
189 ppa
= add_acf_pargs(&npargs
, pa
);
190 if (!parse_common_args(&argc
, argv
, PCA_CAN_VIEW
| PCA_CAN_TIME
,
191 NFILE
, fnm
, npargs
, ppa
, asize(desc
), desc
, asize(bugs
), bugs
,
201 GMX_RELEASE_ASSERT(opt
[0] != NULL
, "Internal option inconsistency; opt[0]==NULL after processing");
218 if (opt2bSet("-or", NFILE
, fnm
))
222 gmx_fatal(FARGS
, "Can not combine angles with trr dump");
226 please_cite(stdout
, "Mu2005a");
230 /* Calculate bin size */
231 maxangstat
= static_cast<int>(maxang
/binwidth
+0.5);
232 binwidth
= maxang
/maxangstat
;
234 rd_index(ftp2fn(efNDX
, NFILE
, fnm
), 1, &isize
, &index
, &grpname
);
235 nangles
= isize
/mult
;
236 if ((isize
% mult
) != 0)
238 gmx_fatal(FARGS
, "number of index elements not multiple of %d, "
239 "these can not be %s\n",
240 mult
, (mult
== 3) ? "angle triplets" : "dihedral quadruplets");
244 /* Check whether specific analysis has to be performed */
245 bCorr
= opt2bSet("-oc", NFILE
, fnm
);
246 bAver
= opt2bSet("-ov", NFILE
, fnm
);
247 bTrans
= opt2bSet("-ot", NFILE
, fnm
);
248 bFrac
= opt2bSet("-of", NFILE
, fnm
);
249 if (bTrans
&& opt
[0][0] != 'd')
251 fprintf(stderr
, "Option -ot should only accompany -type dihedral. Disabling -ot.\n");
255 if (bChandler
&& !bCorr
)
262 fprintf(stderr
, "Warning:"
263 " calculating fractions as defined in this program\n"
264 "makes sense for Ryckaert Bellemans dihs. only. Ignoring -of\n\n");
268 if ( (bTrans
|| bFrac
|| bCorr
) && mult
== 3)
270 gmx_fatal(FARGS
, "Can only do transition, fraction or correlation\n"
271 "on dihedrals. Select -d\n");
275 * We need to know the nr of frames so we can allocate memory for an array
276 * with all dihedral angles at all timesteps. Works for me.
278 if (bTrans
|| bCorr
|| bALL
|| opt2bSet("-or", NFILE
, fnm
))
283 snew(angstat
, maxangstat
);
285 read_ang_dih(ftp2fn(efTRX
, NFILE
, fnm
), (mult
== 3),
286 bALL
|| bCorr
|| bTrans
|| opt2bSet("-or", NFILE
, fnm
),
287 bRb
, bPBC
, maxangstat
, angstat
,
288 &nframes
, &time
, isize
, index
, &trans_frac
, &aver_angle
, dih
,
291 dt
= (time
[nframes
-1]-time
[0])/(nframes
-1);
295 sprintf(title
, "Average Angle: %s", grpname
);
296 out
= xvgropen(opt2fn("-ov", NFILE
, fnm
),
297 title
, "Time (ps)", "Angle (degrees)", oenv
);
298 for (i
= 0; (i
< nframes
); i
++)
300 fprintf(out
, "%10.5f %8.3f", time
[i
], aver_angle
[i
]*RAD2DEG
);
303 for (j
= 0; (j
< nangles
); j
++)
308 fprintf(out
, " %8.3f", std::atan2(std::sin(dd
), std::cos(dd
))*RAD2DEG
);
312 fprintf(out
, " %8.3f", dih
[j
][i
]*RAD2DEG
);
320 if (opt2bSet("-or", NFILE
, fnm
))
322 dump_dih_trr(nframes
, nangles
, dih
, opt2fn("-or", NFILE
, fnm
), time
);
327 sprintf(title
, "Trans fraction: %s", grpname
);
328 out
= xvgropen(opt2fn("-of", NFILE
, fnm
),
329 title
, "Time (ps)", "Fraction", oenv
);
331 for (i
= 0; (i
< nframes
); i
++)
333 fprintf(out
, "%10.5f %10.3f\n", time
[i
], trans_frac
[i
]);
334 tfrac
+= trans_frac
[i
];
339 fprintf(stderr
, "Average trans fraction: %g\n", tfrac
);
345 ana_dih_trans(opt2fn("-ot", NFILE
, fnm
), opt2fn("-oh", NFILE
, fnm
),
346 dih
, nframes
, nangles
, grpname
, time
, bRb
, oenv
);
351 /* Autocorrelation function */
354 fprintf(stderr
, "Not enough frames for correlation function\n");
361 real dval
, sixty
= DEG2RAD
*60;
364 for (i
= 0; (i
< nangles
); i
++)
366 for (j
= 0; (j
< nframes
); j
++)
371 bTest
= (dval
> -sixty
) && (dval
< sixty
);
375 bTest
= (dval
< -sixty
) || (dval
> sixty
);
379 dih
[i
][j
] = dval
-tfrac
;
396 do_autocorr(opt2fn("-oc", NFILE
, fnm
), oenv
,
397 "Dihedral Autocorrelation Function",
398 nframes
, nangles
, dih
, dt
, mode
, bAverCorr
);
403 /* Determine the non-zero part of the distribution */
404 for (first
= 0; (first
< maxangstat
-1) && (angstat
[first
+1] == 0); first
++)
408 for (last
= maxangstat
-1; (last
> 0) && (angstat
[last
-1] == 0); last
--)
414 for (i
= 0; (i
< nframes
); i
++)
416 aver
+= RAD2DEG
*aver_angle
[i
];
417 aver2
+= gmx::square(RAD2DEG
*aver_angle
[i
]);
421 aversig
= std::sqrt(aver2
-gmx::square(aver
));
422 printf("Found points in the range from %d to %d (max %d)\n",
423 first
, last
, maxangstat
);
424 printf(" < angle > = %g\n", aver
);
425 printf("< angle^2 > = %g\n", aver2
);
426 printf("Std. Dev. = %g\n", aversig
);
430 sprintf(title
, "Angle Distribution: %s", grpname
);
434 sprintf(title
, "Dihedral Distribution: %s", grpname
);
436 calc_distribution_props(maxangstat
, angstat
, -180.0, 0, NULL
, &S2
);
437 fprintf(stderr
, "Order parameter S^2 = %g\n", S2
);
440 bPeriodic
= (mult
== 4) && (first
== 0) && (last
== maxangstat
-1);
442 out
= xvgropen(opt2fn("-od", NFILE
, fnm
), title
, "Degrees", "", oenv
);
443 if (output_env_get_print_xvgr_codes(oenv
))
445 fprintf(out
, "@ subtitle \"average angle: %g\\So\\N\"\n", aver
);
447 norm_fac
= 1.0/(nangles
*nframes
*binwidth
);
451 for (i
= first
; (i
<= last
); i
++)
453 maxstat
= std::max(maxstat
, angstat
[i
]*norm_fac
);
455 if (output_env_get_print_xvgr_codes(oenv
))
457 fprintf(out
, "@with g0\n");
458 fprintf(out
, "@ world xmin -180\n");
459 fprintf(out
, "@ world xmax 180\n");
460 fprintf(out
, "@ world ymin 0\n");
461 fprintf(out
, "@ world ymax %g\n", maxstat
*1.05);
462 fprintf(out
, "@ xaxis tick major 60\n");
463 fprintf(out
, "@ xaxis tick minor 30\n");
464 fprintf(out
, "@ yaxis tick major 0.005\n");
465 fprintf(out
, "@ yaxis tick minor 0.0025\n");
468 for (i
= first
; (i
<= last
); i
++)
470 fprintf(out
, "%10g %10f\n", i
*binwidth
+180.0-maxang
, angstat
[i
]*norm_fac
);
474 /* print first bin again as last one */
475 fprintf(out
, "%10g %10f\n", 180.0, angstat
[0]*norm_fac
);
480 do_view(oenv
, opt2fn("-od", NFILE
, fnm
), "-nxy");
483 do_view(oenv
, opt2fn("-ov", NFILE
, fnm
), "-nxy");