docs/reference-manual/special/remove-fast-dgf.rst

   1 .. _rmfast:
   2
   3 Removing fastest degrees of freedom
   4 -----------------------------------
   5
   6 The maximum time step in MD simulations is limited by the smallest
   7 oscillation period that can be found in the simulated system.
   8 Bond-stretching vibrations are in their quantum-mechanical ground state
   9 and are therefore better represented by a constraint instead of a
  10 harmonic potential.
  11
  12 For the remaining degrees of freedom, the shortest oscillation period
  13 (as measured from a simulation) is 13 fs for bond-angle vibrations
  14 involving hydrogen atoms. Taking as a guideline that with a Verlet
  15 (leap-frog) integration scheme a minimum of 5 numerical integration
  16 steps should be performed per period of a harmonic oscillation in order
  17 to integrate it with reasonable accuracy, the maximum time step will be
  18 about 3 fs. Disregarding these very fast oscillations of period 13 fs,
  19 the next shortest periods are around 20 fs, which will allow a maximum
  20 time step of about 4 fs.
  21
  22 Removing the bond-angle degrees of freedom from hydrogen atoms can best
  23 be done by defining them as virtual interaction sites instead of normal
  24 atoms. Whereas a normal atom is connected to the molecule with bonds,
  25 angles and dihedrals, a virtual site’s position is calculated from the
  26 position of three nearby heavy atoms in a predefined manner (see also
  27 sec. :ref:`virtualsites`). For the hydrogens in water and in hydroxyl,
  28 sulfhydryl, or amine groups, no degrees of freedom can be removed,
  29 because rotational freedom should be preserved. The only other option
  30 available to slow down these motions is to increase the mass of the
  31 hydrogen atoms at the expense of the mass of the connected heavy atom.
  32 This will increase the moment of inertia of the water molecules and the
  33 hydroxyl, sulfhydryl, or amine groups, without affecting the equilibrium
  34 properties of the system and without affecting the dynamical properties
  35 too much. These constructions will shortly be described in
  36 sec. :ref:`vsitehydro` and have previously been described in full
  37 detail \ :ref:`148 <reffeenstra99>`.
  38
  39 Using both virtual sites and modified masses, the next bottleneck is
  40 likely to be formed by the improper dihedrals (which are used to
  41 preserve planarity or chirality of molecular groups) and the peptide
  42 dihedrals. The peptide dihedral cannot be changed without affecting the
  43 physical behavior of the protein. The improper dihedrals that preserve
  44 planarity mostly deal with aromatic residues. Bonds, angles, and
  45 dihedrals in these residues can also be replaced with somewhat elaborate
  46 virtual site constructions.
  47
  48 All modifications described in this section can be performed using the
  49 |Gromacs| topology building tool :ref:`pdb2gmx <gmx pdb2gmx>`. Separate options exist to
  50 increase hydrogen masses, virtualize all hydrogen atoms, or also
  51 virtualize all aromatic residues. **Note** that when all hydrogen atoms
  52 are virtualized, those inside the aromatic residues will be virtualized
  53 as well, *i.e.* hydrogens in the aromatic residues are treated
  54 differently depending on the treatment of the aromatic residues.
  55
  56 Parameters for the virtual site constructions for the hydrogen atoms are
  57 inferred from the force-field parameters (*vis*. bond lengths and
  58 angles) directly by :ref:`grompp <gmx grompp>` while processing the topology file. The
  59 constructions for the aromatic residues are based on the bond lengths
  60 and angles for the geometry as described in the force fields, but these
  61 parameters are hard-coded into :ref:`pdb2gmx <gmx pdb2gmx>` due to the complex nature of the
  62 construction needed for a whole aromatic group.
  63
  64 .. _vsitehydro:
  65
  66 Hydrogen bond-angle vibrations
  67 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  68
  69 Construction of virtual sites
  70 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  71
  72 .. _fig-vsitehydro:
  73
  74 .. figure:: plots/dumtypes.*
  75    :width: 11.00000cm
  76
  77    The different types of virtual site constructions used for hydrogen
  78    atoms. The atoms used in the construction of the virtual site(s) are
  79    depicted as black circles, virtual sites as gray ones. Hydrogens are
  80    smaller than heavy atoms. A: fixed bond angle, note
  81    that here the hydrogen is not a virtual site; B: in
  82    the plane of three atoms, with fixed distance; C: in
  83    the plane of three atoms, with fixed angle and distance;
  84    D: construction for amine groups
  85    (-NH:math:`_2` or -NH:math:`_3^+`),
  86    see text for details.
  87
  88 The goal of defining hydrogen atoms as virtual sites is to remove all
  89 high-frequency degrees of freedom from them. In some cases, not all
  90 degrees of freedom of a hydrogen atom should be removed, *e.g.* in the
  91 case of hydroxyl or amine groups the rotational freedom of the hydrogen
  92 atom(s) should be preserved. Care should be taken that no unwanted
  93 correlations are introduced by the construction of virtual sites, *e.g.*
  94 bond-angle vibration between the constructing atoms could translate into
  95 hydrogen bond-length vibration. Additionally, since virtual sites are by
  96 definition massless, in order to preserve total system mass, the mass of
  97 each hydrogen atom that is treated as virtual site should be added to
  98 the bonded heavy atom.
  99
 100 Taking into account these considerations, the hydrogen atoms in a
 101 protein naturally fall into several categories, each requiring a
 102 different approach (see also :numref:`Fig. %s <fig-vsitehydro>`).
 103
 104 -  *hydroxyl (-OH) or sulfhydryl (-SH) hydrogen:*
 105    The only internal degree of freedom in a hydroxyl group
 106    that can be constrained is the bending of the C-O-H
 107    angle. This angle is fixed by defining an additional bond of
 108    appropriate length, see :numref:`Fig. %s A<fig-vsitehydro>`.
 109    Doing so removes the high-frequency angle bending, but leaves the
 110    dihedral rotational freedom. The same goes for a sulfhydryl group.
 111    **Note** that in these cases the hydrogen is not treated as a virtual
 112    site.
 113
 114 -  *single amine or amide (-NH-) and aromatic hydrogens
 115    (-CH-):*
 116    The position of these hydrogens cannot be
 117    constructed from a linear combination of bond vectors, because of the
 118    flexibility of the angle between the heavy atoms. Instead, the
 119    hydrogen atom is positioned at a fixed distance from the bonded heavy
 120    atom on a line going through the bonded heavy atom and a point on the
 121    line through both second bonded atoms, see
 122    :numref:`Fig. %s B<fig-vsitehydro>`.
 123
 124 -  *planar amine (-NH*:math:`_2`) *hydrogens:* The method
 125    used for the single amide hydrogen is not well suited for planar
 126    amine groups, because no suitable two heavy atoms can be found to
 127    define the direction of the hydrogen atoms. Instead, the hydrogen is
 128    constructed at a fixed distance from the nitrogen atom, with a fixed
 129    angle to the carbon atom, in the plane defined by one of the other
 130    heavy atoms, see :numref:`Fig. %s C<fig-vsitehydro>`.
 131
 132 -  *amine group (umbrella -NH*:math:`_2` *or
 133    -NH*:math:`_3^+`)* hydrogens:* Amine hydrogens with
 134    rotational freedom cannot be constructed as virtual sites from the
 135    heavy atoms they are connected to, since this would result in loss of
 136    the rotational freedom of the amine group. To preserve the rotational
 137    freedom while removing the hydrogen bond-angle degrees of freedom,
 138    two “dummy masses” are constructed with the same total mass, moment
 139    of inertia (for rotation around the C-N bond) and
 140    center of mass as the amine group. These dummy masses have no
 141    interaction with any other atom, except for the fact that they are
 142    connected to the carbon and to each other, resulting in a rigid
 143    triangle. From these three particles, the positions of the nitrogen
 144    and hydrogen atoms are constructed as linear combinations of the two
 145    carbon-mass vectors and their outer product, resulting in an amine
 146    group with rotational freedom intact, but without other internal
 147    degrees of freedom. See :numref:`Fig. %s D<fig-vsitehydro>`.
 148
 149 .. figure:: plots/dumaro.*
 150    :width: 15.00000cm
 151
 152    The different types of virtual site constructions used for aromatic
 153    residues. The atoms used in the construction of the virtual site(s)
 154    are depicted as black circles, virtual sites as gray ones. Hydrogens
 155    are smaller than heavy atoms. A: phenylalanine;
 156    B: tyrosine (note that the hydroxyl hydrogen is *not*
 157    a virtual site); C: tryptophan; D:
 158    histidine.
 159
 160 Out-of-plane vibrations in aromatic groups
 161 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 162
 163 The planar arrangements in the side chains of the aromatic residues
 164 lends itself perfectly to a virtual-site construction, giving a
 165 perfectly planar group without the inherently unstable constraints that
 166 are necessary to keep normal atoms in a plane. The basic approach is to
 167 define three atoms or dummy masses with constraints between them to fix
 168 the geometry and create the rest of the atoms as simple virtual sites
 169 type (see sec. :ref:`virtualsites`) from these three. Each of the
 170 aromatic residues require a different approach:
 171
 172 -  *Phenylalanine:* C\ :math:`_\gamma`,
 173    C\ :math:`_{{\epsilon}1}`, and
 174    C\ :math:`_{{\epsilon}2}` are kept as normal atoms,
 175    but with each a mass of one third the total mass of the phenyl group.
 176    See :numref:`Fig. %s A<fig-vsitehydro>`.
 177
 178 -  *Tyrosine:* The ring is treated identically to the phenylalanine
 179    ring. Additionally, constraints are defined between
 180    C\ :math:`_{{\epsilon}1}`,
 181    C\ :math:`_{{\epsilon}2}`, and
 182    O\ :math:`_{\eta}`. The original improper dihedral
 183    angles will keep both triangles (one for the ring and one with
 184    O\ :math:`_{\eta}`) in a plane, but due to the larger
 185    moments of inertia this construction will be much more stable. The
 186    bond-angle in the hydroxyl group will be constrained by a constraint
 187    between C\ :math:`_\gamma` and
 188    H\ :math:`_{\eta}`. **Note** that the hydrogen is not
 189    treated as a virtual site. See
 190    :numref:`Fig. %s B<fig-vsitehydro>`.
 191
 192 -  *Tryptophan:* C\ :math:`_\beta` is kept as a normal
 193    atom and two dummy masses are created at the center of mass of each
 194    of the rings, each with a mass equal to the total mass of the
 195    respective ring (C\ :math:`_{{\delta}2}` and
 196    C\ :math:`_{{\epsilon}2}` are each counted half for
 197    each ring). This keeps the overall center of mass and the moment of
 198    inertia almost (but not quite) equal to what it was. See
 199    :numref:`Fig. %s C<fig-vsitehydro>`.
 200
 201 -  *Histidine:* C\ :math:`_\gamma`,
 202    C\ :math:`_{{\epsilon}1}` and
 203    N\ :math:`_{{\epsilon}2}` are kept as normal atoms,
 204    but with masses redistributed such that the center of mass of the
 205    ring is preserved. See :numref:`Fig. %s D<fig-vsitehydro>`.