docs/reference-manual/covariance-analysis.rst

   1 .. _covanal:
   2
   3 Covariance analysis
   4 -------------------
   5
   6 Covariance analysis, also called principal component analysis or
   7 essential dynamics :ref:`169 <refAmadei93>`\ , can find
   8 correlated motions. It uses the covariance matrix :math:`C` of the
   9 atomic coordinates:
  10
  11 .. math::
  12
  13    C_{ij} = \left \langle
  14    M_{ii}^{\frac{1}{2}} (x_i - \langle x_i \rangle)
  15    M_{jj}^{\frac{1}{2}}  (x_j - \langle x_j \rangle)
  16    \right \rangle
  17
  18 where :math:`M` is a diagonal matrix containing the masses of the atoms
  19 (mass-weighted analysis) or the unit matrix (non-mass weighted
  20 analysis). :math:`C` is a symmetric :math:`3N \times 3N` matrix, which
  21 can be diagonalized with an orthonormal transformation matrix :math:`R`:
  22
  23 .. math::
  24
  25    R^T C R = \mbox{diag}(\lambda_1,\lambda_2,\ldots,\lambda_{3N})
  26    ~~~~\mbox{where}~~\lambda_1 \geq \lambda_2 \geq \ldots \geq \lambda_{3N}
  27
  28 The columns of :math:`R` are the eigenvectors, also called principal or
  29 essential modes. :math:`R` defines a transformation to a new coordinate
  30 system. The trajectory can be projected on the principal modes to give
  31 the principal components :math:`p_i(t)`:
  32
  33 .. math:: {\bf p}(t) = R^T M^{\frac{1}{2}} ({\bf x}(t) - \langle {\bf x} \rangle)
  34
  35 The eigenvalue :math:`\lambda_i` is the mean square fluctuation of
  36 principal component :math:`i`. The first few principal modes often
  37 describe collective, global motions in the system. The trajectory can be
  38 filtered along one (or more) principal modes. For one principal mode
  39 :math:`i` this goes as follows:
  40
  41 .. math::
  42
  43    {\bf x}^f(t) =
  44    \langle {\bf x} \rangle + M^{-\frac{1}{2}} R_{ * i} \, p_i(t)
  45
  46 When the analysis is performed on a macromolecule, one often wants to
  47 remove the overall rotation and translation to look at the internal
  48 motion only. This can be achieved by least square fitting to a reference
  49 structure. Care has to be taken that the reference structure is
  50 representative for the ensemble, since the choice of reference structure
  51 influences the covariance matrix.
  52
  53 One should always check if the principal modes are well defined. If the
  54 first principal component resembles a half cosine and the second
  55 resembles a full cosine, you might be filtering noise (see below). A
  56 good way to check the relevance of the first few principal modes is to
  57 calculate the overlap of the sampling between the first and second half
  58 of the simulation. **Note** that this can only be done when the same
  59 reference structure is used for the two halves.
  60
  61 A good measure for the overlap has been defined in \ :ref:`170 <refHess2002b>`. The
  62 elements of the covariance matrix are proportional to the square of the
  63 displacement, so we need to take the square root of the matrix to
  64 examine the extent of sampling. The square root can be calculated from
  65 the eigenvalues :math:`\lambda_i` and the eigenvectors, which are the
  66 columns of the rotation matrix :math:`R`. For a symmetric and
  67 diagonally-dominant matrix :math:`A` of size :math:`3N \times 3N` the
  68 square root can be calculated as:
  69
  70 .. math::
  71
  72    A^\frac{1}{2} =
  73    R \, \mbox{diag}(\lambda_1^\frac{1}{2},\lambda_2^\frac{1}{2},\ldots,\lambda_{3N}^\frac{1}{2}) \, R^T
  74
  75 It can be verified easily that the product of this matrix with itself
  76 gives :math:`A`. Now we can define a difference :math:`d` between
  77 covariance matrices :math:`A` and :math:`B` as follows:
  78
  79 .. math::
  80
  81    \begin{aligned}
  82    d(A,B) & = & \sqrt{\mbox{tr}\left(\left(A^\frac{1}{2} - B^\frac{1}{2}\right)^2\right)
  83    }
  84    \\ & = &
  85    \sqrt{\mbox{tr}\left(A + B - 2 A^\frac{1}{2} B^\frac{1}{2}\right)}
  86    \\ & = &
  87    \left( \sum_{i=1}^N \left( \lambda_i^A + \lambda_i^B \right)
  88    - 2 \sum_{i=1}^N \sum_{j=1}^N \sqrt{\lambda_i^A \lambda_j^B}
  89    \left(R_i^A \cdot R_j^B\right)^2 \right)^\frac{1}{2}\end{aligned}
  90
  91 where tr is the trace of a matrix. We can now define the overlap
  92 :math:`s` as:
  93
  94 .. math:: s(A,B) = 1 - \frac{d(A,B)}{\sqrt{\mbox{tr}A + \mbox{tr} B}}
  95
  96 The overlap is 1 if and only if matrices :math:`A` and :math:`B` are
  97 identical. It is 0 when the sampled subspaces are completely orthogonal.
  98
  99 A commonly-used measure is the subspace overlap of the first few
 100 eigenvectors of covariance matrices. The overlap of the subspace spanned
 101 by :math:`m` orthonormal vectors :math:`{\bf w}_1,\ldots,{\bf w}_m` with
 102 a reference subspace spanned by :math:`n` orthonormal vectors
 103 :math:`{\bf v}_1,\ldots,{\bf v}_n` can be quantified as follows:
 104
 105 .. math::
 106
 107    \mbox{overlap}({\bf v},{\bf w}) =
 108    \frac{1}{n} \sum_{i=1}^n \sum_{j=1}^m ({\bf v}_i \cdot {\bf w}_j)^2
 109
 110 The overlap will increase with increasing :math:`m` and will be 1 when
 111 set :math:`{\bf v}` is a subspace of set :math:`{\bf w}`. The
 112 disadvantage of this method is that it does not take the eigenvalues
 113 into account. All eigenvectors are weighted equally, and when degenerate
 114 subspaces are present (equal eigenvalues), the calculated overlap will
 115 be too low.
 116
 117 Another useful check is the cosine content. It has been proven that the
 118 the principal components of random diffusion are cosines with the number
 119 of periods equal to half the principal component
 120 index \ :ref:`170 <refHess2002b>`, :ref:`171 <refHess2000>`.
 121 The eigenvalues are proportional to the index to the power
 122 :math:`-2`. The cosine content is defined as:
 123
 124 .. math::
 125
 126    \frac{2}{T}
 127    \left( \int_0^T \cos\left(\frac{i \pi t}{T}\right) \, p_i(t) \mbox{d} t \right)^2
 128    \left( \int_0^T p_i^2(t) \mbox{d} t \right)^{-1}
 129
 130 When the cosine content of the first few principal components is close
 131 to 1, the largest fluctuations are not connected with the potential, but
 132 with random diffusion.
 133
 134 The covariance matrix is built and diagonalized by
 135 :ref:`gmx covar <gmx covar>`. The principal components and
 136 overlap (and many more things) can be plotted and analyzed with
 137 :ref:`gmx anaeig <gmx anaeig>`. The cosine
 138 content can be calculated with
 139 :ref:`gmx analyze <gmx analyze>`.