Repository files navigation

Code to do a binless WHAM (aka multistate Bennett acceptance ratio / MBAR). The code is flexible enough that, given a series of umbrella sampling windows along one reaction coordinate, can project to another reaction coordinate (assuming the reaction coordinates overlap sufficiently). It can also combine different umbrella sampling simulations performed with differen reaction coordnates (perhaps you were trying a few different reaction coordinates) to make a "super" PMF.

The input requires ametadata file, which essentially is the same as used in Alan Grossfield's WHAM code, except there is an additional argument indicating which column of data corresponds to the reaction coordinate the restraining force was applied to.

In other words, each row of themetadata file goes:

# (path to file)     (float)              (float)          (integer)trajectory_file_1 restraint_location_1 force_constant_1 reaction_coordinate_index_1trajectory_file_2 restraint_location_2 force_constant_2 reaction_coordinate_index_2...trajectory_file_n restraint_location_n force_constant_n reaction_coordinate_index_n

and eachtrajectory_file has the following columns:

Time trajectory_along_coordinate_1 trajectory_along_coordinate_2 ...trajectory_along_coordinate_n

It should go without saying that the trajectory file columns should be consistent among each other.

numpy andmatplotlib (if you want to plot)

Given ametadata and trajectory files available like above, here's a sample script:

from wham import WHAMwham = WHAM(T=300.0,metadata='metadata')wham.compute_free()wham.compute_pmf(hmin=-2.0,hmax=2.0,num_bins=40,pmf_crd=2)

which would compute the free energies from the umbrella sampling simulation defined along the reaction coordinate in yourmetadata file, then once those are computed it would build up the histogram from-2.0 to2.0 with 40 bins along the reaction coordinate in the second column of your trajectory data files.

You can also save and load free energies, to either avoid re-computing the unbiasing, or to compare with free energies obtained from another WHAM program.

wham = WHAM(T=300.0,metadata='metadata')wham.compute_free(save_free_energies='./myFreeEnergies.txt')wham.compute_pmf(hmin=-2.0,hmax=2.0,num_bins=40,pmf_crd=2,load_free_energies='./myFreeEnergies.txt')

would save the free energies computed incompute_free() to./myFreeEnergies.txt and then load them for use in plotting the PMF incompute_pmf(). Default save and load location is in the current folder:./free-energies.txt.