alchemical-analysis.py not computing BAR/MBAR results

[mqbpksi2]

Dear Mobley lab,

I am trying to run the alchemical-analysis tool on a set of output files from Amber TI simulations. The tool successfully computes dG, but fails to compute BAR/MBAR results and output any plots.

My command line input:

alchemical_analysis -a AMBER -p l-prod-0.[1-9]-0 -q out -o . -r 5 -u kcal -s 0 -g -w

The output from alchemical_analysis:

Started on Thu Jan 26 11:12:41 2017
Command line was: /usr/local/bin/alchemical_analysis -a AMBER -p l-prod-0.[1-9]-0 -q out -o . -r 5 -u kcal -s 0 -g -w

Loading in data from ./l-prod-0.4-0.out... 1000 data points, 1 DV/DL averages
Loading in data from ./l-prod-0.6-0.out... 1000 data points, 1 DV/DL averages
Loading in data from ./l-prod-0.7-0.out... 1000 data points, 1 DV/DL averages
Loading in data from ./l-prod-0.5-0.out... 1000 data points, 1 DV/DL averages
Loading in data from ./l-prod-0.8-0.out... 1000 data points, 1 DV/DL averages
Loading in data from ./l-prod-0.9-0.out... 1000 data points, 1 DV/DL averages
Loading in data from ./l-prod-0.1-0.out... 1000 data points, 1 DV/DL averages
Loading in data from ./l-prod-0.3-0.out... 1000 data points, 1 DV/DL averages
Loading in data from ./l-prod-0.2-0.out... 1000 data points, 1 DV/DL averages

Sorting input data by starting time

Note: BAR/MBAR results are not computed.

Skipping first 0 steps (= 0.000000 ps)

Note: extrapolating missing gradient for lambda = 0.0: 58.481188
Note: extrapolating missing gradient for lambda = 1.0: -127.351642

The gradients (DV/DL) from the correlated samples (kcal/mol):

Lambda gradient

0.00000 58.481
0.10000 15.875
0.20000 -5.820
0.30000 -10.537
0.40000 -19.724
0.50000 -29.159
0.60000 -36.760
0.70000 -54.858
0.80000 -57.177
0.90000 -71.284
1.00000 -127.352

dG = -30.388

The correlated gradient (DV/DL) components from every_single step (kcal/mol):

Lambda BOND ANGLE DIHED 1-4 NB 1-4 EEL VDWAALS EELEC RESTRAINT

0.10000 2.512 0.854 -0.093 0.426 -16.137 3.537 24.232 0.000
0.20000 2.186 0.694 -0.095 0.420 -16.186 3.245 3.514 0.000
0.30000 1.932 0.580 -0.151 0.401 -16.272 5.021 -2.080 0.000
0.40000 1.835 0.351 -0.116 0.388 -16.119 3.152 -9.142 0.000
0.50000 1.758 0.577 -0.152 0.433 -16.059 -0.498 -15.740 0.000
0.60000 1.534 0.080 -0.400 0.313 -16.274 0.016 -21.656 0.000
0.70000 1.290 0.134 -0.126 0.372 -16.155 -9.769 -30.628 0.000
0.80000 1.085 -0.168 -0.227 0.340 -16.201 -0.396 -42.075 0.000
0.90000 0.898 -0.210 -0.134 0.353 -16.175 5.599 -61.267 0.000

dG = 1.606 0.283 -0.125 0.360 -16.266 -1.093 -14.374 0.000

WARNING: Found all 0 in input data.
Generating results.txt with all 0
Exiting...

The .out files can be found here: https://www.dropbox.com/sh/i4c44d0yrid1a9u/AABjyJjOMBykC9FWefZYCMNGa?dl=0

I have pymbar installed. Could it be that it's because these were generated by sander or because I'm using a lambda interval from 0.1 to 0.9? Any help and advice would be greatly appreciated.

Sincerely,

Stefan

@davidlmobley

[halx]

Yes, this is because sander cannot sample the end points. I could probably change the code such that BAR/MBAR free energies are calculated for the existing lambdas only (I assume pymbar would be fine with this). The user could then plot the free energy of each window vs. lambda and find a suitable extrapolation scheme (see discussion in a different thread of the issue tracker). The practical problem though is that sander provides MBAR energies only for equally spaced lambdas which limits what the minimum and maximum lambdas can be.

…

On 26 January 2017 at 11:29, mqbpksi2 ***@***.***> wrote: @davidlmobley <https://github.com/davidlmobley> — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#97 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AH3aWjbHlcPXBWmi4yP43HFQBHj4DA0Bks5rWIOAgaJpZM4LuiP6> .

[davidlmobley]

Yes, I'm with @halx - a key question in this case is what the expected behavior would be. If you're running a free energy calculation and you're lacking data at the endpoints, what OUGHT BAR and MBAR to give you?

