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