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

[Miguel Caro]

I think the mail server messed up the latex formulas, hopefully this is
displayed better:

Hello,

I want to do thermodynamic integration directly from the ∂H/∂λ 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 and ∂H/∂λ_Cou , and so I am computing
⟨∂H/∂λ⟩=⟨∂H/∂λ_vdW ⟩∂λ_vdW /∂λ+⟨∂H/∂λ_Cou ⟩∂λ_Cou /∂λ "by hand".  I am
then using a spline to generate smoother curves in λ from the discrete
array (21 data points, one for each λ 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


On 2017-07-20 13:47, Miguel Caro wrote:
> 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