docs/reference-manual/algorithms/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           :label: eqnddgHamiltonian
  70
  71 In |Gromacs|, the functional form of the :math:`\lambda`-dependence is
  72 different for the various force-field contributions and is described in
  73 section sec. :ref:`feia`.
  74
  75 The Helmholtz free energy :math:`A` is related to the partition function
  76 :math:`Q` of an :math:`N,V,T` ensemble, which is assumed to be the
  77 equilibrium ensemble generated by a MD simulation at constant volume and
  78 temperature. The generally more useful Gibbs free energy :math:`G` is
  79 related to the partition function :math:`\Delta` of an :math:`N,p,T`
  80 ensemble, which is assumed to be the equilibrium ensemble generated by a
  81 MD simulation at constant pressure and temperature:
  82
  83 .. math:: \begin{aligned}
  84            A(\lambda) &=&  -k_BT \ln Q \\
  85            Q &=& c \int\!\!\int \exp[-\beta H(p,q;\lambda)]\,dp\,dq \\
  86            G(\lambda) &=&  -k_BT \ln \Delta \\
  87            \Delta &=& c \int\!\!\int\!\!\int \exp[-\beta H(p,q;\lambda) -\beta
  88           pV]\,dp\,dq\,dV \\
  89           G &=& A + pV, \end{aligned}
  90           :label: eqnddgGibs
  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:: \frac{dA}{d\lambda} =  \frac{\int\!\!\int (\partial H/ \partial
  98           \lambda) \exp[-\beta H(p,q;\lambda)]\,dp\,dq}{\int\!\!\int \exp[-\beta
  99           H(p,q;\lambda)]\,dp\,dq} =
 100           \left\langle \frac{\partial H}{\partial \lambda} \right\rangle_{NVT;\lambda},
 101           :label: eqnddgensembleave
 102
 103 with a similar relation for :math:`dG/d\lambda` in the :math:`N,p,T`
 104 ensemble. The difference in free energy between A and B can be found by
 105 integrating the derivative over :math:`\lambda`:
 106
 107 .. math::  \begin{aligned}
 108            A{^{\mathrm{B}}}(V,T)-A{^{\mathrm{A}}}(V,T) &=& \int_0^1 \left\langle \frac{\partial
 109            H}{\partial \lambda} \right\rangle_{NVT;\lambda} \,d\lambda
 110            \end{aligned}
 111            :label: eqdelA
 112
 113 .. math:: \begin{aligned}
 114           G{^{\mathrm{B}}}(p,T)-G{^{\mathrm{A}}}(p,T) &=& \int_0^1 \left\langle \frac{\partial
 115           H}{\partial \lambda} \right\rangle_{NpT;\lambda} \,d\lambda.
 116           \end{aligned}
 117           :label: eqdelG
 118
 119 If one wishes to evaluate
 120 :math:`G{^{\mathrm{B}}}(p,T)-G{^{\mathrm{A}}}(p,T)`, the natural choice
 121 is a constant-pressure simulation. However, this quantity can also be
 122 obtained from a slow-growth simulation at constant volume, starting with
 123 system A at pressure :math:`p` and volume :math:`V` and ending with
 124 system B at pressure :math:`p_B`, by applying the following small (but,
 125 in principle, exact) correction:
 126
 127 .. math:: G{^{\mathrm{B}}}(p)-G{^{\mathrm{A}}}(p) =
 128           A{^{\mathrm{B}}}(V)-A{^{\mathrm{A}}}(V) - \int_p^{p{^{\mathrm{B}}}}[V{^{\mathrm{B}}}(p')-V]\,dp'
 129           :label: eqnddgpresscorr
 130
 131 Here we omitted the constant :math:`T` from the notation. This
 132 correction is roughly equal to
 133 :math:`-\frac{1}{2} (p{^{\mathrm{B}}}-p)\Delta V=(\Delta V)^2/(2
 134 \kappa V)`, where :math:`\Delta V` is the volume change at :math:`p` and
 135 :math:`\kappa` is the isothermal compressibility. This is usually small;
 136 for example, the growth of a water molecule from nothing in a bath of
 137 1000 water molecules at constant volume would produce an additional
 138 pressure of as much as 22 bar, but a correction to the Helmholtz free
 139 energy of just -1 kJ mol\ :math:`^{-1}`. In Cartesian coordinates, the
 140 kinetic energy term in the Hamiltonian depends only on the momenta, and
 141 can be separately integrated and, in fact, removed from the equations.
 142 When masses do not change, there is no contribution from the kinetic
 143 energy at all; otherwise the integrated contribution to the free energy
 144 is :math:`-\frac{3}{2} k_BT \ln
 145 (m{^{\mathrm{B}}}/m{^{\mathrm{A}}})`. **Note** that this is only true in
 146 the absence of constraints.
 147
 148 Thermodynamic integration
 149 ~~~~~~~~~~~~~~~~~~~~~~~~~
 150
 151 |Gromacs| offers the possibility to integrate :eq:`eq. %s <eqdelA>` or eq.
 152 :eq:`%s <eqdelG>` in one simulation over the full range from A to B. However, if
 153 the change is large and insufficient sampling can be expected, the user
 154 may prefer to determine the value of :math:`\langle
 155 dG/d\lambda \rangle` accurately at a number of well-chosen intermediate
 156 values of :math:`\lambda`. This can easily be done by setting the
 157 stepsize ``delta_lambda`` to zero. Each simulation can be equilibrated
 158 first, and a proper error estimate can be made for each value of
 159 :math:`dG/d\lambda` from the fluctuation of :math:`\partial H/\partial
 160 \lambda`. The total free energy change is then determined afterward by
 161 an appropriate numerical integration procedure.
 162
 163 |Gromacs| now also supports the use of Bennett’s Acceptance Ratio
 164 \ :ref:`58 <refBennett1976>` for calculating values of :math:`\Delta`\ G for transformations
 165 from state A to state B using the program :ref:`gmx bar`. The same data can
 166 also be used to calculate free energies using MBAR \ :ref:`59 <refShirts2008>`,
 167 though the analysis currently requires external tools from the
 168 external `pymbar package <https://SimTK.org/home/pymbar>`_.
 169
 170 The :math:`\lambda`-dependence for the force-field contributions is
 171 described in detail in section sec. :ref:`feia`.