1. Andrew Ferguson
  2. bayeswham_matlab_v2

Overview

HTTPS SSH

README

Change Log and Usage

BayesWHAM_v2 is a generalization of BayesWHAM_v1 that adds the capacity to (i) incorporate biased simulation data collected at any temperature and pressure in the NPT ensemble, and (ii) reweight these data to estimate the unbiased free energy landscape FES at any T and P.

  • The code is designed to function in the NPT ensemble.

  • Reweighting extrapolatively or interpolatively to [T,P] state points far from those at which data was collected can result in convergence failure of the WHAM iteration.

  • Execution is typcically more CPU and RAM intensive than BayesWHAM_v1 due to necessity to collect histograms in the two additional dimensions U and V to permit reweighting in T and P.

  • For simulation data collected at a single T and P state point, and construction of the landscape at that T and P, the additional functionality of BayesWHAM_v2 is not required. Accordingly, (i) BayesWHAM_v1 may be used or (ii) BayesWHAM_v2 may be used regardless, and the unnecessary cost associated with E and V binning of the data -- since no biasing is introduced in E or V -- eliminated by binning the data into a single U bin and a single V bin. See Example 1 below. (Note that reweighting to a different T or P state with a single U and V bin will introduce large errors due to poor resolution in U and V reweighting.)

BayesReweight_v2 is a generalization of BayesReweight_v1 designed for compatibility with BayesWHAM_v2.

BayesWHAM_plotter and BayesReweight_plotter remain unchanged from v1.

Theory

The theoretical underpinnings of the v2 generalization are described in BayesWHAM_addendum.pdf that should be read as an addendum to the published work 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).

Synopsis

BayesWHAM_v2.m - 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 conducted in NPT ensemble. The FES can be constructed from data collected at different T and P, with biasing in multidimensional collective variables psi, and requires a record of the trajectories in potential energy U, volume V, and biasing variables psi. 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.m - a utility code to plot the 1, 2, and 3-dimensional output of BayesWHAM.m

BayesReweight_v2.m - 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.m - a utility code to plot the 1, 2, and 3-dimensional output of BayesReweight.m

Version

v2 - 06/2017

v1 - 08/2016

Installation

No installation required. These are MATLAB functions that can be called and executed directly from the MATLAB prompt.

I/O

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.

Example 1 - alanine dipeptide 2D umbrella sampling in [phi,psi] dihedrals at fixed T and P

  • Umbrella sampling input data

Ensemble of 36 2D umbrella sampling calculations in [phi,psi] dihedrals at T=300 K and P=1 bar

./diala_phi_psi_EXAMPLE1/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_EXAMPLE1/bias/T_P.txt - list of umbrella simulation ID, temperature / K, and pressure / bar, one entry per line

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

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

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

./diala_phi_psi_EXAMPLE1/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_EXAMPLE1/traj; the projection variables are arbitrary and need not contain the umbrella sampling variables

  • BayesWHAM_v2.m

[1]

Calculation of 2D FES F(phi,psi) at T=300 K and P=1 bar (i.e., at same T & P at which all biased simulations were conducted) where phi and psi are backbone dihedral angles with range [-180 deg, 180 deg] and 360 deg periodicity.

Histogramming has been conducted binning phi and psi every 10 degrees, but E and V have each been histogrammed into a single bin. Since landscape is being computed at the T=300 K and P=1 bar at wich all simulations were conducted, the E and V biasing factors are precisely zero and we can accelerate calculation by effectively declining to bin in E and V.

We solve the generalized WHAM equations to a tolerance of 1E-15 under a uniform prior. Uncertainties are estimated by performing 2E5 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_EXAMPLE1/BayesWHAM_v2_OUTPUT

Execution:

>> dim=2;

>> periodicity=[1,1];

>> T=300;

>> P=1;

>> TPFile='./diala_phi_psi_EXAMPLE1/bias/T_P.txt';

>> harmonicBiasesFile='./diala_phi_psi_EXAMPLE1/bias/harmonic_biases.txt';

>> histBinEdgesFile='./diala_phi_psi_EXAMPLE1/hist/hist_binEdges.txt';

>> histDir='./diala_phi_psi_EXAMPLE1/hist';

