BayesWHAM.py - a Bayesian implementation of the weighted histogram analysis method (WHAM) to make statistically optimal estimates of multidimensional molecular free energy surface (FES) from an ensemble of umbrella sampling simulations. The code functions in n-dimensions, natively supports uniform, symmetric Dirichlet, and Gaussian priors, and rigorously estimates uncertainties by Metropolis-Hastings sampling of the Bayes posterior. The source code may be straightforwardly modified to support arbitrary user-defined priors by augmenting two switch statements defining the MAP expression and the Metropolis acceptance criterion.

BayesWHAM_plotter.py - a utility code to plot the 1, 2, and 3-dimensional output of BayesWHAM.py

BayesReweight.py - a tool to project the unbiased FES calculated by BayesWHAM into any number of arbitrary variables beyond those in which umbrella sampling was conducted.

BayesReweight_plotter.py - a utility code to plot the 1, 2, and 3-dimensional output of BayesReweight.py

v1.0 - 08/2016

No installation required. These are Python 2.7 scripts functions that can be called and executed directly from the Python interpreter.

Full details of required inputs, input files and formatting, and output files are provided in the headers of each script. Example calculations detailing the execution steps and providing all input and output files are detailed below.

• Umbrella sampling input data

Ensemble of 18 1D umbrella sampling calculations in phi dihedral

./diala_phi_EXAMPLE/bias/harmonic_biases.txt - list of umbrella simulation ID, location of harmonic biasing center (degrees), and harmonic force constant (kJ/mol.degrees^2), one entry per line

./diala_phi_EXAMPLE/traj - trajectories of umbrella sampling variable phi in each umbrella sampling calculation

./diala_phi_EXAMPLE/hist/hist_binEdges.txt - bin edges of rectilinear grid used to create phi histograms over each umbrella sampling trajectory

./diala_phi_EXAMPLE/hist - phi histograms compiled over phi trajectories in ./diala_phi_EXAMPLE/traj

./diala_phi_EXAMPLE/traj_aux__phi_psi - trajectories of projection variables [phi,psi] into which we wish to reweight the FES recorded over the course of each umbrella sampling calculation at the same frequency as those in ./diala_phi_EXAMPLE/traj; the projection variables are arbitrary and need not contain the umbrella sampling variable

• BayesWHAM.py

Calculation of 1D FES F(phi) for alanine dipeptide from 18 1D umbrella sampling simulations at values of phi=[-170:20:170]. Phi is a backbone dihedral angle with range [-180 deg, 180 deg] and 360 deg periodicity.

We solve the generalized WHAM equations to a tolerance of 1E-15 under a Dirichlet prior with concentration parameter alpha = 2 (i.e., adding 1 pseudo count to each bin). Uncertainties are estimated by performing 1E6 rounds of Metropolis-Hastings sampling from the Bayes posterior and saving every 1E3 realizations. In practice, many more rounds of sampling should be performed until the log likelihood stabilizes indicating the Markov chain burn-in period has passed and samples are being drawn from the stationary distribution.

Copies of all files generated by this calculation are provided in ./diala_phi_EXAMPLE/BayesWHAM_OUTPUT

Execution:

dim=1;

periodicity=[1];

T=298;

harmonicBiasesFile='./diala_phi_EXAMPLE/bias/harmonic_biases.txt';

histBinEdgesFile='./diala_phi_EXAMPLE/hist/hist_binEdges.txt';

histDir='./diala_phi_EXAMPLE/hist';

tol_WHAM=1E-15;

maxIter_WHAM=1E6;

steps_MH=1E6;

saveMod_MH=1E3;

printMod_MH=1E3;

maxDpStep_MH=5E-4;

seed_MH=200184;

prior='Dirichlet';

alpha=2;

/Users/alf/Library/Enthought/Canopy_64bit/User/bin/python2.7 BayesWHAM.py $dim $periodicity $T $harmonicBiasesFile $histBinEdgesFile $histDir $tol_WHAM $maxIter_WHAM $steps_MH $saveMod_MH $printMod_MH $maxDpStep_MH $seed_MH $prior $alpha

/Users/alf/Library/Enthought/Canopy_64bit/User/bin/python2.7 BayesWHAM_plotter.py

• BayesReweight.py

Reweighting of the biased simulation data and the free energy shifts calculated from solution of the generalized WHAM equations under the Dirichlet prior to estimate F(phi,psi). Phi and psi are backbone dihedral angles with range [-180 deg, 180 deg] and 360 deg periodicity. Within this projection we anticipate good sampling in the driven variable phi but not so in psi where exploration is reliant on spontaneous thermal fluctuations.

Copies of all files generated by this calculation are provided in ./diala_phi_EXAMPLE/BayesReweight_OUTPUT

Execution:

T=298;

dim_UMB=1;

periodicity_UMB=[1];

harmonicBiasesFile_UMB='./diala_phi_EXAMPLE/bias/harmonic_biases.txt';

trajDir_UMB='./diala_phi_EXAMPLE/traj';

histBinEdgesFile_UMB='./diala_phi_EXAMPLE/hist/hist_binEdges.txt';

histDir_UMB='./diala_phi_EXAMPLE/hist';

fMAPFile_UMB='f_MAP.txt';

fMHFile_UMB='f_MH.txt';

trajDir_PROJ='./diala_phi_EXAMPLE/traj_aux__phi_psi';

histBinEdgesFile_PROJ='./diala_phi_EXAMPLE/hist_binEdges_phi_psi.txt';

/Users/alf/Library/Enthought/Canopy_64bit/User/bin/python2.7 BayesReweight.py $T $dim_UMB $periodicity_UMB $harmonicBiasesFile_UMB $trajDir_UMB $histBinEdgesFile_UMB $histDir_UMB $fMAPFile_UMB $fMHFile_UMB $trajDir_PROJ $histBinEdgesFile_PROJ

