Gromacs
2016.6
|
#include "gmxpre.h"
#include "config.h"
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <sstream>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/xvgr.h"
#include "gromacs/gmxana/gmx_ana.h"
#include "gromacs/math/functions.h"
#include "gromacs/math/units.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/pull-params.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/random/tabulatednormaldistribution.h"
#include "gromacs/random/threefry.h"
#include "gromacs/random/uniformintdistribution.h"
#include "gromacs/random/uniformrealdistribution.h"
#include "gromacs/utility/arraysize.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxomp.h"
#include "gromacs/utility/pleasecite.h"
#include "gromacs/utility/smalloc.h"
Implementation of the Weighted Histogram Analysis Method (WHAM)
Classes | |
struct | t_pullcoord |
Parameters of one pull coodinate. More... | |
struct | t_UmbrellaHeader |
Parameters of the umbrella potentials. More... | |
struct | t_UmbrellaWindow |
Data in the umbrella histograms. More... | |
struct | t_coordselection |
Selection of pull coordinates to be used in WHAM (one structure for each tpr file) More... | |
struct | t_UmbrellaOptions |
Parameters of WHAM. More... | |
Macros | |
#define | WHAM_MAXFILELEN 2048 |
longest file names allowed in input files | |
#define | int2YN(a) (((a) == 0) ? ("N") : ("Y")) |
translate 0/1 to N/Y to write pull dimensions | |
#define | WHAM_AC_ZERO_LIMIT 0.05 |
Stop integrating autoccorelation function when ACF drops under this value. | |
#define | WHAMBOOLXOR(a, b) ( ((!(a)) && (b)) || ((a) && (!(b)))) |
Boolean XOR. | |
#define | NFILE asize(fnm) |
Number of elements in fnm (used for command line parsing) | |
Enumerations | |
enum | { enSel, en_kJ, en_kCal, en_kT, enNr } |
enum for energy units | |
enum | { whamin_unknown, whamin_tpr, whamin_pullxf, whamin_pdo } |
enum for type of input files (pdos, tpr, or pullf) | |
enum | { bsMethod_unknown, bsMethod_BayesianHist, bsMethod_hist, bsMethod_traj, bsMethod_trajGauss } |
enum for bootstrapping method More... | |
Functions | |
t_UmbrellaWindow * | initUmbrellaWindows (int nwin) |
Make an umbrella window (may contain several histograms) | |
void | freeUmbrellaWindows (t_UmbrellaWindow *win, int nwin) |
Delete an umbrella window (may contain several histograms) | |
void | setup_tab (const char *fn, t_UmbrellaOptions *opt) |
Read and setup tabulated umbrella potential. | |
void | read_pdo_header (FILE *file, t_UmbrellaHeader *header, t_UmbrellaOptions *opt) |
Read the header of an PDO file (position, force const, nr of groups) | |
static char * | fgets3 (FILE *fp, char ptr[], int *len) |
smarter fgets | |
void | read_pdo_data (FILE *file, t_UmbrellaHeader *header, int fileno, t_UmbrellaWindow *win, t_UmbrellaOptions *opt, gmx_bool bGetMinMax, real *mintmp, real *maxtmp) |
Read the data columns of and PDO file. More... | |
void | enforceEqualWeights (t_UmbrellaWindow *window, int nWindows) |
Set identical weights for all histograms. More... | |
double | tabulated_pot (double dist, t_UmbrellaOptions *opt) |
Simple linear interpolation between two given tabulated points. | |
void | setup_acc_wham (double *profile, t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt) |
Check which bins substiantially contribute (accelerates WHAM) More... | |
void | calc_profile (double *profile, t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt, gmx_bool bExact) |
Compute the PMF (one of the two main WHAM routines) | |
double | calc_z (double *profile, t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt, gmx_bool bExact) |
Compute the free energy offsets z (one of the two main WHAM routines) | |
void | symmetrizeProfile (double *profile, t_UmbrellaOptions *opt) |
Make PMF symmetric around 0 (useful e.g. for membranes) | |
void | prof_normalization_and_unit (double *profile, t_UmbrellaOptions *opt) |
Set energy unit (kJ/mol,kT,kCal/mol) and set it to zero at opt->zProf0. | |
void | getRandomIntArray (int nPull, int blockLength, int *randomArray, gmx::DefaultRandomEngine *rng) |
Make an array of random integers (used for bootstrapping) | |
void | copy_pullgrp_to_synthwindow (t_UmbrellaWindow *synthWindow, t_UmbrellaWindow *thisWindow, int pullid) |
Set pull group information of a synthetic histogram. More... | |
void | calc_cumulatives (t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt, const char *fnhist, const char *xlabel) |
Calculate cumulative distribution function of of all histograms. More... | |
void | searchCumulative (double xx[], int n, double x, int *j) |
Return j such that xx[j] <= x < xx[j+1]. More... | |
void | create_synthetic_histo (t_UmbrellaWindow *synthWindow, t_UmbrellaWindow *thisWindow, int pullid, t_UmbrellaOptions *opt) |
Bootstrap new trajectories and thereby generate new (bootstrapped) histograms. | |
void | print_histograms (const char *fnhist, t_UmbrellaWindow *window, int nWindows, int bs_index, t_UmbrellaOptions *opt, const char *xlabel) |
Write all histograms to a file. More... | |
int | func_wham_is_larger (const void *a, const void *b) |
Used for qsort to sort random numbers. | |
void | setRandomBsWeights (t_UmbrellaWindow *synthwin, int nAllPull, t_UmbrellaOptions *opt) |
Make random weights for histograms for the Bayesian bootstrap of complete histograms) | |
void | do_bootstrapping (const char *fnres, const char *fnprof, const char *fnhist, const char *xlabel, char *ylabel, double *profile, t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt) |
The main bootstrapping routine. | |
int | whaminFileType (char *fn) |
Return type of input file based on file extension (xvg, pdo, or tpr) | |
void | read_wham_in (const char *fn, char ***filenamesRet, int *nfilesRet, t_UmbrellaOptions *opt) |
Read the files names in pdo-files.dat, pullf/x-files.dat, tpr-files.dat. | |
FILE * | open_pdo_pipe (const char *fn, t_UmbrellaOptions *opt, gmx_bool *bPipeOpen) |
Open a file or a pipe to a gzipped file. | |
void | pdo_close_file (FILE *fp) |
Close file or pipe. | |
void | read_pdo_files (char **fn, int nfiles, t_UmbrellaHeader *header, t_UmbrellaWindow *window, t_UmbrellaOptions *opt) |
Reading all pdo files. | |
void | read_tpr_header (const char *fn, t_UmbrellaHeader *header, t_UmbrellaOptions *opt, t_coordselection *coordsel) |
Read pull groups from a tpr file (including position, force const, geometry, number of groups) | |
double | dist_ndim (double **dx, int ndim, int line) |
2-norm in a ndim-dimensional space | |
void | read_pull_xf (const char *fn, t_UmbrellaHeader *header, t_UmbrellaWindow *window, t_UmbrellaOptions *opt, gmx_bool bGetMinMax, real *mintmp, real *maxtmp, t_coordselection *coordsel) |
Read pullx.xvg or pullf.xvg. | |
void | read_tpr_pullxf_files (char **fnTprs, char **fnPull, int nfiles, t_UmbrellaHeader *header, t_UmbrellaWindow *window, t_UmbrellaOptions *opt) |
read pullf-files.dat or pullx-files.dat and tpr-files.dat | |
void | readIntegratedAutocorrelationTimes (t_UmbrellaWindow *window, int nwins, const char *fn) |
Read integrated autocorrelation times from input file (option -iiact) More... | |
void | smoothIact (t_UmbrellaWindow *window, int nwins, t_UmbrellaOptions *opt) |
Smooth autocorreltion times along the reaction coordinate. More... | |
void | calcIntegratedAutocorrelationTimes (t_UmbrellaWindow *window, int nwins, t_UmbrellaOptions *opt, const char *fn, const char *xlabel) |
Try to compute the autocorrelation time for each umbrealla window. | |
void | averageSigma (t_UmbrellaWindow *window, int nwins) |
compute average and sigma of each umbrella histogram | |
void | computeAverageForce (t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt) |
Use histograms to compute average force on pull group. | |
void | checkReactionCoordinateCovered (t_UmbrellaWindow *window, int nwins, t_UmbrellaOptions *opt) |
Check if the complete reaction coordinate is covered by the histograms. | |
void | guessPotByIntegration (t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt, const char *xlabel) |
Compute initial potential by integrating the average force. More... | |
static int | wordcount (char *ptr) |
Count number of words in an line. | |
void | readPullCoordSelection (t_UmbrellaOptions *opt, char **fnTpr, int nTpr) |
Read input file for pull group selection (option -is) More... | |
int | gmx_wham (int argc, char *argv[]) |
The main gmx wham routine. | |
anonymous enum |
enum for bootstrapping method
These bootstrap methods are supported:
FOR MORE DETAILS ON THE BOOTSTRAP METHODS (INCLUDING EXAMPLES), SEE JS Hub, BL de Groot, D van der Spoel g_wham - A free weighted histogram analysis implementation including robust error and autocorrelation estimates, J Chem Theory Comput, 6(12), 3713-3720 (2010)
void calc_cumulatives | ( | t_UmbrellaWindow * | window, |
int | nWindows, | ||
t_UmbrellaOptions * | opt, | ||
const char * | fnhist, | ||
const char * | xlabel | ||
) |
Calculate cumulative distribution function of of all histograms.
This allow to create random number sequences which are distributed according to the histograms. Required to generate the "synthetic" histograms for the Bootstrap method
void copy_pullgrp_to_synthwindow | ( | t_UmbrellaWindow * | synthWindow, |
t_UmbrellaWindow * | thisWindow, | ||
int | pullid | ||
) |
Set pull group information of a synthetic histogram.
This is used when bootstapping new trajectories and thereby create new histogtrams, but it is not required if we bootstrap complete histograms.
void enforceEqualWeights | ( | t_UmbrellaWindow * | window, |
int | nWindows | ||
) |
Set identical weights for all histograms.
Normally, the weight is given by the number data points in each histogram, together with the autocorrelation time. This can be overwritten by this routine (not recommended). Since we now support autocorrelations, it is better to set an appropriate autocorrelation times instead of using this function.
void guessPotByIntegration | ( | t_UmbrellaWindow * | window, |
int | nWindows, | ||
t_UmbrellaOptions * | opt, | ||
const char * | xlabel | ||
) |
Compute initial potential by integrating the average force.
This speeds up the convergence by roughly a factor of 2
void print_histograms | ( | const char * | fnhist, |
t_UmbrellaWindow * | window, | ||
int | nWindows, | ||
int | bs_index, | ||
t_UmbrellaOptions * | opt, | ||
const char * | xlabel | ||
) |
Write all histograms to a file.
If bs_index>=0, a number is added to the output file name to allow the ouput of all sets of bootstrapped histograms.
void read_pdo_data | ( | FILE * | file, |
t_UmbrellaHeader * | header, | ||
int | fileno, | ||
t_UmbrellaWindow * | win, | ||
t_UmbrellaOptions * | opt, | ||
gmx_bool | bGetMinMax, | ||
real * | mintmp, | ||
real * | maxtmp | ||
) |
Read the data columns of and PDO file.
TO DO: Get rid of the scanf function to avoid the clang warning. At the moment, this warning is avoided by hiding the format string the variable fmtlf.
void readIntegratedAutocorrelationTimes | ( | t_UmbrellaWindow * | window, |
int | nwins, | ||
const char * | fn | ||
) |
Read integrated autocorrelation times from input file (option -iiact)
Note: Here we consider tau[int] := int_0^inf ACF(t) as the integrated autocorrelation times. The factor g := 1 + 2*tau[int]
subsequently enters the uncertainty.
void readPullCoordSelection | ( | t_UmbrellaOptions * | opt, |
char ** | fnTpr, | ||
int | nTpr | ||
) |
Read input file for pull group selection (option -is)
TO DO: ptr=fgets(...) is never freed (small memory leak)
void searchCumulative | ( | double | xx[], |
int | n, | ||
double | x, | ||
int * | j | ||
) |
Return j such that xx[j] <= x < xx[j+1].
This is used to generate a random sequence distributed according to a histogram
void setup_acc_wham | ( | double * | profile, |
t_UmbrellaWindow * | window, | ||
int | nWindows, | ||
t_UmbrellaOptions * | opt | ||
) |
Check which bins substiantially contribute (accelerates WHAM)
Don't worry, that routine does not mean we compute the PMF in limited precision. After rapid convergence (using only substiantal contributions), we always switch to full precision.
void smoothIact | ( | t_UmbrellaWindow * | window, |
int | nwins, | ||
t_UmbrellaOptions * | opt | ||
) |
Smooth autocorreltion times along the reaction coordinate.
This is useful if the ACT is subject to high uncertainty in case if limited sampling. Note that -in case of limited sampling- the ACT may be severely underestimated. Note: the g=1+2tau are overwritten. If opt->bAllowReduceIact==FALSE, the ACTs are never reduced, only increased by the smoothing