>> tol_WHAM=1E-15;

>> maxIter_WHAM=1E6;

>> steps_MH=2E5;

>> saveMod_MH=1E3;

>> printMod_MH=1E3;

>> maxDpStep_MH=5E-5;

>> seed_MH=200184;

>> prior='none';

>> alpha=0;

>> BayesWHAM_v2(dim,periodicity,T,P,TPFile,harmonicBiasesFile,histBinEdgesFile,histDir,tol_WHAM,maxIter_WHAM,steps_MH,saveMod_MH,printMod_MH,maxDpStep_MH,seed_MH,prior,alpha)

>> BayesWHAM_plotter

  • BayesWHAM_v1.m

[1]

Since landscape is being computed at the T=300 K and P=1 bar at wich all simulations were conducted and binning in E and V was into a single bin each, we can run the previous version of BayesWHAM on these data to compare results

We solve the generalized WHAM equations to a tolerance of 1E-15 under a uniform prior. Uncertainties are estimated by performing 2E5 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_EXAMPLE1/BayesWHAM_v1_OUTPUT

Output files are identical to those generated by BayesWHAM_v2.m in ./diala_phi_psi_EXAMPLE1/BayesWHAM_v2_OUTPUT

Execution:

>> dim=2;

>> periodicity=[1,1];

>> T=300;

>> harmonicBiasesFile='./diala_phi_psi_EXAMPLE1/bias/harmonic_biases.txt';

>> histBinEdgesFile='./diala_phi_psi_EXAMPLE1/hist/hist_binEdges_noEnoV.txt';

>> histDir='./diala_phi_psi_EXAMPLE1/hist';

>> tol_WHAM=1E-15;

>> maxIter_WHAM=1E6;

>> steps_MH=2E5;

>> saveMod_MH=1E3;

>> printMod_MH=1E3;

>> maxDpStep_MH=5E-5;

>> seed_MH=200184;

>> prior='none';

>> alpha=0;

>> BayesWHAM_v1(dim,periodicity,T,harmonicBiasesFile,histBinEdgesFile,histDir,tol_WHAM,maxIter_WHAM,steps_MH,saveMod_MH,printMod_MH,maxDpStep_MH,seed_MH,prior,alpha)

>> BayesWHAM_plotter

Example 2 - alanine dipeptide 2D umbrella sampling in [phi,psi] dihedrals at fixed T and P

  • Umbrella sampling input data

Ensemble of 36 2D umbrella sampling calculations in [phi,psi] dihedrals at T=300 K and P=1 bar

./diala_phi_psi_EXAMPLE2/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_EXAMPLE2/bias/T_P.txt - list of umbrella simulation ID, temperature / K, and pressure / bar, one entry per line

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

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

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

./diala_phi_psi_EXAMPLE2/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_EXAMPLE2/traj; the projection variables are arbitrary and need not contain the umbrella sampling variables

  • BayesWHAM_v2.m

[1]

Calculation of 2D FES F(phi,psi) at T=300 K and P=1 bar (i.e., at same T & P at which all biased simulations were conducted) where phi and psi are backbone dihedral angles with range [-180 deg, 180 deg] and 360 deg periodicity.

Histogramming has been conducted binning phi and psi every 10 degrees, E every 1000 kJ/mol, and V every 1 nm^3. Since landscape is being constructed at same T and P at which all biased simulations were conducted, the E and V binning is irrelevant since the biasing factors in E and V are identically zero but does affect execution speed.

We solve the generalized WHAM equations to a tolerance of 1E-15 under a uniform prior. Uncertainties are estimated by performing 2E5 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_EXAMPLE2/BayesWHAM_v2_OUTPUT_T300P1

Output files [betaF_MAP.txt, f_MAP.txt, hist_binCenters.txt, hist_binWidths.txt, p_MAP.txt, pdf_MAP.txt] are identical to those generated by BayesWHAM_v2.m in ./diala_phi_psi_EXAMPLE1/BayesWHAM_v2_OUTPUT and ./diala_phi_psi_EXAMPLE1/BayesWHAM_v1_OUTPUT.

Output files [logL_MAP.txt, all _MH, all _augUV] differ due to different histogram binning (i.e., the binning of E and V in the present case) that gives rise to differnt value of likelihood and sequence of random numbers in MH sampling.

