Use MBAR error as uncertainty with a single protocol repeat in RBFE

[frannerin]

When only a single repeat of the RBFE ProtocolUnit is used in a project (which is the approach in e.g., FEP+), the MBAR error estimate should be used as the uncertainty of the free energy estimate. If not, the uncertainty is 0 (applying np.std() to a list of one number returns 0).

Checklist

• Added a news entry

Developers certificate of origin

• I certify that this contribution is covered by the MIT License here and the Developer Certificate of Origin at https://developercertificate.org/.

[pep8speaks]

Hello @frannerin! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

• In the file openfecli/commands/gather.py:

Line 182:45: E261 at least two spaces before inline comment
Line 185:48: E261 at least two spaces before inline comment
Line 195:50: E261 at least two spaces before inline comment
Line 196:62: W291 trailing whitespace
Line 222:48: E261 at least two spaces before inline comment
Line 268:48: E261 at least two spaces before inline comment
Line 340:5: E266 too many leading '#' for block comment
Line 351:1: W293 blank line contains whitespace
Line 363:9: E266 too many leading '#' for block comment
Line 366:1: W293 blank line contains whitespace
Line 380:80: E501 line too long (86 > 79 characters)
Line 381:1: W293 blank line contains whitespace
Line 391:46: E261 at least two spaces before inline comment

Comment last updated at 2024-07-05 09:34:35 UTC

[IAlibay]

@frannerin - I think this requires more discussion. I don't believe this is the right way to do this, as you end up with a case where the return can mean two vastly different things statistically.

The programmatic approach to getting the mbar errors is via get_individual_estimates, although I recognise that it isn't the most intuitive thing nor does it play very well with the CLI. What is your use case?

[frannerin]

I'm running the rbfe tutorial and wanted to do everything, up to using the "gather" cli command and plotting with cinnabar, but making only 1 protocol repeat and getting !=0 errors in the plots. After, I will try my own system and ligand set but that's a different beast

[IAlibay]

Thanks, I think the solution here might be to fix the CLI gather command to fetch mbar errors instead if there is insufficient data to do a standard deviation (rather than changing the behaviour of the results object itself).

[frannerin]

So for gather the column should still be named "uncertainty (kcal/mol)" even if MBAR error est. is used (and e.g. print a warning)? I'm thinking in case there are transformations with different number of protocol repeats that are gathered together, for example if gather is run before all calculations finished, or some other weird extreme cases (or put a check/assert so that all transformations must have the same number of protocolunits)

[IAlibay]

So for gather the column should still be named "uncertainty (kcal/mol)" even if MBAR error est. is used (and e.g. print a warning)? I'm thinking in case there are transformations with different number of protocol repeats that are gathered together, for example if gather is run before all calculations finished, or some other weird extreme cases (or put a check/assert so that all transformations must have the same number of protocolunits)

That's a good question - and where it gets really messy.

This definitely would require a longer chat (would you be interested in joining a call about it at some point?).

A take (and not guaranteed to be the right one - definitely needs more thinking through!):

1. Gather should populate a "stdev (kcal/mol)" and "mbar error (kcal/mol)" column, where the former is just the standard deviation, and the latter is the propagated MBAR error (I think that would be something like sqrt(err1^2 + err2^2.. + errN^2)/ N)?
2. If all the samples have a standard deviation, then use the standard deviation when passing things to cinnabar.
3. If not, raise a UserWarning and use the MBAR errors.

Note: we're in a 1.x release, so renaming "uncertainty" might be a tad bit problematic - that might need some consideration (nothing stops us from starting a 2.0 branch though!).

[IAlibay]

@frannerin may I please suggest that we have that chat before you spend a lot of time working on this? I would hate to see you spend a lot of efforts on a solution just for us to turn around and claim we want something else.

[frannerin]

Hi @IAlibay yes I'm open for a chat, but I'm no expert I have to warn you, and I don't think I could contribute much. The last pushed commits where honestly for me and forgot they would show up here hehe 🐋 What you propose in your previous comment seems reasonable to me