#include "gmxpre.h"
#include "config.h"
#include <cassert>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <filesystem>
#include <memory>
#include <string>
#include <vector>
#include "gromacs/commandline/filenm.h"
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/filetypes.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/math/vectypes.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/pull_params.h"
#include "gromacs/mdtypes/state.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/random/seed.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/basedefinitions.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/path.h"
#include "gromacs/utility/pleasecite.h"
#include "gromacs/utility/real.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/stringutil.h"

Include dependency graph for gmx_wham.cpp:

Description

Implementation of the Weighted Histogram Analysis Method (WHAM)

Author: Jochen Hub jhub@.nosp@m.gwdg.nosp@m..de

Classes
struct	t_pullcoord
	Parameters of one pull coordinate. 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	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)

Typedefs
typedef struct UmbrellaOptions	t_UmbrellaOptions
	Parameters of WHAM.

Enumerations
enum	{ enSel, en_kJ, en_kCal, en_kT, enNr }
	enum for energy units

enum	{ whamin_unknown, whamin_tpr, whamin_pullxf }
	enum for type of input files (tpr or pullx/pullf)

enum	{ bsMethod_unknown, bsMethod_BayesianHist, bsMethod_hist, bsMethod_traj, bsMethod_trajGauss }
	enum for bootstrapping method More...

Functions
static t_UmbrellaWindow *	initUmbrellaWindows (int nwin)
	Make an umbrella window (may contain several histograms)

static void	freeUmbrellaWindows (t_UmbrellaWindow *win, int nwin)
	Delete an umbrella window (may contain several histograms)

static void	setup_tab (const char fn, t_UmbrellaOptions opt)
	Read and setup tabulated umbrella potential.

static char *	fgets3 (FILE fp, char ptr[], int len)
	smarter fgets

static void	enforceEqualWeights (t_UmbrellaWindow *window, int nWindows)
	Set identical weights for all histograms. More...

static double	tabulated_pot (double dist, t_UmbrellaOptions *opt)
	Simple linear interpolation between two given tabulated points.

static void	setup_acc_wham (const double profile, t_UmbrellaWindow window, int nWindows, t_UmbrellaOptions *opt)
	Check which bins substiantially contribute (accelerates WHAM) More...

static 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)

static double	calc_z (const 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)

static void	symmetrizeProfile (double profile, t_UmbrellaOptions opt)
	Make PMF symmetric around 0 (useful e.g. for membranes)

static 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.

static void	getRandomIntArray (int nPull, int blockLength, int randomArray, gmx::DefaultRandomEngine rng)
	Make an array of random integers (used for bootstrapping)

static void	copy_pullgrp_to_synthwindow (t_UmbrellaWindow synthWindow, t_UmbrellaWindow thisWindow, int pullid)
	Set pull group information of a synthetic histogram. More...

static 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...

static void	searchCumulative (const double xx[], int n, double x, int *j)
	Return j such that xx[j] <= x < xx[j+1]. More...

static void	create_synthetic_histo (t_UmbrellaWindow synthWindow, t_UmbrellaWindow thisWindow, int pullid, t_UmbrellaOptions *opt)
	Bootstrap new trajectories and thereby generate new (bootstrapped) histograms.

static 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...

static void	setRandomBsWeights (t_UmbrellaWindow synthwin, int nAllPull, t_UmbrellaOptions opt)
	Make random weights for histograms for the Bayesian bootstrap of complete histograms)

static 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.

static int	whaminFileType (char *fn)
	Return type of input file based on file extension (xvg or tpr)

static void	read_wham_in (const char fn, char *filenamesRet, int nfilesRet, t_UmbrellaOptions *opt)
	Read the files names in pullf/pullx-files.dat, tpr-files.dat.

static 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)

static 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.

static 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

static void	readIntegratedAutocorrelationTimes (t_UmbrellaWindow window, int nwins, const char fn)
	Read integrated autocorrelation times from input file (option -iiact) More...

static void	smoothIact (t_UmbrellaWindow window, int nwins, t_UmbrellaOptions opt)
	Smooth autocorreltion times along the reaction coordinate. More...

static 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.

static void	averageSigma (t_UmbrellaWindow *window, int nwins)
	compute average and sigma of each umbrella histogram

static void	computeAverageForce (t_UmbrellaWindow window, int nWindows, t_UmbrellaOptions opt)
	Use histograms to compute average force on pull group.

static void	checkReactionCoordinateCovered (t_UmbrellaWindow window, int nwins, t_UmbrellaOptions opt)
	Check if the complete reaction coordinate is covered by the histograms.

static 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.

static 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.

Enumeration Type Documentation

anonymous enum

enum for bootstrapping method

These bootstrap methods are supported:

bootstrap complete histograms with continuous weights (Bayesian bootstrap) (bsMethod_BayesianHist)
bootstrap complete histograms (bsMethod_hist)
bootstrap trajectories from given umbrella histograms. This generates new "synthetic" histograms (bsMethod_traj)
bootstrap trajectories from Gaussian with mu/sigma computed from the respective histogram (bsMethod_trajGauss). This gives very similar results compared to bsMethod_traj.

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)

Function Documentation

static void calc_cumulatives	(	t_UmbrellaWindow *	window,
		int	nWindows,
		t_UmbrellaOptions *	opt,
		const char *	fnhist,
		const char *	xlabel
	)

static

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

static void copy_pullgrp_to_synthwindow	(	t_UmbrellaWindow *	synthWindow,
		t_UmbrellaWindow *	thisWindow,
		int	pullid
	)

static

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.

static void enforceEqualWeights	(	t_UmbrellaWindow *	window,
		int	nWindows
	)

static

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.

static void guessPotByIntegration	(	t_UmbrellaWindow *	window,
		int	nWindows,
		t_UmbrellaOptions *	opt,
		const char *	xlabel
	)

static

Compute initial potential by integrating the average force.

This speeds up the convergence by roughly a factor of 2

static void print_histograms	(	const char *	fnhist,
		t_UmbrellaWindow *	window,
		int	nWindows,
		int	bs_index,
		t_UmbrellaOptions *	opt,
		const char *	xlabel
	)

static

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.

static void readIntegratedAutocorrelationTimes	(	t_UmbrellaWindow *	window,
		int	nwins,
		const char *	fn
	)

static

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.

static void readPullCoordSelection	(	t_UmbrellaOptions *	opt,
		char **	fnTpr,
		int	nTpr
	)

static

Read input file for pull group selection (option -is)

TO DO: ptr=fgets(...) is never freed (small memory leak)

static void searchCumulative	(	const double	xx[],
		int	n,
		double	x,
		int *	j
	)

static

Return j such that xx[j] <= x < xx[j+1].

This is used to generate a random sequence distributed according to a histogram

static void setup_acc_wham	(	const double *	profile,
		t_UmbrellaWindow *	window,
		int	nWindows,
		t_UmbrellaOptions *	opt
	)

static

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.

static void smoothIact	(	t_UmbrellaWindow *	window,
		int	nwins,
		t_UmbrellaOptions *	opt
	)

static

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

Description

Classes

Macros

Typedefs

Enumerations

Functions

Enumeration Type Documentation

Function Documentation