Execution:

>> dim=2;

>> periodicity=[1,1];

>> T=300;

>> P=1;

>> TPFile='./diala_phi_psi_EXAMPLE2/bias/T_P.txt';

>> harmonicBiasesFile='./diala_phi_psi_EXAMPLE2/bias/harmonic_biases.txt';

>> histBinEdgesFile='./diala_phi_psi_EXAMPLE2/hist/hist_binEdges.txt';

>> histDir='./diala_phi_psi_EXAMPLE2/hist';

>> tol_WHAM=1E-15;

>> maxIter_WHAM=1E6;

>> steps_MH=2E5;

>> saveMod_MH=1E3;

>> printMod_MH=1E3;

>> maxDpStep_MH=5E-5;

>> seed_MH=200184;

>> prior='none';

>> alpha=0;

>> BayesWHAM_v2(dim,periodicity,T,P,TPFile,harmonicBiasesFile,histBinEdgesFile,histDir,tol_WHAM,maxIter_WHAM,steps_MH,saveMod_MH,printMod_MH,maxDpStep_MH,seed_MH,prior,alpha)

>> BayesWHAM_plotter

[2]

Calculation of 2D FES F(phi,psi) at T=310 K and P=50 bar where phi and psi are backbone dihedral angles with range [-180 deg, 180 deg] and 360 deg periodicity.

Extrapolations / interpolations to state points far from those at which data were collected can lead to large errors.

Histogramming has been conducted binning phi and psi every 10 degrees, E every 1000 kJ/mol, and V every 1 nm^3. Since landscape is being constructed at different T and P from that at which all biased simulations were conducted, the E and V binning affects reweighting accuracy and execution speed. In practice, the E and V binning is probably too coarse, and the robustness of the results to the bin width should be assessed.

We solve the generalized WHAM equations to a tolerance of 1E-15 under a uniform prior. Uncertainties are estimated by performing 2E5 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_EXAMPLE2/BayesWHAM_v2_OUTPUT_T310P50

Execution:

>> dim=2;

>> periodicity=[1,1];

>> T=310;

>> P=50;

>> TPFile='./diala_phi_psi_EXAMPLE2/bias/T_P.txt';

>> harmonicBiasesFile='./diala_phi_psi_EXAMPLE2/bias/harmonic_biases.txt';

>> histBinEdgesFile='./diala_phi_psi_EXAMPLE2/hist/hist_binEdges.txt';

>> histDir='./diala_phi_psi_EXAMPLE2/hist';

>> tol_WHAM=1E-15;

>> maxIter_WHAM=1E6;

>> steps_MH=2E5;

>> saveMod_MH=1E3;

>> printMod_MH=1E3;

>> maxDpStep_MH=5E-5;

>> seed_MH=200184;

>> prior='none';

>> alpha=0;

>> BayesWHAM_v2(dim,periodicity,T,P,TPFile,harmonicBiasesFile,histBinEdgesFile,histDir,tol_WHAM,maxIter_WHAM,steps_MH,saveMod_MH,printMod_MH,maxDpStep_MH,seed_MH,prior,alpha)

>> BayesWHAM_plotter

[3]

Calculation of 2D FES F(phi,psi) at T=290 K and P=1 bar where phi and psi are backbone dihedral angles with range [-180 deg, 180 deg] and 360 deg periodicity.

Extrapolations / interpolations to state points far from those at which data were collected can lead to large errors.

Histogramming has been conducted binning phi and psi every 10 degrees, E every 1000 kJ/mol, and V every 1 nm^3. Since landscape is being constructed at different T and P from that at which all biased simulations were conducted, the E and V binning affects reweighting accuracy and execution speed. In practice, the E and V binning is probably too coarse, and the robustness of the results to the bin width should be assessed.

We solve the generalized WHAM equations to a tolerance of 1E-15 under a uniform prior. Uncertainties are estimated by performing 2E5 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_EXAMPLE2/BayesWHAM_v2_OUTPUT_T290P1

Execution:

>> dim=2;

>> periodicity=[1,1];

>> T=290;

>> P=1;

>> TPFile='./diala_phi_psi_EXAMPLE2/bias/T_P.txt';