In the other thread, we discussed extrapolation, but one view was that this is a bad idea since it will be difficult to know when it will break.

Historically, this has been why AMBER analysis of these types of simulations has been only with TI.

[mqbpksi2]

Many thanks for the feedback.

Best wishes,

Stef

[danmermelstein]

Hi,

I got the same error, but I was using the pmemd version of TI, which does compute endpoints. The analysis worked when I reran my simulations with MBAR on. @halx do you know if this code is only designed to work when MBAR is on? I'm just wondering because right now the GPU version of the code I have only supports TI, and for the time being it would be nice to be able to use this code for analysis/testing

Thanks!

Dan

P.S. Command line was: /net/home/dmermels/SOFTWARE/alchemical-analysis/alchemical_analysis/alchemical_analysis.py -a AMBER -d . -p [01].*/charge.gas_0-5ns -q out -m -mbar-bar -r 5 -o . -u kcal -s 0 -g

[davidlmobley]

Could you provide inputs, @danmermelstein ?

[danmermelstein]

Sure! Let me know if anything else is missing

Thanks!

mermelstein_alchemical_analysis_inputs.tar.gz

[halx]

Hi Dan,

I'm not sure what exactly the problem is. When I run your command with your data I get the output just fine. Only "issue" is that the tool also computes zero results for both BAR and MBAR. If you use "-m ti+ti-cubic" you can avoid that.

You have explicitly switch off MBAR reporting with ifmbar=0. The message you see is just a warning and can safely be ignored.

Cheers,
Hannes.

[danmermelstein]

Hi @halx,

Thanks for looking into this! What version of python and pymbar are you using?

Thanks!

Dan

python 2.7.13

[halx]

Python 2.7.9 and pymbar 3.0.0

…

On 20 April 2017 at 22:26, danmermelstein ***@***.***> wrote: Hi @halx <https://github.com/halx>, Thanks for looking into this! What version of python and pymbar are you using? Thanks! Dan python 2.7.13 — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#97 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AH3aWlR7jO0AhUGjHRCG_w16n-uv5hbKks5rx81tgaJpZM4LuiP6> .

[davidlmobley]

@danmermelstein - LMK if you think it will be helpful for me to try too.

[danmermelstein]

Hi Hannes and David,

I tried pymbar 3.0.0 and no luck, so I just disabled the system exit following the warning like Hannes said and I get the correct results. I'll post here if I figure out why that was happening, but may have just been an unnecessary exit.

Either way, problem solved, thanks for all of your help!

Dan

[mqbpksi2]

Hello all, Thanks for cc-ing me and letting me know about this, much appreciated. Best wishes, Stef

…

________________________________ From: danmermelstein [notifications@github.com] Sent: Friday, April 21, 2017 7:34 PM To: MobleyLab/alchemical-analysis Cc: Stefan Ivanov; Author Subject: Re: [MobleyLab/alchemical-analysis] alchemical-analysis.py not computing BAR/MBAR results (#97) Hi Hannes and David, I tried pymbar 3.0.0 and no luck, so I just disabled the system exit following the warning like Hannes said and I get the correct results. I'll post here if I figure out why that was happening, but may have just been an unnecessary exit. Either way, problem solved, thanks for all of your help! Dan — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub<#97 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AYMI_dYpKrCtUUnAqPpStgt56BCxyXGkks5ryPargaJpZM4LuiP6>.

[jchodera]

MBAR should absolutely be able to extrapolate to the endpoints even if they aren't sampled. Let me know if you folks need any help on the MBAR interface end of things.

Yes, this is because sander cannot sample the end points.

I believe this has was fixed with the introduction of softcore Lennard-Jones, which should be absolutely required for modern alchemical free energy calculations. Neglecting endpoint sampling in non-softcore free energy calculations has been shown to be highly problematic and has not been part of recommended best practices for many years now (nearly a decade).

[halx]

This is really an annoying implementation issue with sander. It does not allow you to set lambda < 0.005 or > 0.995. @danmermelstein can probably add more information. It is not an issue with pmemd though.

In this context, being able to estimate the non-sampled end-points would obviously be very interesting. Except, that I often get wrong results with MBAR (tested with pmemd result) that is the free energy can be very different from TI or BAR.

[jchodera]

Except, that I often get wrong results with MBAR (tested with pmemd result) that is the free energy can be very different from TI or BAR.

Can you elaborate on this? We find this universally indicates an issue with either (1) a problem in the thermodynamic self-consistency of your data, or (2) actual bias in the calculation that is impossible to detect via TI or BAR. TI allows you to miss near-singularities that hide the primary contributions to your integrals at the endpoints beyond or between the lambda values you evaluate, while one-sided BAR reduces to EXP, which has an enormous problem with bias being as large as the statistical error.