Weighted Histogram Analysis Method - WHAM

WHAM is used extensively to calculate the free energy of a system along a properly chosen reaction coordinate in molecular dynamics simulations or monte carlo simulations. To be able to use WHAM you need to perform multiple simulations that cover the entire reaction coordinate to sample possible configurations. The choice of reaction coordinate constitutes a length discussion and is a topic of ongoing debate. So I encourage the reader to look it up elsewhere. After you have chosen a reaction coordinate and performed Umbrella Sampling or Free Energy Perturbation calculations, histograms corresponding to every window need to be extracted.

For the simple case of dragging a methane molecule from water into a lipid bilayer and out into water again, z axis is a good reaction coordinate if the lipid bilayer is on the x-y plane. With this set up, the z axis is split into discrete windows and the center of mass (COM) of the methane molecule is constrained (usually harmonic) to the center of every window. The length of windows and force constant of the constraint are chose such that the distribution of the z coordinate of the COM of methane overlaps with the distribution of the adjacent windows. This is a necessary prerequisite for using WHAM. Following the figure below, for a 60A long reaction coordinate a choice of 121 windows (each 0.5A and including both ends) will be plenty to ensure overlap of z coordinate distributions.

Collect z coordinate distributions and put them 121 files labeled window-001 to window-121.
Create a new file input.in that contains the following lines;

./window-001 -30.0 xx
./window-002 -29.5 xx
...
./window-120 +29.5 xx
./window-121 +30.0 xx

Here the first column is the file name, second is the mean value of that window and the third column needs to be filled with an integer but is irrelevant and not used by Alan Grossfield's implmentation of WHAM. Then run the following command

./wham 1 11 10 0.001 300 0 input.in freefile

where the inputs are

./wham hist.min hist.max num.bins tol Temp numpad infile outfile

Figure Source: http://www.utdallas.edu/~son051000/comp/FreeE.pdf

The output contains a plot for the potential of mean force (PMF). Pay particular attention to units since they depend on your input. WHAM's documentation can be found here [http://membrane.urmc.rochester.edu/sites/default/files/wham/doc.html].

DNA simulation set-up with AMBER

I have been working on the structure and dynamics of single and double stranded DNA for some time now and below is a write up of an easy way to setup the simulations using the AMBER suite of software. Along with primary executables sander, pmemd, tleap and ptraj, AMBER comes with a number of other useful tools. One such tool is the nucleic acid builder (or nab). To create a 10mer poly-Adenosine poly-Thymidine, put the following text in a file called nuc.nab;

molecule m;
m = fd_helix( "abdna", "aaaaaaaaaa", "dna" );
putpdb( "nuc.pdb", m, "-wwpdb");

Then run these in the terminal

nab nuc.nab
./a.out

You should get a nuc.pdb file which contains the poly A-T double stranded DNA. Everything so far is available in more detail at http://ambermd.org/tutorials/basic/tutorial1/section2.htm

To get a single stranded poly-A DNA simply delete the lines corresponding to poly-T. To create stem-loop type configurations commonly found in RNA, do the following. Say you want the ssDNA

AAATAGAGCTAGTGAGGATTT

with the first 4 bp [AAAT] paired with the last 4 [TTTA or ATTT if you want to read from above]. For this I simply generate a double stranded DNA like so...

AAATAGAGCTAGTGAGGATTT
TTTATCTCGATCACTCCTAAA

Then delete the bases shown in red from the PDB file.

AAATAGAGCTAGTGAGGATTT
TTTATCTCGATCACTCCTAAA

Now you should have

AAATAGAGCTAGTGAGG
TTTA

This is the same as AAATAGAGCTAGTGAGGATTT but with the ends paired up. For the stem-loop conformation you should be able to use the parameter and topology files generated for the fully stretched single stranded DNA. You will have a long bond between  AAATAGAGCTAGTGAGG--and--ATTT. But this can be easily fixed with a quick minimization.