>> harmonicBiasesFile='./diala_phi_psi_EXAMPLE2/bias/harmonic_biases.txt';

>> histBinEdgesFile='./diala_phi_psi_EXAMPLE2/hist/hist_binEdges.txt';

>> histDir='./diala_phi_psi_EXAMPLE2/hist';

>> tol_WHAM=1E-15;

>> maxIter_WHAM=1E6;

>> steps_MH=2E5;

>> saveMod_MH=1E3;

>> printMod_MH=1E3;

>> maxDpStep_MH=5E-5;

>> seed_MH=200184;

>> prior='none';

>> alpha=0;

>> BayesWHAM_v2(dim,periodicity,T,P,TPFile,harmonicBiasesFile,histBinEdgesFile,histDir,tol_WHAM,maxIter_WHAM,steps_MH,saveMod_MH,printMod_MH,maxDpStep_MH,seed_MH,prior,alpha)

>> BayesWHAM_plotter

Example 3 - alanine dipeptide 2D umbrella sampling in [phi,psi] dihedrals at various T and P

  • Umbrella sampling input data

Ensemble of 108 2D umbrella sampling calculations in [phi,psi] dihedrals at (i) T=300 K and P=1 bar (1-36), (ii) T=310 K and P=50 bar (37-72), and (iii) T=320 K and P=100 bar (73-108).

./diala_phi_psi_EXAMPLE3/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_EXAMPLE3/bias/T_P.txt - list of umbrella simulation ID, temperature / K, and pressure / bar, one entry per line

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

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

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

./diala_phi_psi_EXAMPLE3/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_EXAMPLE3/traj; the projection variables are arbitrary and need not contain the umbrella sampling variables

  • BayesWHAM_v2.m

[1]

Phi and psi are backbone dihedral angles with range [-180 deg, 180 deg] and 360 deg periodicity.

Histogramming has been conducted binning phi and psi every 20 degrees, E every 200 kJ/mol, and V every 0.5 nm^3. E and V binning affects accuracy and speed of reweighting.

We solve the generalized WHAM equations to a tolerance of 1E-15 under a uniform prior. Uncertainties are estimated by performing 2E5 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.

Calculation of 2D FES F(phi,psi) at T = 300 K and P = 1 bar from data collected at T = 300 K and P = 1 bar (1-36).

Copies of all files generated by this calculation are provided in ./diala_phi_psi_EXAMPLE3/BayesWHAM_v2_OUTPUT_T300P1_1to36

Execution:

>> dim=2;

>> periodicity=[1,1];

>> T=300;

>> P=1;

>> TPFile='./diala_phi_psi_EXAMPLE3/bias/T_P_1to36.txt';

>> harmonicBiasesFile='./diala_phi_psi_EXAMPLE3/bias/harmonic_biases_1to36.txt';

>> histBinEdgesFile='./diala_phi_psi_EXAMPLE3/hist/hist_binEdges.txt';

>> histDir='./diala_phi_psi_EXAMPLE3/hist';

>> tol_WHAM=1E-15;

>> maxIter_WHAM=1E6;

>> steps_MH=2E5;

>> saveMod_MH=1E3;

>> printMod_MH=1E3;

>> maxDpStep_MH=5E-5;

>> seed_MH=200184;

>> prior='none';

>> alpha=0;

>> BayesWHAM_v2(dim,periodicity,T,P,TPFile,harmonicBiasesFile,histBinEdgesFile,histDir,tol_WHAM,maxIter_WHAM,steps_MH,saveMod_MH,printMod_MH,maxDpStep_MH,seed_MH,prior,alpha)

>> BayesWHAM_plotter

[2]

Calculation of 2D FES F(phi,psi) at T = 300 K and P = 1 bar from data collected at T = 300 K and P = 1 bar (1-36) and T = 310 K and P = 50 bar (37-72).

Copies of all files generated by this calculation are provided in ./diala_phi_psi_EXAMPLE3/BayesWHAM_v2_OUTPUT_T300P1_1to72

Execution:

>> dim=2;

>> periodicity=[1,1];

>> T=300;

>> P=1;

>> TPFile='./diala_phi_psi_EXAMPLE3/bias/T_P_1to72.txt';

