Visualisation of the results¶
It is quite often that the user want to visualise the results to gain confidence on the computed free energy. alchemlyb provides various visualisation tools to help user to judge the estimate.
Plotting Functions¶
|
Plot the MBAR overlap matrix. |
|
Plot the dhdl of TI. |
|
Plot the dhdl of TI. |
|
Plot the forward and backward convergence. |
Overlap Matrix of the MBAR¶
The accuracy of the MBAR estimator depends on
the overlap between different lambda states. The overlap matrix from the
MBAR estimator could be plotted using
plot_mbar_overlap_matrix() to check
the degree of overlap. It is recommended that there should be at least
0.03 [Klimovich2015] overlap between neighboring states.
>>> import pandas as pd
>>> from alchemtest.gmx import load_benzene
>>> from alchemlyb.parsing.gmx import extract_u_nk
>>> from alchemlyb.estimators import MBAR
>>> bz = load_benzene().data
>>> u_nk_coul = alchemlyb.concat([extract_u_nk(xvg, T=300) for xvg in bz['Coulomb']])
>>> mbar_coul = MBAR()
>>> mbar_coul.fit(u_nk_coul)
>>> from alchemlyb.visualisation import plot_mbar_overlap_matrix
>>> ax = plot_mbar_overlap_matrix(mbar_coul.overlap_matrix)
>>> ax.figure.savefig('O_MBAR.pdf', bbox_inches='tight', pad_inches=0.0)
Will give a plot looks like this
Overlap between the distributions of potential energy differences is essential for accurate free energy calculations and can be quantified by computing the overlap matrix 𝐎. Its elements 𝑂𝑖𝑗 are the probabilities of observing a sample from state i (𝑖 th row) in state j (𝑗 th column).¶
dhdl Plot of the TI¶
In order for the TI estimator to work reliably,
the change in the dhdl between lambda state 0 and lambda state 1 should be
adequately sampled. The function plot_ti_dhdl()
can be used to assess the change of the dhdl across the lambda states.
More than one TI estimators can be plotted
together as well.
>>> import pandas as pd
>>> from alchemtest.gmx import load_benzene
>>> from alchemlyb.parsing.gmx import extract_dHdl
>>> from alchemlyb.estimators import TI
>>> bz = load_benzene().data
>>> dHdl_coul = alchemlyb.concat([extract_dHdl(xvg, T=300) for xvg in bz['Coulomb']])
>>> ti_coul = TI().fit(dHdl_coul)
>>> dHdl_vdw = alchemlyb.concat([extract_dHdl(xvg, T=300) for xvg in bz['VDW']])
>>> ti_vdw = TI().fit(dHdl_vdw)
>>> from alchemlyb.visualisation import plot_ti_dhdl
>>> ax = plot_ti_dhdl([ti_coul, ti_vdw], labels=['Coul', 'VDW'], colors=['r', 'g'])
>>> ax.figure.savefig('dhdl_TI.pdf')
Will give a plot looks like this
A plot of ⟨∂𝑈/∂𝜆⟩ versus 𝜆 for thermodynamic integration, with filled areas indicating free energy estimates from the trapezoid rule. Different 𝛥𝐺 components are shown in distinct colors: in red is the electrostatic 𝛥𝐺 component (𝜆 indices 0–4), while in green is the van der Waals 𝛥𝐺 component (𝜆 indices 5–19). Color intensity alternates with increasing 𝜆 index.¶
dF States Plots between Different estimators¶
Another way of assessing the quality of free energy estimate would be comparing
the free energy difference between adjacent lambda states (dF) using different
estimators [Klimovich2015]. The function plot_dF_state() can
be used, for example, to compare the dF of both Coulombic and VDW
transformations using TI,
BAR and MBAR
estimators.
>>> from alchemtest.gmx import load_benzene
>>> from alchemlyb.parsing.gmx import extract_u_nk, extract_dHdl
>>> from alchemlyb.estimators import MBAR, TI, BAR
>>> import matplotlib.pyplot as plt
>>> import pandas as pd
>>> from alchemlyb.visualisation.dF_state import plot_dF_state
>>> bz = load_benzene().data
>>> u_nk_coul = alchemlyb.concat([extract_u_nk(xvg, T=300) for xvg in bz['Coulomb']])
>>> dHdl_coul = alchemlyb.concat([extract_dHdl(xvg, T=300) for xvg in bz['Coulomb']])
>>> u_nk_vdw = alchemlyb.concat([extract_u_nk(xvg, T=300) for xvg in bz['VDW']])
>>> dHdl_vdw = alchemlyb.concat([extract_dHdl(xvg, T=300) for xvg in bz['VDW']])
>>> ti_coul = TI().fit(dHdl_coul)
>>> ti_vdw = TI().fit(dHdl_vdw)
>>> bar_coul = BAR().fit(u_nk_coul)
>>> bar_vdw = BAR().fit(u_nk_vdw)
>>> mbar_coul = MBAR().fit(u_nk_coul)
>>> mbar_vdw = MBAR().fit(u_nk_vdw)
>>> estimators = [(ti_coul, ti_vdw),
(bar_coul, bar_vdw),
(mbar_coul, mbar_vdw),]
>>> fig = plot_dF_state(estimators, orientation='portrait')
>>> fig.savefig('dF_state.pdf', bbox_inches='tight')
Will give a plot looks like this
A bar plot of the free energy differences evaluated between pairs of adjacent states via several methods, with corresponding error estimates for each method.¶
Forward and Backward Convergence¶
One way of determining the simulation end point is to plot the forward and
backward convergence of the estimate using
plot_convergence().
Note that this is just a plotting function to plot [Klimovich2015] style convergence plot. The user need to provide the forward and backward data list and the corresponding error.
>>> import pandas as pd
>>> from alchemtest.gmx import load_benzene
>>> from alchemlyb.parsing.gmx import extract_u_nk
>>> from alchemlyb.estimators import MBAR
>>> bz = load_benzene().data
>>> data_list = [extract_u_nk(xvg, T=300) for xvg in bz['Coulomb']]
>>> forward = []
>>> forward_error = []
>>> backward = []
>>> backward_error = []
>>> num_points = 10
>>> for i in range(1, num_points+1):
>>> # Do the forward
>>> slice = int(len(data_list[0])/num_points*i)
>>> u_nk_coul = alchemlyb.concat([data[:slice] for data in data_list])
>>> estimate = MBAR().fit(u_nk_coul)
>>> forward.append(estimate.delta_f_.iloc[0,-1])
>>> forward_error.append(estimate.d_delta_f_.iloc[0,-1])
>>> # Do the backward
>>> u_nk_coul = alchemlyb.concat([data[-slice:] for data in data_list])
>>> estimate = MBAR().fit(u_nk_coul)
>>> backward.append(estimate.delta_f_.iloc[0,-1])
>>> backward_error.append(estimate.d_delta_f_.iloc[0,-1])
>>> from alchemlyb.visualisation import plot_convergence
>>> ax = plot_convergence(forward, forward_error, backward, backward_error)
>>> ax.figure.savefig('dF_t.pdf')
Will give a plot looks like this
A convergence plot of showing that the forward and backward has converged fully.¶