For a result to be valid, we need to ensure that longer simulation time would not result in different results, i.e., that our results are converged. The
alchemlyb.convergence module provides functions to assess the convergence of free energy estimates or other quantities.
One way of determining the simulation end point is to compute and plot the
forward and backward convergence of the estimate using
>>> from alchemtest.gmx import load_benzene >>> from alchemlyb.parsing.gmx import extract_u_nk >>> from alchemlyb.visualisation import plot_convergence >>> from alchemlyb.convergence import forward_backward_convergence >>> bz = load_benzene().data >>> data_list = [extract_u_nk(xvg, T=300) for xvg in bz['Coulomb']] >>> df = forward_backward_convergence(data_list, 'mbar') >>> ax = plot_convergence(df) >>> ax.figure.savefig('dF_t.pdf')
Will give a plot looks like this
Fractional equilibration time¶
Another way of assessing whether the simulation has converged is to check the
energy files. In [Fan2021] (and [Fan2020]), \(R_c\) and
\(A_c\) are two criteria of checking the
fwdrev_cumavg_Rc() takes a decorrelated
pandas.Series as input and gives the metric
\(R_c\), which is 0 for fully-equilibrated
simulation and 1 for fully-unequilibrated simulation.
>>> from alchemtest.gmx import load_ABFE >>> from alchemlyb.parsing.gmx import extract_dHdl >>> from alchemlyb.preprocessing import decorrelate_dhdl, dhdl2series >>> from alchemlyb.convergence import fwdrev_cumavg_Rc >>> from alchemlyb.visualisation import plot_convergence >>> file = load_ABFE().data['ligand'] >>> dhdl = extract_dHdl(file, T=300) >>> decorrelated = decorrelate_dhdl(dhdl, remove_burnin=True) >>> R_c, running_average = fwdrev_cumavg_Rc(dhdl2series(decorrelated), tol=2) >>> print(R_c) 0.04 >>> ax = plot_convergence(running_average, final_error=2) >>> ax.set_ylabel("$\partial H/\partial\lambda$ (in kT)")
Will give a plot like this.
>>> from alchemlyb.convergence import A_c >>> dhdl_list =  >>> for file in load_ABFE().data['ligand']: >>> dhdl = extract_dHdl(file, T=300) >>> decorrelated = decorrelate_dhdl(dhdl, remove_burnin=True) >>> decorrelated = dhdl2series(decorrelated) >>> dhdl_list.append(decorrelated) >>> value = A_c(dhdl_list, tol=2) 0.7085
Convergence functions are available from
alchemlyb.convergence. Internally, they are imported from submodules, as documented below.
Functions for assessing convergence of free energy estimates and raw data.