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 For example, if you wanted to first to change the Coulombic terms, then
  51 the van der Waals terms, changing bonded at the same time rate as the
  52 van der Waals, but changing the restraints throughout the first
  53 two-thirds of the simulation, then you could use this :math:`\lambda`
  54 vector:
  55
  56 ::
  57
  58     coul-lambdas           = 0.0 0.2 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0
  59     vdw-lambdas            = 0.0 0.0 0.0 0.0 0.4 0.5 0.6 0.7 0.8 1.0
  60     bonded-lambdas         = 0.0 0.0 0.0 0.0 0.4 0.5 0.6 0.7 0.8 1.0
  61     restraint-lambdas      = 0.0 0.0 0.1 0.2 0.3 0.5 0.7 1.0 1.0 1.0
  62
  63 This is also equivalent to:
  64
  65 ::
  66
  67     fep-lambdas            = 0.0 0.0 0.0 0.0 0.4 0.5 0.6 0.7 0.8 1.0
  68     coul-lambdas           = 0.0 0.2 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0
  69     restraint-lambdas      = 0.0 0.0 0.1 0.2 0.3 0.5 0.7 1.0 1.0 1.0
  70
  71 The ``fep-lambda`` array, in this case, is being used as the
  72 default to fill in the bonded and van der Waals :math:`\lambda` arrays.
  73 Usually, it’s best to fill in all arrays explicitly, just to make sure
  74 things are properly assigned.
  75
  76 If you want to turn on only restraints going from :math:`A` to
  77 :math:`B`, then it would be:
  78
  79 ::
  80
  81     restraint-lambdas      = 0.0 0.1 0.2 0.4 0.6 1.0
  82
  83 and all of the other components of the :math:`\lambda` vector would be
  84 left in the :math:`A` state.
  85
  86 To compute free energies with a vector :math:`\lambda` using
  87 thermodynamic integration, then the TI equation becomes vector equation:
  88
  89 .. math:: \Delta F = \int \langle \nabla H \rangle \cdot d\vec{\lambda}
  90
  91 or for finite differences:
  92
  93 .. math:: \Delta F \approx \int \sum \langle \nabla H \rangle \cdot \Delta\lambda
  94
  95 The external `pymbar script <https://SimTK.org/home/pymbar>`_
  96 can compute this integral automatically
  97 from the |Gromacs| ``dhdl.xvg`` output.
  98
  99 Potential of mean force
 100 -----------------------
 101
 102 A potential of mean force (PMF) is a potential that is obtained by
 103 integrating the mean force from an ensemble of configurations. In
 104 |Gromacs|, there are several different methods to calculate the mean
 105 force. Each method has its limitations, which are listed below.
 106
 107 -  **pull code:** between the centers of mass of molecules or groups of
 108    molecules.
 109
 110 -  **AWH code:** currently acts on coordinates provided by the pull
 111    code.
 112
 113 -  **free-energy code with harmonic bonds or constraints:** between
 114    single atoms.
 115
 116 -  **free-energy code with position restraints:** changing the
 117    conformation of a relatively immobile group of atoms.
 118
 119 -  **pull code in limited cases:** between groups of atoms that are part
 120    of a larger molecule for which the bonds are constrained with SHAKE
 121    or LINCS. If the pull group if relatively large, the pull code can be
 122    used.
 123
 124 The pull and free-energy code a described in more detail in the
 125 following two sections.
 126
 127 Entropic effects
 128 ^^^^^^^^^^^^^^^^
 129
 130 When a distance between two atoms or the centers of mass of two groups
 131 is constrained or restrained, there will be a purely entropic
 132 contribution to the PMF due to the rotation of the two
 133 groups \ :ref:`134 <refRMNeumann1980a>`. For a system of two
 134 non-interacting masses the potential of mean force is:
 135
 136 .. math:: V_{pmf}(r) = -(n_c - 1) k_B T \log(r)
 137
 138 where :math:`n_c` is the number of dimensions in which the constraint
 139 works (i.e. :math:`n_c=3` for a normal constraint and :math:`n_c=1` when
 140 only the :math:`z`-direction is constrained). Whether one needs to
 141 correct for this contribution depends on what the PMF should represent.
 142 When one wants to pull a substrate into a protein, this entropic term
 143 indeed contributes to the work to get the substrate into the protein.
 144 But when calculating a PMF between two solutes in a solvent, for the
 145 purpose of simulating without solvent, the entropic contribution should
 146 be removed. **Note** that this term can be significant; when at 300K the
 147 distance is halved, the contribution is 3.5 kJ mol\ :math:`^{-1}`.