docs/reference-manual/special/free-energy-implementation.rst

   1 .. _dgimplement:
   2
   3 Free energy implementation
   4 --------------------------
   5
   6 For free energy calculations, there are two things that must be
   7 specified; the end states, and the pathway connecting the end states.
   8 The end states can be specified in two ways. The most straightforward is
   9 through the specification of end states in the topology file. Most
  10 potential forms support both an :math:`A` state and a :math:`B` state.
  11 Whenever both states are specified, the :math:`A` state corresponds
  12 to the initial free energy state, and the :math:`B` state corresponds to
  13 the final state.
  14
  15 In some cases, the end state can also be defined in some cases without
  16 altering the topology, solely through the :ref:`mdp` file,
  17 through the use of the
  18 ``couple-moltype``,
  19 ``couple-lambda0``,
  20 ``couple-lambda1``, and ``couple-intramol`` :ref:`mdp`
  21 keywords. Any molecule type selected in ``couple-moltype``
  22 will automatically have a :math:`B` state implicitly constructed (and
  23 the :math:`A` state redefined) according to the
  24 ``couple-lambda`` keywords. ``couple-lambda0``
  25 and ``couple-lambda1`` define the non-bonded parameters that
  26 are present in the :math:`A` state (``couple-lambda0``) and
  27 the :math:`B` state (``couple-lambda1``). The choices are
  28 ``q``,
  29 ``vdw``, and ``vdw-q``; these indicate the Coulombic, van der Waals, or
  30 both parameters that are turned on in the respective state.
  31
  32 Once the end states are defined, then the path between the end states
  33 has to be defined. This path is defined solely in the .mdp file.
  34 Starting in 4.6, :math:`\lambda` is a vector of components, with
  35 Coulombic, van der Waals, bonded, restraint, and mass components all
  36 able to be adjusted independently. This makes it possible to turn off
  37 the Coulombic term linearly, and then the van der Waals using soft core,
  38 all in the same simulation. This is especially useful for replica
  39 exchange or expanded ensemble simulations, where it is important to
  40 sample all the way from interacting to non-interacting states in the
  41 same simulation to improve sampling.
  42
  43 ``fep-lambdas`` is the default array of :math:`\lambda`
  44 values ranging from 0 to 1. All of the other lambda arrays use the
  45 values in this array if they are not specified. The previous behavior,
  46 where the pathway is controlled by a single :math:`\lambda` variable,
  47 can be preserved by using only ``fep-lambdas`` to define the
  48 pathway.
  49
  50
  51 .. _figlambdaval:
  52
  53 .. figure:: plots/lambda-values.*
  54    :width: 12.00000cm
  55
  56    Separate :math:`\lambda` values for Coulomb, van-der-Waals and restraint interactions.
  57
  58 :numref:`Fig. %s <figlambdaval>` shows an example of different lambda arrays.
  59 There, first the Coulombic terms are reduced, then
  60 the van der Waals terms, changing bonded at the same time rate as the
  61 van der Waals, but changing the restraints throughout the first
  62 two-thirds of the simulation. The corresponding :math:`\lambda`
  63 vector is given here:
  64
  65 ::
  66
  67     coul-lambdas           = 0.0 0.2 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0
  68     vdw-lambdas            = 0.0 0.0 0.0 0.0 0.4 0.5 0.6 0.7 0.8 1.0
  69     bonded-lambdas         = 0.0 0.0 0.0 0.0 0.4 0.5 0.6 0.7 0.8 1.0
  70     restraint-lambdas      = 0.0 0.0 0.1 0.2 0.3 0.5 0.7 1.0 1.0 1.0
  71
  72 This is also equivalent to:
  73
  74 ::
  75
  76     fep-lambdas            = 0.0 0.0 0.0 0.0 0.4 0.5 0.6 0.7 0.8 1.0
  77     coul-lambdas           = 0.0 0.2 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0
  78     restraint-lambdas      = 0.0 0.0 0.1 0.2 0.3 0.5 0.7 1.0 1.0 1.0
  79
  80 The ``fep-lambda`` array, in this case, is being used as the
  81 default to fill in the bonded and van der Waals :math:`\lambda` arrays.
  82 Usually, it’s best to fill in all arrays explicitly, just to make sure
  83 things are properly assigned.
  84
  85 If you want to turn on only restraints going from :math:`A` to
  86 :math:`B`, then it would be:
  87
  88 ::
  89
  90     restraint-lambdas      = 0.0 0.1 0.2 0.4 0.6 1.0
  91
  92 and all of the other components of the :math:`\lambda` vector would be
  93 left in the :math:`A` state.
  94
  95 To compute free energies with a vector :math:`\lambda` using
  96 thermodynamic integration, then the TI equation becomes vector equation:
  97
  98 .. math:: \Delta F = \int \langle \nabla H \rangle \cdot d\vec{\lambda}
  99           :label: eqnfepti
 100
 101 or for finite differences:
 102
 103 .. math:: \Delta F \approx \int \sum \langle \nabla H \rangle \cdot \Delta\lambda
 104           :label: eqnfepfinitediff
 105
 106 The external `pymbar script <https://SimTK.org/home/pymbar>`_
 107 can compute this integral automatically
 108 from the |Gromacs| ``dhdl.xvg`` output.
 109
 110 Potential of mean force
 111 -----------------------
 112
 113 A potential of mean force (PMF) is a potential that is obtained by
 114 integrating the mean force from an ensemble of configurations. In
 115 |Gromacs|, there are several different methods to calculate the mean
 116 force. Each method has its limitations, which are listed below.
 117
 118 -  **pull code:** between the centers of mass of molecules or groups of
 119    molecules.
 120
 121 -  **AWH code:** currently acts on coordinates provided by the pull
 122    code or the free-energy lambda parameter.
 123
 124 -  **free-energy code with harmonic bonds or constraints:** between
 125    single atoms.
 126
 127 -  **free-energy code with position restraints:** changing the
 128    conformation of a relatively immobile group of atoms.
 129
 130 -  **pull code in limited cases:** between groups of atoms that are part
 131    of a larger molecule for which the bonds are constrained with SHAKE
 132    or LINCS. If the pull group if relatively large, the pull code can be
 133    used.
 134
 135 The pull and free-energy code a described in more detail in the
 136 following two sections.
 137
 138 Entropic effects
 139 ^^^^^^^^^^^^^^^^
 140
 141 When a distance between two atoms or the centers of mass of two groups
 142 is constrained or restrained, there will be a purely entropic
 143 contribution to the PMF due to the rotation of the two
 144 groups \ :ref:`134 <refRMNeumann1980a>`. For a system of two
 145 non-interacting masses the potential of mean force is:
 146
 147 .. math:: V_{pmf}(r) = -(n_c - 1) k_B T \log(r)
 148           :label: eqnfepentropy
 149
 150 where :math:`n_c` is the number of dimensions in which the constraint
 151 works (i.e. :math:`n_c=3` for a normal constraint and :math:`n_c=1` when
 152 only the :math:`z`-direction is constrained). Whether one needs to
 153 correct for this contribution depends on what the PMF should represent.
 154 When one wants to pull a substrate into a protein, this entropic term
 155 indeed contributes to the work to get the substrate into the protein.
 156 But when calculating a PMF between two solutes in a solvent, for the
 157 purpose of simulating without solvent, the entropic contribution should
 158 be removed. **Note** that this term can be significant; when at 300K the
 159 distance is halved, the contribution is 3.5 kJ mol\ :math:`^{-1}`.