Tuesday, October 8, 2013

How to use WHAM to obtain a free energy curve

WHAM is a program to be used to generate an relative free energy (FE) curve. Relative means you need specify a reference as zero.
1. Install WHAM.
2. Run your umbrella sampling (US) molecular dynamics simulations with PLUMED add-in. You should have a series of results. For example, if you want to obtain FE vs distances, you need manually specify a series of distances and run US around these distances.
3. Once you have all the results, you can use following script (a bash file):
#!/bin/bash

usfile=wham_usfile
fesout=fes.out

rm -rf $usfile $fesout
points=500  #there are 500 points/COLVAR.

for file in shortp? longp? longp??; do

tail -n $points $file/COLVAR | awk '{print $1,$2}' > $file/whamcv;
##get k constant
kappa=`grep UMBRELLA $file/plumed.dat  | head -n 1 | awk '{print $5}'`
## get umbrella center
PV=`grep UMBRELLA $file/plumed.dat  | head -n 1 | awk '{print $7}'`
echo $file/whamcv $PV $kappa >> wham_usfile
done;

first=`sort -n -k2 wham_usfile | awk '{print $2}' | head -n 1 `
last=`sort -n -k2 wham_usfile | awk '{print $2}' | tail -n 1 `
num=`wc -l wham_usfile | awk '{print $1}'`

wham $first $last $num 0.0001 500 0 $usfile $fesout
## above line shows command, min, max, number of points, tolerance, points of sampling around one distance, US input file, free energy outputfile. 

4. The output file will look like below.
#Coor           Free    +/-             Prob            +/-
9.410909        0.000000        nan     0.112293        nan
9.592727        0.228147        nan     0.076523        nan
9.774545        0.130469        nan     0.090179        nan
9.956364        0.022609        nan     0.108106        nan
10.138182       0.309783        nan     0.066711        nan
10.320000       0.200620        nan     0.080147        nan
10.501818       0.169518        nan     0.084449        nan
10.683636       0.125105        nan     0.090995        nan
10.865455       0.268023        nan     0.071562        nan
11.047273       0.013260        nan     0.109818        nan
11.229091       0.016527        nan     0.109217        nan
#Window         Free    +/-  
#0      0.000000        nan
#1      0.140395        nan
#2      -0.095729       nan
#3      -0.155083       nan
#4      -0.186455       nan
#5      -0.196517       nan
#6      -0.188120       nan
#7      -0.159148       nan
#8      -0.103205       nan
#9      -0.011637       nan
#10     0.124413        nan

We can have a FE curve using the first part. At the same time, we can obtain the probability distribution function based on the Boltzmann distribution equation (the fourth column).
The last part: these are the final F values from the wham calculation, and can be used
for computing weighted averages for properties other than the free energy.

The default unit is kcal/mol. 

No comments:

Post a Comment