.. _fecalc:
Free energy calculations
------------------------
Slow-growth methods
~~~~~~~~~~~~~~~~~~~
Free energy calculations can be performed in |Gromacs| using a number of
methods, including “slow-growth.†An example problem might be
calculating the difference in free energy of binding of an inhibitor
**I** to an enzyme **E** and to a mutated enzyme
**E**\ :math:`^{\prime}`. It is not feasible with computer simulations
to perform a docking calculation for such a large complex, or even
releasing the inhibitor from the enzyme in a reasonable amount of
computer time with reasonable accuracy. However, if we consider the free
energy cycle in :numref:`Fig. %s A<fig-free1>` we can write:
.. math:: \Delta G_1 - \Delta G_2 = \Delta G_3 - \Delta G_4
:label: eqnddg
If we are interested in the left-hand term we can equally well compute
the right-hand term.
.. _fig-free1:
.. figure:: plots/free1.*
:width: 6.00000cm
Free energy cycles. **A:** to calculate :math:`\Delta G_{12}`, the free
energy difference between the binding of inhibitor **I** to enzymes
**E** respectively **E**\ :math:`^{\prime}`.
.. _fig-free2:
.. figure:: plots/free2.*
:width: 6.00000cm
Free energy cycles. **B:** to calculate
:math:`\Delta G_{12}`, the free energy difference for binding of
inhibitors **I** respectively **I**\ :math:`^{\prime}` to enzyme
**E**.
If we want to compute the difference in free energy of binding of two
inhibitors **I** and **I**\ :math:`^{\prime}` to an enzyme **E**
(:numref:`Fig. %s <fig-free2>`) we can again use
:eq:`eqn. %s <eqnddg>` to compute the desired property.
Free energy differences between two molecular species can be calculated
in |Gromacs| using the “slow-growth†method. Such free energy differences
between different molecular species are physically meaningless, but they
can be used to obtain meaningful quantities employing a thermodynamic
cycle. The method requires a simulation during which the Hamiltonian of
the system changes slowly from that describing one system (A) to that
describing the other system (B). The change must be so slow that the
system remains in equilibrium during the process; if that requirement is
fulfilled, the change is reversible and a slow-growth simulation from B
to A will yield the same results (but with a different sign) as a
slow-growth simulation from A to B. This is a useful check, but the user
should be aware of the danger that equality of forward and backward
growth results does not guarantee correctness of the results.
The required modification of the Hamiltonian :math:`H` is realized by
making :math:`H` a function of a *coupling parameter* :math:`\lambda:
H=H(p,q;\lambda)` in such a way that :math:`\lambda=0` describes system
A and :math:`\lambda=1` describes system B:
.. math:: H(p,q;0)=H{^{\mathrm{A}}}(p,q);~~~~ H(p,q;1)=H{^{\mathrm{B}}}(p,q).
:label: eqnddgHamiltonian
In |Gromacs|, the functional form of the :math:`\lambda`-dependence is
different for the various force-field contributions and is described in
section sec. :ref:`feia`.
The Helmholtz free energy :math:`A` is related to the partition function
:math:`Q` of an :math:`N,V,T` ensemble, which is assumed to be the
equilibrium ensemble generated by a MD simulation at constant volume and
temperature. The generally more useful Gibbs free energy :math:`G` is
related to the partition function :math:`\Delta` of an :math:`N,p,T`
ensemble, which is assumed to be the equilibrium ensemble generated by a
MD simulation at constant pressure and temperature:
.. math:: \begin{aligned}
A(\lambda) &=& -k_BT \ln Q \\
Q &=& c \int\!\!\int \exp[-\beta H(p,q;\lambda)]\,dp\,dq \\
G(\lambda) &=& -k_BT \ln \Delta \\
\Delta &=& c \int\!\!\int\!\!\int \exp[-\beta H(p,q;\lambda) -\beta
pV]\,dp\,dq\,dV \\
G &=& A + pV, \end{aligned}
:label: eqnddgGibs
where :math:`\beta = 1/(k_BT)` and :math:`c = (N! h^{3N})^{-1}`. These
integrals over phase space cannot be evaluated from a simulation, but it
is possible to evaluate the derivative with respect to :math:`\lambda`
as an ensemble average:
.. math:: \frac{dA}{d\lambda} = \frac{\int\!\!\int (\partial H/ \partial
\lambda) \exp[-\beta H(p,q;\lambda)]\,dp\,dq}{\int\!\!\int \exp[-\beta
H(p,q;\lambda)]\,dp\,dq} =
\left\langle \frac{\partial H}{\partial \lambda} \right\rangle_{NVT;\lambda},
:label: eqnddgensembleave
with a similar relation for :math:`dG/d\lambda` in the :math:`N,p,T`
ensemble. The difference in free energy between A and B can be found by
integrating the derivative over :math:`\lambda`:
.. math:: \begin{aligned}
A{^{\mathrm{B}}}(V,T)-A{^{\mathrm{A}}}(V,T) &=& \int_0^1 \left\langle \frac{\partial
H}{\partial \lambda} \right\rangle_{NVT;\lambda} \,d\lambda
\end{aligned}
:label: eqdelA
.. math:: \begin{aligned}
G{^{\mathrm{B}}}(p,T)-G{^{\mathrm{A}}}(p,T) &=& \int_0^1 \left\langle \frac{\partial
H}{\partial \lambda} \right\rangle_{NpT;\lambda} \,d\lambda.
\end{aligned}
:label: eqdelG
If one wishes to evaluate
:math:`G{^{\mathrm{B}}}(p,T)-G{^{\mathrm{A}}}(p,T)`, the natural choice
is a constant-pressure simulation. However, this quantity can also be
obtained from a slow-growth simulation at constant volume, starting with
system A at pressure :math:`p` and volume :math:`V` and ending with
system B at pressure :math:`p_B`, by applying the following small (but,
in principle, exact) correction:
.. math:: G{^{\mathrm{B}}}(p)-G{^{\mathrm{A}}}(p) =
A{^{\mathrm{B}}}(V)-A{^{\mathrm{A}}}(V) - \int_p^{p{^{\mathrm{B}}}}[V{^{\mathrm{B}}}(p')-V]\,dp'
:label: eqnddgpresscorr
Here we omitted the constant :math:`T` from the notation. This
correction is roughly equal to
:math:`-\frac{1}{2} (p{^{\mathrm{B}}}-p)\Delta V=(\Delta V)^2/(2
\kappa V)`, where :math:`\Delta V` is the volume change at :math:`p` and
:math:`\kappa` is the isothermal compressibility. This is usually small;
for example, the growth of a water molecule from nothing in a bath of
1000 water molecules at constant volume would produce an additional
pressure of as much as 22 bar, but a correction to the Helmholtz free
energy of just -1 kJ mol\ :math:`^{-1}`. In Cartesian coordinates, the
kinetic energy term in the Hamiltonian depends only on the momenta, and
can be separately integrated and, in fact, removed from the equations.
When masses do not change, there is no contribution from the kinetic
energy at all; otherwise the integrated contribution to the free energy
is :math:`-\frac{3}{2} k_BT \ln
(m{^{\mathrm{B}}}/m{^{\mathrm{A}}})`. **Note** that this is only true in
the absence of constraints.
Thermodynamic integration
~~~~~~~~~~~~~~~~~~~~~~~~~
|Gromacs| offers the possibility to integrate :eq:`eq. %s <eqdelA>` or eq.
:eq:`%s <eqdelG>` in one simulation over the full range from A to B. However, if
the change is large and insufficient sampling can be expected, the user
may prefer to determine the value of :math:`\langle
dG/d\lambda \rangle` accurately at a number of well-chosen intermediate
values of :math:`\lambda`. This can easily be done by setting the
stepsize ``delta_lambda`` to zero. Each simulation can be equilibrated
first, and a proper error estimate can be made for each value of
:math:`dG/d\lambda` from the fluctuation of :math:`\partial H/\partial
\lambda`. The total free energy change is then determined afterward by
an appropriate numerical integration procedure.
|Gromacs| now also supports the use of Bennett’s Acceptance Ratio
\ :ref:`58 <refBennett1976>` for calculating values of :math:`\Delta`\ G for transformations
from state A to state B using the program :ref:`gmx bar`. The same data can
also be used to calculate free energies using MBARÂ \ :ref:`59 <refShirts2008>`,
though the analysis currently requires external tools from the
external `pymbar package <https://SimTK.org/home/pymbar>`_.
The :math:`\lambda`-dependence for the force-field contributions is
described in detail in section sec. :ref:`feia`.