docs/reference-manual/analysis/dihedral-pca.rst

   1 Dihedral principal component analysis
   2 -------------------------------------
   3
   4 | :ref:`gmx angle <gmx angle>`, :ref:`gmx covar <gmx covar>`,
   5   :ref:`gmx anaeig <gmx anaeig>`
   6 | Principal component analysis can be performed in dihedral
   7   space \ :ref:`172 <refMu2005a>` using |Gromacs|. You start by defining the
   8   dihedral angles of interest in an index file, either using
   9   :ref:`gmx mk_angndx <gmx mk_angndx>` or otherwise. Then you use the
  10   :ref:`gmx angle <gmx angle>` program with the ``-or`` flag to
  11   produce a new :ref:`trr` file containing the cosine and sine
  12   of each dihedral angle in two coordinates, respectively. That is, in
  13   the :ref:`trr` file you will have a series of numbers
  14   corresponding to: cos(\ :math:`\phi_1`), sin(\ :math:`\phi_1`),
  15   cos(\ :math:`\phi_2`), sin(\ :math:`\phi_2`), ...,
  16   cos(\ :math:`\phi_n`), sin(\ :math:`\phi_n`), and the array is padded
  17   with zeros, if necessary. Then you can use this :ref:`trr`
  18   file as input for the :ref:`gmx covar <gmx covar>` program and perform
  19   principal component analysis as usual. For this to work you will need
  20   to generate a reference file (:ref:`tpr`,
  21   :ref:`gro`, :ref:`pdb` etc.) containing the same
  22   number of “atoms” as the new :ref:`trr` file, that is for
  23   :math:`n` dihedrals you need 2\ :math:`n`/3 atoms (rounded up if not
  24   an integer number). You should use the ``-nofit`` option
  25   for :ref:`gmx covar <gmx covar>` since the coordinates in the dummy
  26   reference file do not correspond in any way to the information in the
  27   :ref:`trr` file. Analysis of the results is done using
  28   :ref:`gmx anaeig <gmx anaeig>`.
  29
  30 Hydrogen bonds
  31 --------------
  32
  33 | :ref:`gmx hbond <gmx hbond>`
  34 | The program :ref:`gmx hbond <gmx hbond>`
  35   analyzes the *hydrogen bonds* (H-bonds) between all possible donors D
  36   and acceptors A. To determine if an H-bond exists, a geometrical
  37   criterion is used, see also :numref:`Fig. %s <fig-hbond>`:
  38
  39   .. math:: \begin{array}{rclcl}
  40             r       & \leq  & r_{HB}        & = & 0.35~\mbox{nm}    \\
  41             \alpha  & \leq  & \alpha_{HB}   & = & 30^o              \\
  42             \end{array}
  43             :label: eqnhbondgeomtric
  44
  45 .. _fig-hbond:
  46
  47 .. figure:: plots/hbond.*
  48    :width: 7.50000cm
  49
  50    Geometrical Hydrogen bond criterion.
  51
  52 The value of :math:`r_{HB} = 0.35 \mathrm{nm}` corresponds to the first minimum
  53 of the RDF of SPC water (see also :numref:`Fig. %s <fig-hbondinsert>`).
  54
  55 The program :ref:`gmx hbond <gmx hbond>` analyzes all hydrogen bonds
  56 existing between two groups of atoms (which must be either identical or
  57 non-overlapping) or in specified donor-hydrogen-acceptor triplets, in
  58 the following ways:
  59
  60 .. _fig-hbondinsert:
  61
  62 .. figure:: plots/hbond-insert.*
  63     :width: 7.50000cm
  64
  65     Insertion of water into an H-bond. (1) Normal H-bond between two
  66     residues. (2) H-bonding bridge via a water molecule.
  67
  68 -  Donor-Acceptor distance (:math:`r`) distribution of all H-bonds
  69
  70 -  Hydrogen-Donor-Acceptor angle (:math:`\alpha`) distribution of all
  71    H-bonds
  72
  73 -  The total number of H-bonds in each time frame
  74
  75 -  The number of H-bonds in time between residues, divided into groups
  76    :math:`n`-:math:`n`\ +\ :math:`i` where :math:`n` and
  77    :math:`n`\ +\ :math:`i` stand for residue numbers and :math:`i` goes
  78    from 0 to 6. The group for :math:`i=6` also includes all H-bonds for
  79    :math:`i>6`. These groups include the
  80    :math:`n`-:math:`n`\ +\ :math:`3`, :math:`n`-:math:`n`\ +\ :math:`4`
  81    and :math:`n`-:math:`n`\ +\ :math:`5` H-bonds, which provide a
  82    measure for the formation of :math:`\alpha`-helices or
  83    :math:`\beta`-turns or strands.
  84
  85 -  The lifetime of the H-bonds is calculated from the average over all
  86    autocorrelation functions of the existence functions (either 0 or 1)
  87    of all H-bonds:
  88
  89    .. math:: C(\tau) ~=~ \langle s_i(t)~s_i (t + \tau) \rangle
  90              :label: eqnhbcorr
  91
  92 -  with :math:`s_i(t) = \{0,1\}` for H-bond :math:`i` at time
  93    :math:`t`. The integral of :math:`C(\tau)` gives a rough estimate of
  94    the average H-bond lifetime :math:`\tau_{HB}`:
  95
  96    .. math::  \tau_{HB} ~=~ \int_{0}^{\infty} C(\tau) d\tau
  97               :label: eqnhblife
  98
  99 -  Both the integral and the complete autocorrelation function
 100    :math:`C(\tau)` will be output, so that more sophisticated analysis
 101    (*e.g.* using multi-exponential fits) can be used to get better
 102    estimates for :math:`\tau_{HB}`. A more complete analysis is given in
 103    ref. \ :ref:`173 <refSpoel2006b>`; one of the more fancy option is the Luzar
 104    and Chandler analysis of hydrogen bond kinetics \ :ref:`174 <refLuzar96b>`, :ref:`175 <refLuzar2000a>`.
 105
 106 -  An H-bond existence map can be generated of dimensions
 107    *# H-bonds*\ :math:`\times`\ *# frames*. The ordering is identical to
 108    the index file (see below), but reversed, meaning that the last
 109    triplet in the index file corresponds to the first row of the
 110    existence map.
 111
 112 -  Index groups are output containing the analyzed groups, all
 113    donor-hydrogen atom pairs and acceptor atoms in these groups,
 114    donor-hydrogen-acceptor triplets involved in hydrogen bonds between
 115    the analyzed groups and all solvent atoms involved in insertion.
 116