.. _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` 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 `) we can again use :eq:`eqn. %s ` 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 ` or eq. :eq:`%s ` 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 ` 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 `, though the analysis currently requires external tools from the external `pymbar package `_. The :math:`\lambda`-dependence for the force-field contributions is described in detail in section sec. :ref:`feia`.