docs/reference-manual/analysis/correlation-function.rst

   1 Correlation functions
   2 ---------------------
   3
   4 Theory of correlation functions
   5 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   6
   7 The theory of correlation functions is well established \ :ref:`108 <refAllen87>`.
   8 We describe here the implementation of the various
   9 correlation function flavors in the |Gromacs| code. The definition of the
  10 autocorrelation function (ACF) :math:`C_f(t)` for a property
  11 :math:`f(t)` is:
  12
  13 .. math:: C_f(t)  ~=~     \left\langle f(\xi) f(\xi+t)\right\rangle_{\xi}
  14           :label: eqncorr
  15
  16 where the notation on the right hand side indicates averaging over
  17 :math:`\xi`, *i.e.* over time origins. It is also possible to compute
  18 cross-correlation function from two properties :math:`f(t)` and
  19 :math:`g(t)`:
  20
  21 .. math:: C_{fg}(t) ~=~   \left\langle f(\xi) g(\xi+t)\right\rangle_{\xi}
  22
  23 however, in |Gromacs| there is no standard mechanism to do this
  24 (**note:** you can use the ``xmgr`` program to compute cross correlations).
  25 The integral of the correlation function over time is the correlation
  26 time :math:`\tau_f`:
  27
  28 .. math:: \tau_f  ~=~     \int_0^{\infty} C_f(t) {\rm d} t
  29           :label: eqncorrtime
  30
  31 In practice, correlation functions are calculated based on data points
  32 with discrete time intervals :math:`\Delta`\ t, so that the ACF from an
  33 MD simulation is:
  34
  35 .. math:: C_f(j\Delta t)  ~=~     \frac{1}{N-j}\sum_{i=0}^{N-1-j} f(i\Delta t) f((i+j)\Delta t)
  36           :label: eqncorrmd
  37
  38 where :math:`N` is the number of available time frames for the
  39 calculation. The resulting ACF is obviously only available at time
  40 points with the same interval :math:`\Delta`\ t. Since, for many
  41 applications, it is necessary to know the short time behavior of the ACF
  42 (*e.g.* the first 10 ps) this often means that we have to save the data
  43 with intervals much shorter than the time scale of interest. Another
  44 implication of :eq:`eqn. %s <eqncorrmd>` is that in principle we can not compute
  45 all points of the ACF with the same accuracy, since we have :math:`N-1`
  46 data points for :math:`C_f(\Delta t)` but only 1 for
  47 :math:`C_f((N-1)\Delta t)`. However, if we decide to compute only an ACF
  48 of length :math:`M\Delta t`, where :math:`M \leq N/2` we can compute all
  49 points with the same statistical accuracy:
  50
  51 .. math:: C_f(j\Delta t)  ~=~ \frac{1}{M}\sum_{i=0}^{N-1-M} f(i\Delta t)f((i+j)\Delta t)
  52
  53 Here of course :math:`j < M`. :math:`M` is sometimes referred to as the
  54 time lag of the correlation function. When we decide to do this, we
  55 intentionally do not use all the available points for very short time
  56 intervals (:math:`j << M`), but it makes it easier to interpret the
  57 results. Another aspect that may not be neglected when computing ACFs
  58 from simulation is that usually the time origins :math:`\xi`
  59 (:eq:`eqn. %s <eqncorr>`) are not statistically independent, which may introduce
  60 a bias in the results. This can be tested using a block-averaging
  61 procedure, where only time origins with a spacing at least the length of
  62 the time lag are included, *e.g.* using :math:`k` time origins with
  63 spacing of :math:`M\Delta t` (where :math:`kM \leq N`):
  64
  65 .. math:: C_f(j\Delta t)  ~=~ \frac{1}{k}\sum_{i=0}^{k-1} f(iM\Delta t)f((iM+j)\Delta t)
  66
  67 However, one needs very long simulations to get good accuracy this way,
  68 because there are many fewer points that contribute to the ACF.
  69
  70 Using FFT for computation of the ACF
  71 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  72
  73 The computational cost for calculating an ACF according to
  74 :eq:`eqn. %s <eqncorrmd>` is proportional to :math:`N^2`, which is considerable.
  75 However, this can be improved by using fast Fourier transforms to do the
  76 convolution \ :ref:`108 <refAllen87>`.
  77
  78 Special forms of the ACF
  79 ~~~~~~~~~~~~~~~~~~~~~~~~
  80
  81 There are some important varieties on the ACF, *e.g.* the ACF of a
  82 vector :math:`\mathbf{p}`:
  83
  84 .. math:: C_{\mathbf{p}}(t) ~=~       \int_0^{\infty} P_n(\cos\angle\left(\mathbf{p}(\xi),\mathbf{p}(\xi+t)\right) {\rm d} \xi
  85           :label: eqncorrleg
  86
  87 where :math:`P_n(x)` is the :math:`n^{th}` order Legendre
  88 polynomial. [1]_ Such correlation times can actually be obtained
  89 experimentally using *e.g.* NMR or other relaxation experiments. |Gromacs|
  90 can compute correlations using the 1\ :math:`^{st}` and 2\ :math:`^{nd}`
  91 order Legendre polynomial (:eq:`eqn. %s <eqncorrleg>`). This can also be used
  92 for rotational autocorrelation (:ref:`gmx rotacf`) and dipole autocorrelation
  93 (:ref:`gmx dipoles <gmx dipoles>`).
  94
  95 In order to study torsion angle dynamics, we define a dihedral
  96 autocorrelation function as \ :ref:`159 <refSpoel97a>`:
  97
  98 .. math:: C(t)    ~=~     \left\langle \cos(\theta(\tau)-\theta(\tau+t))\right\rangle_{\tau}
  99           :label: eqncoenk
 100
 101 **Note** that this is not a product of two functions as is generally
 102 used for correlation functions, but it may be rewritten as the sum of
 103 two products:
 104
 105 .. math:: C(t)    ~=~     \left\langle\cos(\theta(\tau))\cos(\theta(\tau+t))\,+\,\sin(\theta(\tau))\sin(\theta(\tau+t))\right\rangle_{\tau}
 106           :label: eqncot
 107
 108 Some Applications
 109 ~~~~~~~~~~~~~~~~~
 110
 111 The program :ref:`gmx velacc <gmx velacc>`
 112 calculates the *velocity autocorrelation function*.
 113
 114 .. math:: C_{\mathbf{v}} (\tau) ~=~ \langle {\mathbf{v}}_i(\tau) \cdot {\mathbf{v}}_i(0) \rangle_{i \in A}
 115
 116 The self diffusion coefficient can be calculated using the Green-Kubo
 117 relation \ :ref:`108 <refAllen87>`:
 118
 119 .. math:: D_A ~=~ {1\over 3} \int_0^{\infty} \langle {\bf v}_i(t) \cdot {\bf v}_i(0) \rangle_{i \in A} \; dt
 120
 121 which is just the integral of the velocity autocorrelation function.
 122 There is a widely-held belief that the velocity ACF converges faster
 123 than the mean square displacement (sec. :ref:`msd`), which can also be
 124 used for the computation of diffusion constants. However, Allen &
 125 Tildesley \ :ref:`108 <refAllen87>` warn us that the long-time
 126 contribution to the velocity ACF can not be ignored, so care must be
 127 taken.
 128
 129 Another important quantity is the dipole correlation time. The *dipole
 130 correlation function* for particles of type :math:`A` is calculated as
 131 follows by :ref:`gmx dipoles <gmx dipoles>`:
 132
 133 .. math::
 134
 135    C_{\mu} (\tau) ~=~
 136    \langle {\bf \mu}_i(\tau) \cdot {\bf \mu}_i(0) \rangle_{i \in A}
 137
 138 with :math:`{\bf \mu}_i = \sum_{j \in i} {\bf r}_j q_j`. The dipole
 139 correlation time can be computed using :eq:`eqn. %s <eqncorrtime>`.
 140 For some applications
 141 see (**???**).
 142
 143 The viscosity of a liquid can be related to the correlation time of the
 144 Pressure tensor
 145 :math:`\mathbf{P}` :ref:`160 <refPSmith93c>`,
 146 :ref:`161 <refBalasubramanian96>`. :ref:`gmx energy` can compute the viscosity,
 147 but this is not very accurate \ :ref:`149 <refHess2002a>`, and actually
 148 the values do not converge.
 149
 150
 151 .. [1]
 152    :math:`P_0(x) = 1`, :math:`P_1(x) = x`, :math:`P_2(x) = (3x^2-1)/2`
 153