docs/reference-manual/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::
  40
  41      \begin{array}{rclcl}
  42      r       & \leq  & r_{HB}        & = & 0.35~\mbox{nm}    \\
  43      \alpha  & \leq  & \alpha_{HB}   & = & 30^o              \\
  44      \end{array}
  45
  46 .. _fig-hbond:
  47
  48 .. figure:: plots/hbond.*
  49    :width: 2.50000cm
  50
  51    Geometrical Hydrogen bond criterion.
  52
  53 The value of :math:`r_{HB} = 0.35 \mathrm{nm}` corresponds to the first minimum
  54 of the RDF of SPC water (see also :numref:`Fig. %s <fig-hbondinsert>`).
  55
  56 The program :ref:`gmx hbond <gmx hbond>` analyzes all hydrogen bonds
  57 existing between two groups of atoms (which must be either identical or
  58 non-overlapping) or in specified donor-hydrogen-acceptor triplets, in
  59 the following ways:
  60
  61 .. _fig-hbondinsert:
  62
  63 .. figure:: plots/hbond-insert.*
  64     :width: 3.50000cm
  65
  66     Insertion of water into an H-bond. (1) Normal H-bond between two
  67     residues. (2) H-bonding bridge via a water molecule.
  68
  69 -  Donor-Acceptor distance (:math:`r`) distribution of all H-bonds
  70
  71 -  Hydrogen-Donor-Acceptor angle (:math:`\alpha`) distribution of all
  72    H-bonds
  73
  74 -  The total number of H-bonds in each time frame
  75
  76 -  The number of H-bonds in time between residues, divided into groups
  77    :math:`n`-:math:`n`\ +\ :math:`i` where :math:`n` and
  78    :math:`n`\ +\ :math:`i` stand for residue numbers and :math:`i` goes
  79    from 0 to 6. The group for :math:`i=6` also includes all H-bonds for
  80    :math:`i>6`. These groups include the
  81    :math:`n`-:math:`n`\ +\ :math:`3`, :math:`n`-:math:`n`\ +\ :math:`4`
  82    and :math:`n`-:math:`n`\ +\ :math:`5` H-bonds, which provide a
  83    measure for the formation of :math:`\alpha`-helices or
  84    :math:`\beta`-turns or strands.
  85
  86 -  The lifetime of the H-bonds is calculated from the average over all
  87    autocorrelation functions of the existence functions (either 0 or 1)
  88    of all H-bonds:
  89
  90    .. math:: C(\tau) ~=~ \langle s_i(t)~s_i (t + \tau) \rangle
  91              :label: eqnhbcorr
  92
  93 -  with :math:`s_i(t) = \{0,1\}` for H-bond :math:`i` at time
  94    :math:`t`. The integral of :math:`C(\tau)` gives a rough estimate of
  95    the average H-bond lifetime :math:`\tau_{HB}`:
  96
  97    .. math::  \tau_{HB} ~=~ \int_{0}^{\infty} C(\tau) d\tau
  98               :label: eqnhblife
  99
 100 -  Both the integral and the complete autocorrelation function
 101    :math:`C(\tau)` will be output, so that more sophisticated analysis
 102    (*e.g.* using multi-exponential fits) can be used to get better
 103    estimates for :math:`\tau_{HB}`. A more complete analysis is given in
 104    ref. \ :ref:`173 <refSpoel2006b>`; one of the more fancy option is the Luzar
 105    and Chandler analysis of hydrogen bond kinetics \ :ref:`174 <refLuzar96b>`, :ref:`175 <refLuzar2000a>`.
 106
 107 -  An H-bond existence map can be generated of dimensions
 108    *# H-bonds*\ :math:`\times`\ *# frames*. The ordering is identical to
 109    the index file (see below), but reversed, meaning that the last
 110    triplet in the index file corresponds to the first row of the
 111    existence map.
 112
 113 -  Index groups are output containing the analyzed groups, all
 114    donor-hydrogen atom pairs and acceptor atoms in these groups,
 115    donor-hydrogen-acceptor triplets involved in hydrogen bonds between
 116    the analyzed groups and all solvent atoms involved in insertion.
 117