About the reproducing dF_state

[xiki-tempula]

I'm interested in porting the dF_state.pdf from the alchemical analysis. However, I'm not quite sure about the best interface of doing it.

My current interface is plot_dF_state(dhdl_data, orientation='portrait', nb=10), where orientation='portrait' and orientation='landscape' corresponds to the original dF_state.pdf and dF_state_long.pdf, nb is a parameter for the portrait cutoff.

For a typical ABFE calculation, we would have the bonded free energy, coulomb free energy and vdw free energy. To plot the TI, BAR and MBAR estimate, the dhdl_data would be:

from alchemtest.gmx import load_benzene
from alchemlyb.parsing.gmx import extract_u_nk, extract_dHdl
from alchemlyb.estimators import MBAR, TI, BAR

bz = load_benzene().data
u_nk_bond = pd.concat([extract_u_nk(xvg, T=300) for xvg in bz['Bonded']])
dHdl_bond = pd.concat([extract_dHdl(xvg, T=300) for xvg in bz['Bonded']])
u_nk_coul = pd.concat([extract_u_nk(xvg, T=300) for xvg in bz['Coulomb']])
dHdl_coul = pd.concat([extract_dHdl(xvg, T=300) for xvg in bz['Coulomb']])
u_nk_vdw = pd.concat([extract_u_nk(xvg, T=300) for xvg in bz['VDW']])
dHdl_vdw = pd.concat([extract_dHdl(xvg, T=300) for xvg in bz['VDW']])

ti_bond = TI().fit(dHdl_bond)
ti_coul = TI().fit(dHdl_coul)
ti_vdw = TI().fit(dHdl_vdw)
bar_bond = BAR().fit(u_nk_bond)
bar_coul = BAR().fit(u_nk_coul)
bar_vdw = BAR().fit(u_nk_vdw)
mbar_bond = MBAR().fit(u_nk_bond)
mbar_coul = MBAR().fit(u_nk_coul)
mbar_vdw = MBAR().fit(u_nk_vdw)

dhdl_data = [(ti_bond, ti_coul, ti_vdw),
                      (bar_bond, bar_coul, bar_vdw),
                      (mbar_bond, mbar_coul, mbar_vdw),]
fig = plot_dF_state(dhdl_data, orientation='portrait', nb=10)
fig.savefig('dF_state.pdf')

So I was also thinking of adding a wrapper function auto for this relatively tedious operation.

al = auto((bz['Bonded'], bz['Coulomb'], bz['Coulomb']), method=['TI', 'BAR', 'MBAR'])
al.subsampling(statistical_inefficiency, **kwargs) # I'm not quite sure of the subsampling part myself Issue #106 
al.run()
al.plot_dF_state(orientation='portrait', filename='dF_state.pdf')

I was wondering what is the thought of the community?

[orbeckst]

This looks like a complete analysis workflow. I can think of two useful approaches

1. Make it a tutorial and show it as a real-live example.
2. Add a new top-level module to alchemlyb where we collect common analysis workflows. In this case, alchemlyb.workflows.abfe. (I think this would be equivalent to MDAnalysis.analysis).

The first option is easy (and would also help with the many requests for better documentation). The second is harder because it's not always clear how general the workflows really are. On the other hand, if they are properly tested and documented then I would see that as a big plus.