Potential of mean force (PMF) calculations with AMBER and WHAM

This tutorial shows how to use the AMBER programs to do umbrella sampling simulations and input the simulation data to WHAM to compute potential of mean force (PMF), a free energy map along a defined reaction coordinate. A brief description on the umbrella sampling technique can be found in the AMBER/sander manuals and the theories of the WHAM analysis (weighted histogram analysis method) can be found on Dr. Alan Grossfield's web page as well as in the original papers. This tutorial was designed to provide a simple practice to understand the basic ideas in PMF calculations, and so we chose the example shown in Dr. Grossfield's wham talk, the rotational free energy of the C-C-C-C dihedral angle in butane. Users of this tutorial are recommended to read the talk for discussions.

all-atom butane structure

1. Umbrella sampling of the torsion angle using sander in AMBER

First, we need to do a number of molecular dynamics (MD) simulations and in each simulation, a energy restraint is added to favor a particular C-C-C-C torsion value. We know a torsion angle has a range of -180 to 180 degree and so we decide to do 19 simulations by varying the favored angle from -180 to 180 with an increment of 20 degree. One of the simulation input file looks like:

This input file specifies a 500ps MD simulation with the GB (generalized Born) solvent model at temperature 300K. nmropt=1 specifies this is a simulation with a restraint and the restraint input file is chi.RST_180. In chi.RST_180, iat=1,5,8,11 specifies the atom numbers of the 4 C atoms in butane.pdb. The minimum of the restraint is placed at torsion angle = 180 (r2 = r3 = 180.0), i.e. no restraint at this angle value, and at other values, the restraint is on with a constant of rk2 = rk3 = 32.83 kcal/mol.rad2 (= 0.02 kcal/mol.deg2). r1 and r4 should be far from r2 to make sure a harmonic restraint, and we set they are 170 and 190 degree far away, respectively. The time development of the torsion angle is recorded in file chi_vs_t_180. We should prepare 19 such restraint input files by varying the values of r1, r2, r3, and r4 (zip file restraint.zip).

We run the 19 simulations by executing script run_script, and the files needed for running this script include topology file btn.top, coordinate file btn.crd, temporary input file umbrella.in, and 19 restraint files (unzip restraint.zip ).

>run_script

This outputs 19 files (zip file chi_vs_t.zip) containing the time development of the torsion angle in these simulations, which will be input to WHAM for PMF calculations.

2. Individual histogram check

Before we input the simulation data into WHAM, we should check the histograms of the torsion angles from all simulation trajectories. There should be no major shift from the expected values and no bare patches (see more discussions in Dr. Grossfield's wham talk ).

We use the program histogram_multi (source histogram_multi.c ) to compute histograms and list them in file reaction.histogram . File hist_metadatafilespecifies the number of simulations (19), the minimal (-180.0) and maximal (180.0) values of the histograms, the bin size ( 4.0), the first few lines to skip before starting counting (100), and the original data files.

>./histogram_multi hist_metadatafile reaction.histogram

The histograms are plotted as:

3. PMF calculations with WHAM

The command line to run WHAM is

>$WHAM/wham P -180.0 180.0 90 0.0001 300.0 0 metadatafile butane.pmf > wham.log

The arguments specify: the reaction coordinate has a periodicity of 360, the boundaries of the histograms are -180.0 and 180.0, the number of bins is 90, the convergence tolerance for the WHAM calculations is 0.0001, the temperature is 300K, no padding point will be printed, the umbrella sampling data are specified in file metadatafile. PMF is printed in file butane.pmf. The temporary output is recorded in wham.log.

The PMF is plotted as:

4. Effects of simulations on PMF results

The PMF results from WHAM strongly depend on how the umbrella sampling simulations are done. Even if all the individual histograms look ok (smooth, no major shift from expected, and no bare patch), you may get artifacts in the computed PMF. To give you a sense on this, we change the umbrella sampling simulations from 19 to 18 by replacing 4 simulations with restraint minima at -180.0, -160.0, 180.0 and 160.0 with 3 simulations with restraint minima at -170.0, 170.0, 10.0. The histograms and the computed PMF are plotted as:

The small turn around 170 degree is an artifact due to missing data. Solutions to obtain accurate PMFs and avoid artifacts include: know your system as much as possible, e.g. in the butane case,sampling at -180 and 180 should not be missing, do simulations as many as possible and as long as possible.

You may send your comments/questions to twang at ucdavis dot edu