TI-based estimators¶
TI-based estimators such as TI take as input dHdl gradients for the calculation of free energy differences.
All TI-based estimators integrate these gradients with respect to \(\lambda\), differing only in how they numerically perform this integration.
As a usage example, we’ll use TI to calculate the free energy of solvation of benzene in water.
We’ll use the benzene-in-water dataset from alchemtest.gmx:
>>> from alchemtest.gmx import load_benzene
>>> bz = load_benzene().data
and parse the datafiles separately for each alchemical leg using alchemlyb.parsing.gmx.extract_dHdl() to obtain dHdl gradients:
>>> from alchemlyb.parsing.gmx import extract_dHdl
>>> import pandas as pd
>>> import alchemlyb
>>> dHdl_coul = alchemlyb.concat([extract_dHdl(xvg, T=300) for xvg in bz['Coulomb']])
>>> dHdl_vdw = alchemlyb.concat([extract_dHdl(xvg, T=300) for xvg in bz['VDW']])
We can now use the TI estimator to obtain the free energy differences between each \(\lambda\) window sampled.
The fit() method is used to perform the free energy estimate, given the gradient data:
>>> from alchemlyb.estimators import TI
>>> ti_coul = TI()
>>> ti_coul.fit(dHdl_coul)
TI(verbose=False)
# we could also just call the `fit` method
# directly, since it returns the `TI` object
>>> ti_vdw = TI().fit(dHdl_vdw)
The sum of the endpoint free energy differences will be the free energy of solvation for benzene in water.
The free energy differences (in units of \(k_B T\)) between each \(\lambda\) window can be accessed via the delta_f_ attribute:
>>> ti_coul.delta_f_
0.00 0.25 0.50 0.75 1.00
0.00 0.000000 1.620328 2.573337 3.022170 3.089027
0.25 -1.620328 0.000000 0.953009 1.401842 1.468699
0.50 -2.573337 -0.953009 0.000000 0.448832 0.515690
0.75 -3.022170 -1.401842 -0.448832 0.000000 0.066857
1.00 -3.089027 -1.468699 -0.515690 -0.066857 0.000000
So we can get the endpoint differences (free energy difference between \(\lambda = 0\) and \(\lambda = 1\)) of each with:
>>> ti_coul.delta_f_.loc[0.00, 1.00]
3.0890270218676896
>>> ti_vdw.delta_f_.loc[0.00, 1.00]
-3.0558175199846058
giving us a solvation free energy in units of \(k_B T\) for benzene of:
>>> ti_coul.delta_f_.loc[0.00, 1.00] + ti_vdw.delta_f_.loc[0.00, 1.00]
0.033209501883083803
In addition to the free energy differences, we also have access to the errors on these differences via the d_delta_f_ attribute:
>>> ti_coul.d_delta_f_
0.00 0.25 0.50 0.75 1.00
0.00 0.000000 0.009706 0.013058 0.015038 0.016362
0.25 0.009706 0.000000 0.008736 0.011486 0.013172
0.50 0.013058 0.008736 0.000000 0.007458 0.009858
0.75 0.015038 0.011486 0.007458 0.000000 0.006447
1.00 0.016362 0.013172 0.009858 0.006447 0.000000
Estimate free energy with gaussian quadrature¶
If the simulations are performed at certain gaussian quadrature points, TI_GQ can be used to estimate the free energy.
The usage of TI_GQ is similar to TI, but instead of the free energy differences between
\(\lambda\) windows, the values in the delta_f_ and d_delta_f_ tables are cumulative addition of estimation from one \(\lambda\) window to another.
To be consistent with TI and other estimators, the diagonal values are set to zeros and two end states at \(\lambda\)
0 and 1 are added, although the simulation may not be performed at \(\lambda\) 0 and 1. The value at \(\lambda\) 0 is set to zero and the value at \(\lambda\) 1
is the same as the previous gaussian quadrature point.
List of TI-based estimators¶
|
Thermodynamic integration (TI). |
|
Thermodynamic integration (TI) with gaussian quadrature estimation. |