[gmx-users] Computing free energy differences from the dH/dl values directly

[Miguel Caro]

Hello,

I want to do thermodynamic integration directly from the ∂H/∂λ\partial H
/ \partial \lambda values in the md.xvg files, rather than using gmx
bar, because I want to do several manipulations (eventually). On the
tutorial on solvated methane how to use gmx bar is explained, however
this is a black box tool so I would like to understand a bit better what
it's doing and how it handles the raw derivatives.

I have a system (one methanol solvated in water) where I decouple
Coulomb interactions before decoupling vdW interactions. The md.xvg
files give me ∂H/∂λvdW\partial H / \partial \lambda_{vdW} and
∂H/∂λcou\partial H / \partial \lambda_{cou}, and so I am
computing ⟨∂H/∂λ⟩=⟨∂H/∂λvdW⟩∂λvdW/∂λ+⟨∂H/∂λc⟩∂λc/∂λ\langle \partial H /
\partial \lambda \rangle = \langle \partial H / \partial \lambda_{vdW}
\rangle \partial \lambda_{vdW} / \partial \lambda + \langle \partial H /
\partial \lambda_{c} \rangle \partial \lambda_{c} / \partial \lambda "by
hand".  I am then using a spline to generate smoother curves in λ\lambda
from the discrete array (21 data points, one for each λ\lambda value) of
the expectation values. I then integrate this smooth curve, although I
guess I could also use a quadrature rule applied to the original data
before smoothing. The final free energy difference I obtain like this is
about 5% larger than the value given by gmx bar.

I would like to know if what I am doing makes sense and if the
difference with gmx bar is to be expected (for instance because of how
the integral is performed or some other fundamental difference).

Many thanks,

Miguel

-- 

*Dr. Miguel Caro*
/Postdoctoral researcher/
Department of Electrical Engineering and Automation,
and COMP Centre of Excellence in Computational Nanoscience
Aalto University, Finland
Personal email: *mcaroba at gmail.com*
Work: *miguel.caro at aalto.fi*
Website: http://mcaroba.dyndns.org