Assessing convergence¶
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.
Time Convergence¶
One way of determining the simulation end point is to compute and plot the
forward and backward convergence of the estimate using
forward_backward_convergence()
and
plot_convergence()
[Klimovich2015].
>>> 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
convergence. 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'][0]
>>> 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.
The A_c()
on the other hand, takes in a list of
decorrelated pandas.Series
and gives a metric of how converged the set
is, where 0 fully-unequilibrated and 1.0 is fully-equilibrated.
>>> 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¶
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. |