>> harmonicBiasesFile='./diala_phi_psi_EXAMPLE3/bias/harmonic_biases_1to72.txt';

>> histBinEdgesFile='./diala_phi_psi_EXAMPLE3/hist/hist_binEdges.txt';

>> histDir='./diala_phi_psi_EXAMPLE3/hist';

>> tol_WHAM=1E-15;

>> maxIter_WHAM=1E6;

>> steps_MH=2E5;

>> saveMod_MH=1E3;

>> printMod_MH=1E3;

>> maxDpStep_MH=5E-5;

>> seed_MH=200184;

>> prior='none';

>> alpha=0;

>> BayesWHAM_v2(dim,periodicity,T,P,TPFile,harmonicBiasesFile,histBinEdgesFile,histDir,tol_WHAM,maxIter_WHAM,steps_MH,saveMod_MH,printMod_MH,maxDpStep_MH,seed_MH,prior,alpha)

>> BayesWHAM_plotter

[3]

Calculation of 2D FES F(phi,psi) at T = 300 K and P = 1 bar from data collected at T = 300 K and P = 1 bar (1-36), T = 310 K and P = 50 bar (37-72), and T = 320 K and P = 100 bar (73-108).

Copies of all files generated by this calculation are provided in ./diala_phi_psi_EXAMPLE3/BayesWHAM_v2_OUTPUT_T300P1_1to108

Execution:

>> dim=2;

>> periodicity=[1,1];

>> T=300;

>> P=1;

>> TPFile='./diala_phi_psi_EXAMPLE3/bias/T_P_1to108.txt';

>> harmonicBiasesFile='./diala_phi_psi_EXAMPLE3/bias/harmonic_biases_1to108.txt';

>> histBinEdgesFile='./diala_phi_psi_EXAMPLE3/hist/hist_binEdges.txt';

>> histDir='./diala_phi_psi_EXAMPLE3/hist';

>> tol_WHAM=1E-15;

>> maxIter_WHAM=1E6;

>> steps_MH=2E5;

>> saveMod_MH=1E3;

>> printMod_MH=1E3;

>> maxDpStep_MH=5E-5;

>> seed_MH=200184;

>> prior='none';

>> alpha=0;

>> BayesWHAM_v2(dim,periodicity,T,P,TPFile,harmonicBiasesFile,histBinEdgesFile,histDir,tol_WHAM,maxIter_WHAM,steps_MH,saveMod_MH,printMod_MH,maxDpStep_MH,seed_MH,prior,alpha)

>> BayesWHAM_plotter

[4]

Calculation of 2D FES F(phi,psi) at T = 320 K and P = 50 bar from data collected at T = 300 K and P = 1 bar (1-36).

Copies of all files generated by this calculation are provided in ./diala_phi_psi_EXAMPLE3/BayesWHAM_v2_OUTPUT_T320P50_1to36

Execution:

>> dim=2;

>> periodicity=[1,1];

>> T=320;

>> P=50;

>> TPFile='./diala_phi_psi_EXAMPLE3/bias/T_P_1to36.txt';

>> harmonicBiasesFile='./diala_phi_psi_EXAMPLE3/bias/harmonic_biases_1to36.txt';

>> histBinEdgesFile='./diala_phi_psi_EXAMPLE3/hist/hist_binEdges.txt';

>> histDir='./diala_phi_psi_EXAMPLE3/hist';

>> tol_WHAM=1E-15;

>> maxIter_WHAM=1E6;

>> steps_MH=2E5;

>> saveMod_MH=1E3;

>> printMod_MH=1E3;

>> maxDpStep_MH=5E-5;

>> seed_MH=200184;

>> prior='none';

>> alpha=0;

>> BayesWHAM_v2(dim,periodicity,T,P,TPFile,harmonicBiasesFile,histBinEdgesFile,histDir,tol_WHAM,maxIter_WHAM,steps_MH,saveMod_MH,printMod_MH,maxDpStep_MH,seed_MH,prior,alpha)

>> BayesWHAM_plotter

[5]

Calculation of 2D FES F(phi,psi) at T = 320 K and P = 50 bar from data collected at T = 300 K and P = 1 bar (1-36) and T = 310 K and P = 50 bar (37-72).

