docs/reference-manual/replica-exchange.rst

   1 Replica exchange
   2 ----------------
   3
   4 Replica exchange molecular dynamics (REMD) is a method that can be used
   5 to speed up the sampling of any type of simulation, especially if
   6 conformations are separated by relatively high energy barriers. It
   7 involves simulating multiple replicas of the same system at different
   8 temperatures and randomly exchanging the complete state of two replicas
   9 at regular intervals with the probability:
  10
  11 .. math::
  12
  13    P(1 \leftrightarrow 2)=\min\left(1,\exp\left[
  14    \left(\frac{1}{k_B T_1} - \frac{1}{k_B T_2}\right)(U_1 - U_2)
  15     \right] \right)
  16
  17 where :math:`T_1` and :math:`T_2` are the reference temperatures and
  18 :math:`U_1` and :math:`U_2` are the instantaneous potential energies of
  19 replicas 1 and 2 respectively. After exchange the velocities are scaled
  20 by :math:`(T_1/T_2)^{\pm0.5}` and a neighbor search is performed the
  21 next step. This combines the fast sampling and frequent barrier-crossing
  22 of the highest temperature with correct Boltzmann sampling at all the
  23 different temperatures \ :ref:`60 <refHukushima96a>`,
  24 :ref:`61 <refSugita99>`. We only attempt exchanges for neighboring temperatures as the
  25 probability decreases very rapidly with the temperature difference. One
  26 should not attempt exchanges for all possible pairs in one step. If, for
  27 instance, replicas 1 and 2 would exchange, the chance of exchange for
  28 replicas 2 and 3 not only depends on the energies of replicas 2 and 3,
  29 but also on the energy of replica 1. In |Gromacs| this is solved by
  30 attempting exchange for all *odd* pairs on *odd* attempts and for all
  31 *even* pairs on *even* attempts. If we have four replicas: 0, 1, 2 and
  32 3, ordered in temperature and we attempt exchange every 1000 steps,
  33 pairs 0-1 and 2-3 will be tried at steps 1000, 3000 etc. and pair 1-2 at
  34 steps 2000, 4000 etc.
  35
  36 How should one choose the temperatures? The energy difference can be
  37 written as:
  38
  39 .. math:: U_1 - U_2 =  N_{df} \frac{c}{2} k_B (T_1 - T_2)
  40
  41 where :math:`N_{df}` is the total number of degrees of freedom of one
  42 replica and :math:`c` is 1 for harmonic potentials and around 2 for
  43 protein/water systems. If :math:`T_2 = (1+\epsilon) T_1` the probability
  44 becomes:
  45
  46 .. math::
  47
  48    P(1 \leftrightarrow 2)
  49      = \exp\left( -\frac{\epsilon^2 c\,N_{df}}{2 (1+\epsilon)} \right)
  50    \approx \exp\left(-\epsilon^2 \frac{c}{2} N_{df} \right)
  51
  52 Thus for a probability of :math:`e^{-2}\approx 0.135` one obtains
  53 :math:`\epsilon \approx 2/\sqrt{c\,N_{df}}`. With all bonds constrained
  54 one has :math:`N_{df} \approx 2\, N_{atoms}` and thus for :math:`c` = 2
  55 one should choose :math:`\epsilon` as :math:`1/\sqrt{N_{atoms}}`.
  56 However there is one problem when using pressure coupling. The density
  57 at higher temperatures will decrease, leading to higher energy
  58 \ :ref:`62 <refSeibert2005a>`, which should be taken into account. The |Gromacs| website
  59 features a so-called ``REMD calculator``, that lets you type in the
  60 temperature range and the number of atoms, and based on that proposes a
  61 set of temperatures.
  62
  63 An extension to the REMD for the isobaric-isothermal ensemble was
  64 proposed by Okabe et al. :ref:`63 <refOkabe2001a>`. In this work the
  65 exchange probability is modified to:
  66
  67 .. math::
  68
  69    P(1 \leftrightarrow 2)=\min\left(1,\exp\left[
  70    \left(\frac{1}{k_B T_1} - \frac{1}{k_B T_2}\right)(U_1 - U_2) +
  71    \left(\frac{P_1}{k_B T_1} - \frac{P_2}{k_B T_2}\right)\left(V_1-V_2\right)
  72     \right] \right)
  73
  74 where :math:`P_1` and :math:`P_2` are the respective reference
  75 pressures and :math:`V_1` and :math:`V_2` are the respective
  76 instantaneous volumes in the simulations. In most cases the differences
  77 in volume are so small that the second term is negligible. It only plays
  78 a role when the difference between :math:`P_1` and :math:`P_2` is large
  79 or in phase transitions.
  80
  81 Hamiltonian replica exchange is also supported in |Gromacs|. In
  82 Hamiltonian replica exchange, each replica has a different Hamiltonian,
  83 defined by the free energy pathway specified for the simulation. The
  84 exchange probability to maintain the correct ensemble probabilities is:
  85
  86 .. math::
  87
  88    P(1 \leftrightarrow 2)=\min\left(1,\exp\left[
  89        \left(\frac{1}{k_B T} - \frac{1}{k_B T}\right)((U_1(x_2) - U_1(x_1)) + (U_2(x_1) - U_2(x_2)))
  90    \right]
  91    \right)
  92
  93 The separate Hamiltonians are defined by the free energy functionality
  94 of |Gromacs|, with swaps made between the different values of
  95 :math:`\lambda` defined in the mdp file.
  96
  97 Hamiltonian and temperature replica exchange can also be performed
  98 simultaneously, using the acceptance criteria:
  99
 100 .. math::
 101
 102    P(1 \leftrightarrow 2)=\min\left(1,\exp\left[
 103    \left(\frac{1}{k_B T} - \right)(\frac{U_1(x_2) - U_1(x_1)}{k_B T_1} + \frac{U_2(x_1) - U_2(x_2)}{k_B T_2})
 104     \right] \right)
 105
 106 Gibbs sampling replica exchange has also been implemented in
 107 |Gromacs| :ref:`64 <refChodera2011>`. In Gibbs sampling replica exchange,
 108 all possible pairs are tested for exchange, allowing swaps between
 109 replicas that are not neighbors.
 110
 111 Gibbs sampling replica exchange requires no additional potential energy
 112 calculations. However there is an additional communication cost in Gibbs
 113 sampling replica exchange, as for some permutations, more than one round
 114 of swaps must take place. In some cases, this extra communication cost
 115 might affect the efficiency.
 116
 117 All replica exchange variants are options of the :ref:`mdrun <gmx mdrun>` program. It will
 118 only work when MPI is installed, due to the inherent parallelism in the
 119 algorithm. For efficiency each replica can run on a separate rank. See
 120 the manual page of :ref:`mdrun <gmx mdrun>` on how to use these multinode features.