docs/reference-manual/free-energy-calculations.rst

   1 .. _fecalc:
   2
   3 Free energy calculations
   4 ------------------------
   5
   6 Slow-growth methods
   7 ~~~~~~~~~~~~~~~~~~~
   8
   9 Free energy calculations can be performed in |Gromacs| using a number of
  10 methods, including “slow-growth.” An example problem might be
  11 calculating the difference in free energy of binding of an inhibitor
  12 **I** to an enzyme **E** and to a mutated enzyme
  13 **E**\ :math:`^{\prime}`. It is not feasible with computer simulations
  14 to perform a docking calculation for such a large complex, or even
  15 releasing the inhibitor from the enzyme in a reasonable amount of
  16 computer time with reasonable accuracy. However, if we consider the free
  17 energy cycle in :numref:`Fig. %s A<fig-free1>` we can write:
  18
  19 .. math:: \Delta G_1 - \Delta G_2 =       \Delta G_3 - \Delta G_4
  20    :label: eqnddg
  21
  22 If we are interested in the left-hand term we can equally well compute
  23 the right-hand term.
  24
  25 .. _fig-free1:
  26
  27 .. figure:: plots/free1.*
  28             :width: 6.00000cm
  29
  30             Free energy cycles. **A:** to calculate :math:`\Delta G_{12}`, the free
  31             energy difference between the binding of inhibitor **I** to enzymes
  32             **E** respectively **E**\ :math:`^{\prime}`.
  33
  34 .. _fig-free2:
  35
  36 .. figure:: plots/free2.*
  37             :width: 6.00000cm
  38
  39             Free energy cycles. **B:** to calculate
  40             :math:`\Delta G_{12}`, the free energy difference for binding of
  41             inhibitors **I** respectively **I**\ :math:`^{\prime}` to enzyme
  42             **E**.
  43
  44 If we want to compute the difference in free energy of binding of two
  45 inhibitors **I** and **I**\ :math:`^{\prime}` to an enzyme **E**
  46 (:numref:`Fig. %s <fig-free2>`) we can again use
  47 :eq:`eqn. %s <eqnddg>` to compute the desired property.
  48
  49 Free energy differences between two molecular species can be calculated
  50 in |Gromacs| using the “slow-growth” method. Such free energy differences
  51 between different molecular species are physically meaningless, but they
  52 can be used to obtain meaningful quantities employing a thermodynamic
  53 cycle. The method requires a simulation during which the Hamiltonian of
  54 the system changes slowly from that describing one system (A) to that
  55 describing the other system (B). The change must be so slow that the
  56 system remains in equilibrium during the process; if that requirement is
  57 fulfilled, the change is reversible and a slow-growth simulation from B
  58 to A will yield the same results (but with a different sign) as a
  59 slow-growth simulation from A to B. This is a useful check, but the user
  60 should be aware of the danger that equality of forward and backward
  61 growth results does not guarantee correctness of the results.
  62
  63 The required modification of the Hamiltonian :math:`H` is realized by
  64 making :math:`H` a function of a *coupling parameter* :math:`\lambda:
  65 H=H(p,q;\lambda)` in such a way that :math:`\lambda=0` describes system
  66 A and :math:`\lambda=1` describes system B:
  67
  68 .. math:: H(p,q;0)=H{^{\mathrm{A}}}(p,q);~~~~ H(p,q;1)=H{^{\mathrm{B}}}(p,q).
  69
  70 In |Gromacs|, the functional form of the :math:`\lambda`-dependence is
  71 different for the various force-field contributions and is described in
  72 section sec. :ref:`feia`.
  73
  74 The Helmholtz free energy :math:`A` is related to the partition function
  75 :math:`Q` of an :math:`N,V,T` ensemble, which is assumed to be the
  76 equilibrium ensemble generated by a MD simulation at constant volume and
  77 temperature. The generally more useful Gibbs free energy :math:`G` is
  78 related to the partition function :math:`\Delta` of an :math:`N,p,T`
  79 ensemble, which is assumed to be the equilibrium ensemble generated by a
  80 MD simulation at constant pressure and temperature:
  81
  82 .. math::
  83
  84    \begin{aligned}
  85     A(\lambda) &=&  -k_BT \ln Q \\
  86     Q &=& c \int\!\!\int \exp[-\beta H(p,q;\lambda)]\,dp\,dq \\
  87     G(\lambda) &=&  -k_BT \ln \Delta \\
  88     \Delta &=& c \int\!\!\int\!\!\int \exp[-\beta H(p,q;\lambda) -\beta
  89    pV]\,dp\,dq\,dV \\
  90    G &=& A + pV, \end{aligned}
  91
  92 where :math:`\beta = 1/(k_BT)` and :math:`c = (N! h^{3N})^{-1}`. These
  93 integrals over phase space cannot be evaluated from a simulation, but it
  94 is possible to evaluate the derivative with respect to :math:`\lambda`
  95 as an ensemble average:
  96
  97 .. math::
  98
  99    \frac{dA}{d\lambda} =  \frac{\int\!\!\int (\partial H/ \partial
 100    \lambda) \exp[-\beta H(p,q;\lambda)]\,dp\,dq}{\int\!\!\int \exp[-\beta
 101    H(p,q;\lambda)]\,dp\,dq} =
 102    \left\langle \frac{\partial H}{\partial \lambda} \right\rangle_{NVT;\lambda},
 103
 104 with a similar relation for :math:`dG/d\lambda` in the :math:`N,p,T`
 105 ensemble. The difference in free energy between A and B can be found by
 106 integrating the derivative over :math:`\lambda`:
 107
 108 .. math::  \begin{aligned}
 109            A{^{\mathrm{B}}}(V,T)-A{^{\mathrm{A}}}(V,T) &=& \int_0^1 \left\langle \frac{\partial
 110            H}{\partial \lambda} \right\rangle_{NVT;\lambda} \,d\lambda
 111            \end{aligned}
 112            :label: eqdelA
 113
 114 .. math:: \begin{aligned}
 115           G{^{\mathrm{B}}}(p,T)-G{^{\mathrm{A}}}(p,T) &=& \int_0^1 \left\langle \frac{\partial
 116           H}{\partial \lambda} \right\rangle_{NpT;\lambda} \,d\lambda.
 117           \end{aligned}
 118           :label: eqdelG
 119
 120 If one wishes to evaluate
 121 :math:`G{^{\mathrm{B}}}(p,T)-G{^{\mathrm{A}}}(p,T)`, the natural choice
 122 is a constant-pressure simulation. However, this quantity can also be
 123 obtained from a slow-growth simulation at constant volume, starting with
 124 system A at pressure :math:`p` and volume :math:`V` and ending with
 125 system B at pressure :math:`p_B`, by applying the following small (but,
 126 in principle, exact) correction:
 127
 128 .. math::
 129
 130    G{^{\mathrm{B}}}(p)-G{^{\mathrm{A}}}(p) =
 131    A{^{\mathrm{B}}}(V)-A{^{\mathrm{A}}}(V) - \int_p^{p{^{\mathrm{B}}}}[V{^{\mathrm{B}}}(p')-V]\,dp'
 132
 133 Here we omitted the constant :math:`T` from the notation. This
 134 correction is roughly equal to
 135 :math:`-\frac{1}{2} (p{^{\mathrm{B}}}-p)\Delta V=(\Delta V)^2/(2
 136 \kappa V)`, where :math:`\Delta V` is the volume change at :math:`p` and
 137 :math:`\kappa` is the isothermal compressibility. This is usually small;
 138 for example, the growth of a water molecule from nothing in a bath of
 139 1000 water molecules at constant volume would produce an additional
 140 pressure of as much as 22 bar, but a correction to the Helmholtz free
 141 energy of just -1 kJ mol\ :math:`^{-1}`. In Cartesian coordinates, the
 142 kinetic energy term in the Hamiltonian depends only on the momenta, and
 143 can be separately integrated and, in fact, removed from the equations.
 144 When masses do not change, there is no contribution from the kinetic
 145 energy at all; otherwise the integrated contribution to the free energy
 146 is :math:`-\frac{3}{2} k_BT \ln
 147 (m{^{\mathrm{B}}}/m{^{\mathrm{A}}})`. **Note** that this is only true in
 148 the absence of constraints.
 149
 150 Thermodynamic integration
 151 ~~~~~~~~~~~~~~~~~~~~~~~~~
 152
 153 |Gromacs| offers the possibility to integrate :eq:`eq. %s <eqdelA>` or eq.
 154 :eq:`%s <eqdelG>` in one simulation over the full range from A to B. However, if
 155 the change is large and insufficient sampling can be expected, the user
 156 may prefer to determine the value of :math:`\langle
 157 dG/d\lambda \rangle` accurately at a number of well-chosen intermediate
 158 values of :math:`\lambda`. This can easily be done by setting the
 159 stepsize ``delta_lambda`` to zero. Each simulation can be equilibrated
 160 first, and a proper error estimate can be made for each value of
 161 :math:`dG/d\lambda` from the fluctuation of :math:`\partial H/\partial
 162 \lambda`. The total free energy change is then determined afterward by
 163 an appropriate numerical integration procedure.
 164
 165 |Gromacs| now also supports the use of Bennett’s Acceptance Ratio
 166 \ :ref:`58 <refBennett1976>` for calculating values of :math:`\Delta`\ G for transformations
 167 from state A to state B using the program :ref:`gmx bar`. The same data can
 168 also be used to calculate free energies using MBAR \ :ref:`59 <refShirts2008>`,
 169 though the analysis currently requires external tools from the
 170 external `pymbar package <https://SimTK.org/home/pymbar>`_.
 171
 172 The :math:`\lambda`-dependence for the force-field contributions is
 173 described in detail in section sec. :ref:`feia`.