alchemlyb
alchemlyb copied to clipboard
alchemistry
documentation for processing AMBER18 TI runs
Hello All, Would it be possible to have some documentation on how to analyse my TI runs ( using AMBER18) using alchemlyb. I was trying to use Alchemical_analysis ( It is no longer maintained, I know, but I found tutorials on how to work with it). If some documentation is provided on how to use alchemlyb to "analyse my runs" and not "setup runs", it will be good. I have the data for 12 lambda windows ( gaussian quadrature) me for disappearing a small molecule . I am trying to perform absolute binding free energy calculation. My protein is ABL Kinase and my ligand is acetonitrile. Thanks everyone, any input will be deeply appreciated.
Hi @DDGmichigan , we are in great need of exactly the kind of documentation that you are mentioning. Unfortunately, nobody has committed the time to write them.
As a very poor substitute I can point you to the tests that show how to use the parts of the library in practice; they implement the workflow outlined in https://alchemlyb.readthedocs.io/en/latest/estimators-fep.html
First parse data files (from the AMBER BACE example) and create a u_nk dataframe https://github.com/alchemistry/alchemlyb/blob/c3d240b4bf0dc861249769594b514ac733545cc8/src/alchemlyb/tests/test_fep_estimators.py#L80-L85
Then run your estimator
est = alchemlyb.estimators.MBAR.fit(u_nk)
dG = est.delta_f_.iloc[0, -1]
sigma_dG = est.d_delta_f_.iloc[0, -1]
to get the free energy difference between first and last window together with the error estimate (which will only make sense if you subsampled the data before to only have independent datapoints).
Oops, you asked for TI and not MBAR. So then look at https://github.com/alchemistry/alchemlyb/blob/c3d240b4bf0dc861249769594b514ac733545cc8/src/alchemlyb/tests/test_ti_estimators.py#L79-L85 and use
est = alchemlyb.estimators.TI.fit(dHdl)
dG = est.delta_f_.iloc[0, -1]
sigma_dG = est.d_delta_f_.iloc[0, -1]
est = alchemlyb.estimators.MBAR.fit(u_nk)
Hello @orbeckst, I tried as what you suggested for MBAR estimation, but it failed. I could get the u_nk, but the estimation was failed due to: < est = alchemlyb.estimators.MBAR.fit(u_nk) < TypeError: fit() missing 1 required positional argument: 'u_nk'
I guess I still have to give a parameter to feed the data, any tips? Thanks a lot!
Not sure why your MBAR would fail with the command you show. Compare to the example https://alchemlyb.readthedocs.io/en/latest/visualisation.html#overlap-matrix-of-the-mbar — can you get the overlap matrix if you follow the steps starting at
>>> mbar_coul = MBAR()
Also, have you updated alchemlyb to the latest version?
Thanks @orbeckst, I tried, but I could not plot the overlap matrix due to:
pymbar.utils.ParameterError: Warning: Should have \sum_n W_nk = 1. Actual column sum for state 0 was 0.002835. 23 other columns have similar problems
I guess that the program did not read the Amber output correctly.
alchemlyb was newly installed, the version 0.5.1.
I put the files in the Dropbox, could you please have a check? Thanks a lot! https://www.dropbox.com/s/gzrl5wd11nfdm33/test_alchylb.zip?dl=0
Not sure why your MBAR would fail with the command you show. Compare to the example https://alchemlyb.readthedocs.io/en/latest/visualisation.html#overlap-matrix-of-the-mbar — can you get the overlap matrix if you follow the steps starting at
>>> mbar_coul = MBAR()Also, have you updated alchemlyb to the latest version?
I think this is what I would like to address in https://github.com/alchemistry/alchemlyb/issues/170, I think it is probably better to get this done than wait for a new pymbar release.
The new ABFE workflow gives a easy to use script of performing all the aforementioned analysis https://alchemlyb.readthedocs.io/en/latest/workflows/alchemlyb.workflows.ABFE.html