WHAM icon indicating copy to clipboard operation
WHAM copied to clipboard

Regarding PMF calculation with error bar (bootsrapping)

Open SoumyaCYZ opened this issue 1 year ago • 5 comments

Hope you are fine and doing well. I have performed an umbrella sampling simulation using the center of mass distance (dz) between membrane and peptide. I have pulled the peptide 4.4 nm to 0nm. Windows are nicely overlapped. I have tried to calculate the PMF using the following command--> "wham --min 0.0 --max 4.4 --bins 101 --temperature 310 --file metadata -o pmf.out --bt 100 --seed 1234". Can I fix the bins=101 (the total number of windows=41)? How should we fix the bin, tolerance, iterations, and bt parameters for WHAM? In some PMF calculations, I have seen bins=bt. It will be a great help if you tell me the protocol in detail. Apart from this what is the unit of free energy in the outputfile?

SoumyaCYZ avatar Dec 15 '24 18:12 SoumyaCYZ

Seems like your question is not really related to the tool, but to the method itself. I don't see a problem in setting bins greater than the number of windows. You just have to make sure there are still "enough" samples in each bin if the number of bins > number of windows.

See https://github.com/dnlbauer/WHAM/issues/16 for the free energy unit.

dnlbauer avatar Dec 17 '24 19:12 dnlbauer

Thanks, sir for your suggestion. I have checked all of my COLVAR files. Each window is nicely overlapped. I have a small doubt. When I used 10 or 20 bootstrapping steps (bt) the error was approx +/- 1.5 kJ/mol. But when I increased bt to 200 the error was +/- 0.3 kJ/mol. Is there any correct protocol to choose bt?

SoumyaCYZ avatar Dec 25 '24 18:12 SoumyaCYZ

In theory, more BT runs will improve the estimation. I would suggest to run it with an increasing number of BT steps and plot the approximated error with the steps. At some point, the error should not change with more steps anymore. That was the process I followed at least.

On 25 Dec 2024, 19:00, at 19:00, Soumya Mondal @.***> wrote:

Thanks, sir for your suggestion. I have checked all of my COLVAR files. Each window is nicely overlapped. I have a small doubt. When I used 10 or 20 bootstrapping steps (bt) the error was approx +/- 1.5 kJ/mol. But when I increased bt to 200 the error was +/- 0.3 kJ/mol. Is there any correct protocol to choose bt?

-- Reply to this email directly or view it on GitHub: https://github.com/dnlbauer/WHAM/issues/17#issuecomment-2561965171 You are receiving this because you commented.

Message ID: @.***>

dnlbauer avatar Dec 25 '24 18:12 dnlbauer

Yes sir I observed the same trend. After certain BT steps the errors are not changing. But in Grossfield's code bootstrapping broke because the very high free energies caused underflows on the probability, which in turn ruins the bootstrapping. I think the process involves a large energy barrier, bootstrap is not a good option. So in order to check the convergence we should plot the free energy profiles (for e.g. 0-10, 0-20,0-30, and 0-40 ns) and compare (block averaging).

SoumyaCYZ avatar Dec 28 '24 11:12 SoumyaCYZ

The two methods are serving different purposes. Bootstrapping is for quantifying the error, not for checking convergence since it takes random samples while convergence is happing along the time coordinate.

Yes, to check convergence, plotting blocks with increasing time intervals is a valid option.

The bootstrapping approach can be used afterwards to quantify the uncertainty of the calculated energy from the simulation. However, for bootstrapping to give meaningful results, you should first check for convergence and apply it to the converged trajectory.

On 28 Dec 2024, 12:48, at 12:48, Soumya Mondal @.***> wrote:

Yes sir I observed the same trend. After certain BT steps the errors are not changing. But in Grossfield's code bootstrapping broke because the very high free energies caused underflows on the probability, which in turn ruins the bootstrapping. I think the process involves a large energy barrier, bootstrap is not a good option. So in order to check the convergence we should plot the free energy profiles (for e.g. 0-10, 0-20,0-30, and 0-40 ns) and compare (block averaging).

-- Reply to this email directly or view it on GitHub: https://github.com/dnlbauer/WHAM/issues/17#issuecomment-2564305194 You are receiving this because you commented.

Message ID: @.***>

dnlbauer avatar Dec 28 '24 12:12 dnlbauer