Copies of all files generated by this calculation are provided in ./diala_phi_psi_EXAMPLE3/BayesWHAM_v2_OUTPUT_T320P50_1to72

Execution:

>> dim=2;

>> periodicity=[1,1];

>> T=320;

>> P=50;

>> TPFile='./diala_phi_psi_EXAMPLE3/bias/T_P_1to72.txt';

>> harmonicBiasesFile='./diala_phi_psi_EXAMPLE3/bias/harmonic_biases_1to72.txt';

>> histBinEdgesFile='./diala_phi_psi_EXAMPLE3/hist/hist_binEdges.txt';

>> histDir='./diala_phi_psi_EXAMPLE3/hist';

>> tol_WHAM=1E-15;

>> maxIter_WHAM=1E6;

>> steps_MH=2E5;

>> saveMod_MH=1E3;

>> printMod_MH=1E3;

>> maxDpStep_MH=5E-5;

>> seed_MH=200184;

>> prior='none';

>> alpha=0;

>> BayesWHAM_v2(dim,periodicity,T,P,TPFile,harmonicBiasesFile,histBinEdgesFile,histDir,tol_WHAM,maxIter_WHAM,steps_MH,saveMod_MH,printMod_MH,maxDpStep_MH,seed_MH,prior,alpha)

>> BayesWHAM_plotter

[6]

Calculation of 2D FES F(phi,psi) at T = 320 K and P = 50 bar from data collected at T = 300 K and P = 1 bar (1-36), T = 310 K and P = 50 bar (37-72), and T = 320 K and P = 100 bar (73-108).

Copies of all files generated by this calculation are provided in ./diala_phi_psi_EXAMPLE3/BayesWHAM_v2_OUTPUT_T320P50_1to108

Execution:

>> dim=2;

>> periodicity=[1,1];

>> T=320;

>> P=50;

>> TPFile='./diala_phi_psi_EXAMPLE3/bias/T_P_1to108.txt';

>> harmonicBiasesFile='./diala_phi_psi_EXAMPLE3/bias/harmonic_biases_1to108.txt';

>> histBinEdgesFile='./diala_phi_psi_EXAMPLE3/hist/hist_binEdges.txt';

>> histDir='./diala_phi_psi_EXAMPLE3/hist';

>> tol_WHAM=1E-15;

>> maxIter_WHAM=1E6;

>> steps_MH=2E5;

>> saveMod_MH=1E3;

>> printMod_MH=1E3;

>> maxDpStep_MH=5E-5;

>> seed_MH=200184;

>> prior='none';

>> alpha=0;

>> BayesWHAM_v2(dim,periodicity,T,P,TPFile,harmonicBiasesFile,histBinEdgesFile,histDir,tol_WHAM,maxIter_WHAM,steps_MH,saveMod_MH,printMod_MH,maxDpStep_MH,seed_MH,prior,alpha)

>> BayesWHAM_plotter

[7]

Calculation of 2D FES F(phi,psi) at T = 285 K and P = 300 bar from data collected at T = 300 K and P = 1 bar (1-36), T = 310 K and P = 50 bar (37-72), and T = 320 K and P = 100 bar (73-108), under a Dirichlet prior with concentration parameter alpha = 10.

Copies of all files generated by this calculation are provided in ./diala_phi_psi_EXAMPLE3/BayesWHAM_v2_OUTPUT_T295P75Dirichlet10_1to108

Execution:

>> dim=2;

>> periodicity=[1,1];

>> T=295;

>> P=75;

>> TPFile='./diala_phi_psi_EXAMPLE3/bias/T_P_1to108.txt';

>> harmonicBiasesFile='./diala_phi_psi_EXAMPLE3/bias/harmonic_biases_1to108.txt';

>> histBinEdgesFile='./diala_phi_psi_EXAMPLE3/hist/hist_binEdges.txt';

>> histDir='./diala_phi_psi_EXAMPLE3/hist';

>> tol_WHAM=1E-15;

>> maxIter_WHAM=1E6;

>> steps_MH=2E5;

>> saveMod_MH=1E3;

>> printMod_MH=1E3;

>> maxDpStep_MH=5E-5;

>> seed_MH=200184;

>> prior='Dirichlet';

