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