Weighted Histogram Analysis Method

Download: www.seas.upenn.edu/~biophys/dissemination/wham2d.f
Download: www.seas.upenn.edu/~biophys/dissemination/wham3d.f

The code to perform the 2D-WHAM is the fortran program, wham2d.f.  There are some files that need to be prepared before hand before the code can be run.  The program was set up to do the wham analysis upon a CHARMM simulation, I altered it so should work for a NAMD simulation with VMD analysis.

The MD simulation needs to have two sets of harmonic constraints upon the simulation.  With a set of constraints is along each reaction coordinate.  For the conformational RMSD constraints using an initial and final structure, one set of constraints to the initial reference structure and one set of constraints to the final reference structure.  However the constraints are applied, the RMSD constraints need to be stepped through multiple simulations to sample the entire pathway: i.e. sim1 to sim15.

Following the simulation the RMSD needs to taken with regards to the initial and final structures at each simulation and stored in separate files, sim1.dat1 contains the RMSD to the initial structure for the first simulation, sim1.dat2 contains RMSD to the final structure for the first simulation and so on to sim15.dat1 and sim15.dat2.  wham2d.f is set to read in dat files with only RMSD values contained with in them.  If you want to change that see the comments in wham2d.f

A file called sampling needs to be created with the filenames of the dat1 and dat2 files, the value harmonic constraint and the length of the constraint to the initial structure, the similar values for the final structure and the number of data points in each simulation run.  For the example 
	sim01	20.0	0.0	20.0	7.0	2000
	sim02	20.0	0.5	20.0	6.5	2000	
	sim03	20.0	1.0	20.0	6.0	2000
	…	
	sim15   20.0 	7.0	20.0	0.0	2000

One note is that since the program is fortran, it is very specific in the number of characters in each field.  So in the sampling file, you must keep the filename the same number of characters.  

The final set of alterations are within wham2d.f to set the number of characters in each dat1 and dat2 file, as well as the max and min along each dimension.  

The number of characters in each file is set in the declaration 
       character(NUM) :: filename, chdummy
       character(NUM+5) :: filename1, filename2, filename3,
			     filename4, filename5
where NUM is the number of characters of the names in the dat1 and dat2 files.  The following line must have 5 more characters than the previous line to account for the full filename.

This is done in the top of the program with the variables dim1_min, dim1_max, dim2_min, dim2_max.  These variables need to be set to contain the maximum and minimum values of the RMSD from the dat1 and dat2 files respectively.  If these are not set correctly, the algorithm will not converge correctly.  

To compile wham2d.f you need to use a fortran 90 compiler.  When the program is executed it will ask for the number of windows and the number of bins.  The number of windows are the number of simulations while the number of bins will provide the energy resolution.  After running the program there will be three output files: dfe2.dat, axis2d.dat and scatter.dat.

Scatter.dat holds the value of the 2D-rmsd plots and plotting scatter.dat in the tool of you choice is an easy way to determine if there are enough sampling points along the reaction coordinate.  If there are not enough data points rerun simulations to fill in the data points and tack all this information at the end of the sampling file, it doesn’t need to be in order.

Dfe2.dat holds the energetic map calculated by the WHAM algorithm, while axis2d.dat holds the axis marking for plotting purposes.