/Users/alf/Library/Enthought/Canopy_64bit/User/bin/python2.7 BayesReweight_plotter.py

• Umbrella sampling input data

Ensemble of 36 2D umbrella sampling calculations in [phi,psi] dihedrals

./diala_phi_psi_EXAMPLE/bias/harmonic_biases.txt - list of umbrella simulation ID, location of harmonic biasing centers (degrees), and harmonic force constants (kJ/mol.degrees^2), one entry per line

./diala_phi_psi_EXAMPLE/traj - trajectories of umbrella sampling variables [phi,psi] in each umbrella sampling calculation

./diala_phi_psi_EXAMPLE/hist/hist_binEdges.txt - bin edges of rectilinear grid used to create [phi,psi] histograms over each umbrella sampling trajectory

./diala_phi_psi_EXAMPLE/hist - [phi,psi] histograms compiled over [phi,psi] trajectories in ./diala_phi_psi_EXAMPLE/traj in row major order

./diala_phi_psi_EXAMPLE/traj_aux__phi_psi_theta - trajectories of projection variables [phi,psi,theta] into which we wish to reweight the FES recorded over the course of each umbrella sampling calculation at the same frequency as those in ./diala_phi_psi_EXAMPLE/traj; the projection variables are arbitrary and need not contain the umbrella sampling variables

• BayesWHAM.py

Calculation of 2D FES F(phi,psi) for alanine dipeptide from 36 2D umbrella sampling simulations at values of phi=[-170:20:170] and psi=[150:20:130]. Phi and psi are backbone dihedral angles with range [-180 deg, 180 deg] and 360 deg periodicity.

We solve the generalized WHAM equations to a tolerance of 1E-15 under a uniform prior. Uncertainties are estimated by performing 1E6 rounds of Metropolis-Hastings sampling from the Bayes posterior and saving every 1E3 realizations. In practice, many more rounds of sampling should be performed until the log likelihood stabilizes indicating the Markov chain burn-in period has passed and samples are being drawn from the stationary distribution.

Copies of all files generated by this calculation are provided in ./diala_phi_psi_EXAMPLE/BayesWHAM_OUTPUT

Execution:

dim=2;

periodicity=[1,1];

T=298;

harmonicBiasesFile='./diala_phi_psi_EXAMPLE/bias/harmonic_biases.txt';

histBinEdgesFile='./diala_phi_psi_EXAMPLE/hist/hist_binEdges.txt';

histDir='./diala_phi_psi_EXAMPLE/hist';

tol_WHAM=1E-15;

maxIter_WHAM=1E6;

steps_MH=1E6;

saveMod_MH=1E3;

printMod_MH=1E3;

maxDpStep_MH=5E-5;

seed_MH=200184;

prior='none';

alpha=0;

/Users/alf/Library/Enthought/Canopy_64bit/User/bin/python2.7 BayesWHAM.py $dim $periodicity $T $harmonicBiasesFile $histBinEdgesFile $histDir $tol_WHAM $maxIter_WHAM $steps_MH $saveMod_MH $printMod_MH $maxDpStep_MH $seed_MH $prior $alpha

/Users/alf/Library/Enthought/Canopy_64bit/User/bin/python2.7 BayesWHAM_plotter.py

• BayesReweight.py

Reweighting of the biased simulation data and the free energy shifts calculated from solution of the generalized WHAM equations under the uniform prior to estimate F(phi,psi,theta). Phi, psi, and theta are backbone dihedral angles with range [-180 deg, 180 deg] and 360 deg periodicity. Within this projection we anticipate good sampling in the driven variables phi and psi but not so in theta where exploration is reliant on spontaneous thermal fluctuations.

Copies of all files generated by this calculation are provided in ./diala_phi_psi_EXAMPLE/BayesReweight_OUTPUT

Execution:

T=298;

dim_UMB=2;

periodicity_UMB=[1,1];

harmonicBiasesFile_UMB='./diala_phi_psi_EXAMPLE/bias/harmonic_biases.txt';

trajDir_UMB='./diala_phi_psi_EXAMPLE/traj';

histBinEdgesFile_UMB='./diala_phi_psi_EXAMPLE/hist/hist_binEdges.txt';

histDir_UMB='./diala_phi_psi_EXAMPLE/hist';

fMAPFile_UMB='f_MAP.txt';

fMHFile_UMB='f_MH.txt';

trajDir_PROJ='./diala_phi_psi_EXAMPLE/traj_aux__phi_psi_theta';

histBinEdgesFile_PROJ='./diala_phi_psi_EXAMPLE/hist_binEdges_phi_psi_theta.txt';

/Users/alf/Library/Enthought/Canopy_64bit/User/bin/python2.7 BayesReweight.py $T $dim_UMB $periodicity_UMB $harmonicBiasesFile_UMB $trajDir_UMB $histBinEdgesFile_UMB $histDir_UMB $fMAPFile_UMB $fMHFile_UMB $trajDir_PROJ $histBinEdgesFile_PROJ

/Users/alf/Library/Enthought/Canopy_64bit/User/bin/python2.7 BayesReweight_plotter.py

If you use the BayesWHAM.py and BayesReweight.py codes, please read and cite:

A.L. Ferguson "BayesWHAM: A Bayesian approach for free energy estimation, reweighting, and uncertainty quantification in the weighted histogram analysis method" J. Comput. Chem. 38 18 1583-1605 (2017)

Andrew Ferguson, PhD

Assistant Professor of Materials Science and Engineering

University of Illinois at Urbana-Champaign

[email protected]

http://ferguson.matse.illinois.edu