Gromacs  2020.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes | Macros | Enumerations | Functions
gmx_wham.cpp File Reference
#include "gmxpre.h"
#include "config.h"
#include <cassert>
#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/mdtypes/state.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/path.h"
#include "gromacs/utility/pleasecite.h"
#include "gromacs/utility/smalloc.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 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

static t_UmbrellaWindowinitUmbrellaWindows (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 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
 
static 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...
 
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, pdo, or tpr)
 
static 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.
 
static FILE * open_pdo_pipe (const char *fn, t_UmbrellaOptions *opt, gmx_bool *bPipeOpen)
 Open a file or a pipe to a gzipped file.
 
static void pdo_close_file (FILE *fp)
 Close file or pipe.
 
static void read_pdo_files (char **fn, int nfiles, t_UmbrellaHeader *header, t_UmbrellaWindow *window, t_UmbrellaOptions *opt)
 Reading all pdo files.
 
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 read_pdo_data ( FILE *  file,
t_UmbrellaHeader header,
int  fileno,
t_UmbrellaWindow win,
t_UmbrellaOptions opt,
gmx_bool  bGetMinMax,
real mintmp,
real maxtmp 
)
static

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.

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