>> alpha=10;

>> BayesWHAM_v2(dim,periodicity,T,P,TPFile,harmonicBiasesFile,histBinEdgesFile,histDir,tol_WHAM,maxIter_WHAM,steps_MH,saveMod_MH,printMod_MH,maxDpStep_MH,seed_MH,prior,alpha)

>> BayesWHAM_plotter

Example 4 - alanine dipeptide 2D umbrella sampling in [phi,psi] dihedrals at various T and P plus reweighting

  • Umbrella sampling input data

Ensemble of 72 2D umbrella sampling calculations in [phi,psi] dihedrals at T=300 K and P=1 bar (1-36) and T=310 K and P=50 bar (37-72).

./diala_phi_psi_EXAMPLE4/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_EXAMPLE4/bias/T_P.txt - list of umbrella simulation ID, temperature / K, and pressure / bar, one entry per line

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

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

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

./diala_phi_psi_EXAMPLE4/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_EXAMPLE4/traj; the projection variables are arbitrary and need not contain the umbrella sampling variables

  • BayesWHAM_v2.m

[1]

Phi and psi are backbone dihedral angles with range [-180 deg, 180 deg] and 360 deg periodicity.

Histogramming has been conducted binning phi and psi every 20 degrees, E every 200 kJ/mol, and V every 0.5 nm^3. E and V binning affects accuracy and speed of reweighting.

We solve the generalized WHAM equations to a tolerance of 1E-15 under a uniform prior. Uncertainties are estimated by performing 2E5 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.

Calculation of 2D FES F(phi,psi) at T = 305 K and P = 25 bar from data collected at T = 300 K and P = 1 bar (1-36) and T = 310 K and P = 50 bar (37-72).

Copies of all files generated by this calculation are provided in ./diala_phi_psi_EXAMPLE4/BayesWHAM_v2_OUTPUT

Execution:

>> dim=2;

>> periodicity=[1,1];

>> T=305;

>> P=25;

>> TPFile='./diala_phi_psi_EXAMPLE4/bias/T_P.txt';

>> harmonicBiasesFile='./diala_phi_psi_EXAMPLE4/bias/harmonic_biases.txt';

>> histBinEdgesFile='./diala_phi_psi_EXAMPLE4/hist/hist_binEdges.txt';

>> histDir='./diala_phi_psi_EXAMPLE4/hist';

>> tol_WHAM=1E-15;

>> maxIter_WHAM=1E6;

>> steps_MH=2E5;

>> saveMod_MH=1E3;

>> printMod_MH=1E3;

>> maxDpStep_MH=5E-5;

>> seed_MH=200184;

>> prior='none';

>> alpha=0;

>> BayesWHAM_v2(dim,periodicity,T,P,TPFile,harmonicBiasesFile,histBinEdgesFile,histDir,tol_WHAM,maxIter_WHAM,steps_MH,saveMod_MH,printMod_MH,maxDpStep_MH,seed_MH,prior,alpha)

>> BayesWHAM_plotter

  • BayesReweight_v2.m

[1]

Reweighting then performed into 2D FES F(phi,psi) at T = 305 K and P = 25 bar with [phi,psi] bins at 20 deg intervals.

Copies of all files generated by this calculation are provided in ./diala_phi_psi_EXAMPLE4/BayesReweight_v2_T305P25_phiPsi_OUTPUT

Execution:

>> T=305;

>> P=25;

>> dim_UMB=2;

>> periodicity_UMB=[1,1];

>> TPFile_UMB='./diala_phi_psi_EXAMPLE4/bias/T_P.txt';

>> harmonicBiasesFile_UMB='./diala_phi_psi_EXAMPLE4/bias/harmonic_biases.txt';

>> trajDir_UMB='./diala_phi_psi_EXAMPLE4/traj';

>> histBinEdgesFile_UMB='./diala_phi_psi_EXAMPLE4/hist/hist_binEdges.txt';

>> histDir_UMB='./diala_phi_psi_EXAMPLE4/hist';

>> fMAPFile_UMB='f_MAP.txt';

>> fMHFile_UMB='f_MH.txt';

>> trajDir_PROJ='./diala_phi_psi_EXAMPLE4/traj_aux__phi_psi';

>> histBinEdgesFile_PROJ='./diala_phi_psi_EXAMPLE4/hist_binEdges_phi_psi.txt';

>> BayesReweight_v2(T,P,dim_UMB,periodicity_UMB,TPFile_UMB,harmonicBiasesFile_UMB,trajDir_UMB,histBinEdgesFile_UMB,histDir_UMB,fMAPFile_UMB,fMHFile_UMB,trajDir_PROJ,histBinEdgesFile_PROJ)

>> BayesReweight_plotter

[2]

Reweighting then performed into 2D FES F(phi,psi) at T = 315 K and P = 40 bar with [phi,psi] bins at 10 deg intervals.

Copies of all files generated by this calculation are provided in ./diala_phi_psi_EXAMPLE4/BayesReweight_v2_T315P40_phiPsi_OUTPUT

Execution:

>> T=315;

>> P=40;

>> dim_UMB=2;

>> periodicity_UMB=[1,1];

>> TPFile_UMB='./diala_phi_psi_EXAMPLE4/bias/T_P.txt';

>> harmonicBiasesFile_UMB='./diala_phi_psi_EXAMPLE4/bias/harmonic_biases.txt';

>> trajDir_UMB='./diala_phi_psi_EXAMPLE4/traj';

>> histBinEdgesFile_UMB='./diala_phi_psi_EXAMPLE4/hist/hist_binEdges.txt';

>> histDir_UMB='./diala_phi_psi_EXAMPLE4/hist';

>> fMAPFile_UMB='f_MAP.txt';

>> fMHFile_UMB='f_MH.txt';

>> trajDir_PROJ='./diala_phi_psi_EXAMPLE4/traj_aux__phi_psi';

>> histBinEdgesFile_PROJ='./diala_phi_psi_EXAMPLE4/hist_binEdges_phi_psi_10.txt';

>> BayesReweight_v2(T,P,dim_UMB,periodicity_UMB,TPFile_UMB,harmonicBiasesFile_UMB,trajDir_UMB,histBinEdgesFile_UMB,histDir_UMB,fMAPFile_UMB,fMHFile_UMB,trajDir_PROJ,histBinEdgesFile_PROJ)

>> BayesReweight_plotter

[3]

Reweighting then performed into 3D FES F(phi,psi,theta) at T = 305 K and P = 25 bar with [phi,psi,theta] bins at 20 deg intervals.

Copies of all files generated by this calculation are provided in ./diala_phi_psi_EXAMPLE4/BayesReweight_v2_T305P25_phiPsiTheta_OUTPUT

Execution:

>> T=305;

>> P=25;

>> dim_UMB=2;

>> periodicity_UMB=[1,1];

>> TPFile_UMB='./diala_phi_psi_EXAMPLE4/bias/T_P.txt';

>> harmonicBiasesFile_UMB='./diala_phi_psi_EXAMPLE4/bias/harmonic_biases.txt';

>> trajDir_UMB='./diala_phi_psi_EXAMPLE4/traj';

>> histBinEdgesFile_UMB='./diala_phi_psi_EXAMPLE4/hist/hist_binEdges.txt';

>> histDir_UMB='./diala_phi_psi_EXAMPLE4/hist';

>> fMAPFile_UMB='f_MAP.txt';

>> fMHFile_UMB='f_MH.txt';

>> trajDir_PROJ='./diala_phi_psi_EXAMPLE4/traj_aux__phi_psi_theta';

>> histBinEdgesFile_PROJ='./diala_phi_psi_EXAMPLE4/hist_binEdges_phi_psi_theta.txt';

>> BayesReweight_v2(T,P,dim_UMB,periodicity_UMB,TPFile_UMB,harmonicBiasesFile_UMB,trajDir_UMB,histBinEdgesFile_UMB,histDir_UMB,fMAPFile_UMB,fMHFile_UMB,trajDir_PROJ,histBinEdgesFile_PROJ)

>> BayesReweight_plotter

Citing

If you use the BayesWHAM and BayesReweight 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)

Contact

Andrew Ferguson, PhD

Assistant Professor of Materials Science and Engineering

University of Illinois at Urbana-Champaign

alf@illinois.edu

http://ferguson.matse.illinois.edu