Fixed manual for dihedral restraints.

[gromacs.git] / docs / reference-manual / functions / restraints.rst

Restraints
----------

Special potentials are used for imposing restraints on the motion of the
system, either to avoid disastrous deviations, or to include knowledge
from experimental data. In either case they are not really part of the
force field and the reliability of the parameters is not important. The
potential forms, as implemented in |Gromacs|, are mentioned just for the
sake of completeness. Restraints and constraints refer to quite
different algorithms in |Gromacs|.

.. _positionrestraint:

Position restraints
~~~~~~~~~~~~~~~~~~~

These are used to restrain particles to fixed reference positions
:math:`\mathbf{R}_i`. They can be used during
equilibration in order to avoid drastic rearrangements of critical parts
(*e.g.* to restrain motion in a protein that is subjected to large
solvent forces when the solvent is not yet equilibrated). Another
application is the restraining of particles in a shell around a region
that is simulated in detail, while the shell is only approximated
because it lacks proper interaction from missing particles outside the
shell. Restraining will then maintain the integrity of the inner part.
For spherical shells, it is a wise procedure to make the force constant
depend on the radius, increasing from zero at the inner boundary to a
large value at the outer boundary. This feature has not, however, been
implemented in |Gromacs|.

The following form is used:

.. math:: V_{pr}(\mathbf{r}_i) = {\frac{1}{2}}k_{pr}|\mathbf{r}_i-\mathbf{R}_i|^2
          :label: eqnposrestform

The potential is plotted in :numref:`Fig. %s <fig-positionrestraint>`.

.. _fig-positionrestraint:

.. figure:: plots/f-pr.*
   :width: 8.00000cm

   Position restraint potential.

The potential form can be rewritten without loss of generality as:

.. math:: V_{pr}(\mathbf{r}_i) = {\frac{1}{2}} \left[ k_{pr}^x (x_i-X_i)^2 ~{\hat{\bf x}} + k_{pr}^y (y_i-Y_i)^2 ~{\hat{\bf y}} + k_{pr}^z (z_i-Z_i)^2 ~{\hat{\bf z}}\right]
         :label: eqnposrestgeneral

Now the forces are:

.. math:: \begin{array}{rcl}
          F_i^x &=& -k_{pr}^x~(x_i - X_i) \\
          F_i^y &=& -k_{pr}^y~(y_i - Y_i) \\
          F_i^z &=& -k_{pr}^z~(z_i - Z_i)
          \end{array}
          :label: eqnposrestforce

Using three different force constants the position restraints can be
turned on or off in each spatial dimension; this means that atoms can be
harmonically restrained to a plane or a line. Position restraints are
applied to a special fixed list of atoms. Such a list is usually
generated by the :ref:`pdb2gmx <gmx pdb2gmx>` program.

Flat-bottomed position restraints
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Flat-bottomed position restraints can be used to restrain particles to
part of the simulation volume. No force acts on the restrained particle
within the flat-bottomed region of the potential, however a harmonic
force acts to move the particle to the flat-bottomed region if it is
outside it. It is possible to apply normal and flat-bottomed position
restraints on the same particle (however, only with the same reference
position :math:`\mathbf{R}_i`). The following general
potential is used (:numref:`Figure %s <fig-fbposres>` A):

.. math:: V_\mathrm{fb}(\mathbf{r}_i) = \frac{1}{2}k_\mathrm{fb} [d_g(\mathbf{r}_i;\mathbf{R}_i) - r_\mathrm{fb}]^2\,H[d_g(\mathbf{r}_i;\mathbf{R}_i) - r_\mathrm{fb}],
          :label: eqnflatbottomposrest

where :math:`\mathbf{R}_i` is the reference position,
:math:`r_\mathrm{fb}` is the distance from the center with a flat
potential, :math:`k_\mathrm{fb}` the force constant, and :math:`H` is
the Heaviside step function. The distance
:math:`d_g(\mathbf{r}_i;\mathbf{R}_i)` from
the reference position depends on the geometry :math:`g` of the
flat-bottomed potential.

.. _fig-fbposres:

.. figure:: plots/fbposres.*
   :width: 10.00000cm

   Flat-bottomed position restraint potential. (A) Not inverted, (B)
   inverted.

| The following geometries for the flat-bottomed potential are
  supported:

| **Sphere** (:math:`g =1`): The
  particle is kept in a sphere of given radius. The force acts towards
  the center of the sphere. The following distance calculation is used:

  .. math:: d_g(\mathbf{r}_i;\mathbf{R}_i) = | \mathbf{r}_i-\mathbf{R}_i |
            :label: eqnfbsphereposrest

| **Cylinder** (:math:`g=6,7,8`): The particle is kept in a cylinder of
  given radius parallel to the :math:`x` (:math:`g=6`), :math:`y`
  (:math:`g=7`), or :math:`z`-axis (:math:`g=8`). For backwards
  compatibility, setting :math:`g=2` is mapped to :math:`g=8` in the
  code so that old :ref:`tpr` files and topologies work. The
  force from the flat-bottomed potential acts towards the axis of the
  cylinder. The component of the force parallel to the cylinder axis is
  zero. For a cylinder aligned along the :math:`z`-axis:

  .. math:: d_g(\mathbf{r}_i;\mathbf{R}_i) = \sqrt{ (x_i-X_i)^2 + (y_i - Y_i)^2 }
            :label: eqnfbcylinderposrest

| **Layer** (:math:`g=3,4,5`): The particle is kept in a layer defined
  by the thickness and the normal of the layer. The layer normal can be
  parallel to the :math:`x`, :math:`y`, or :math:`z`-axis. The force
  acts parallel to the layer normal.

  .. math:: d_g(\mathbf{r}_i;\mathbf{R}_i) = |x_i-X_i|, \;\;\;\mbox{or}\;\;\; 
            d_g(\mathbf{r}_i;\mathbf{R}_i) = |y_i-Y_i|, \;\;\;\mbox{or}\;\;\; 
            d_g(\mathbf{r}_i;\mathbf{R}_i) = |z_i-Z_i|.
            :label: eqnfblayerposrest

It is possible to apply multiple independent flat-bottomed position
restraints of different geometry on one particle. For example, applying
a cylinder and a layer in :math:`z` keeps a particle within a disk.
Applying three layers in :math:`x`, :math:`y`, and :math:`z` keeps the
particle within a cuboid.

In addition, it is possible to invert the restrained region with the
unrestrained region, leading to a potential that acts to keep the
particle *outside* of the volume defined by
:math:`\mathbf{R}_i`, :math:`g`, and
:math:`r_\mathrm{fb}`. That feature is switched on by defining a
negative :math:`r_\mathrm{fb}` in the topology. The following potential
is used (:numref:`Figure %s <fig-fbposres>` B):

.. math:: V_\mathrm{fb}^{\mathrm{inv}}(\mathbf{r}_i) = \frac{1}{2}k_\mathrm{fb}
          [d_g(\mathbf{r}_i;\mathbf{R}_i) - | r_\mathrm{fb} | ]^2\,
          H[ -(d_g(\mathbf{r}_i;\mathbf{R}_i) - | r_\mathrm{fb} | )].
          :label: eqninvertrest

Angle restraints
~~~~~~~~~~~~~~~~

These are used to restrain the angle between two pairs of particles or
between one pair of particles and the :math:`z`-axis. The functional
form is similar to that of a proper dihedral. For two pairs of atoms:

.. math:: V_{ar}(\mathbf{r}_i,\mathbf{r}_j,\mathbf{r}_k,\mathbf{r}_l)
                  = k_{ar}(1 - \cos(n (\theta - \theta_0))
                  )
          ,~~~~\mbox{where}~~
          \theta = \arccos\left(\frac{\mathbf{r}_j -\mathbf{r}_i}{\|\mathbf{r}_j -\mathbf{r}_i\|}
          \cdot \frac{\mathbf{r}_l -\mathbf{r}_k}{\|\mathbf{r}_l -\mathbf{r}_k\|} \right)
          :label: eqnanglerest

For one pair of atoms and the :math:`z`-axis:

.. math:: V_{ar}(\mathbf{r}_i,\mathbf{r}_j) = k_{ar}(1 - \cos(n (\theta - \theta_0))
                  )
          ,~~~~\mbox{where}~~
          \theta = \arccos\left(\frac{\mathbf{r}_j -\mathbf{r}_i}{\|\mathbf{r}_j -\mathbf{r}_i\|}
          \cdot \left( \begin{array}{c} 0 \\ 0 \\ 1 \\ \end{array} \right) \right)
          :label: eqnanglerestzaxis

A multiplicity (:math:`n`) of 2 is useful when you do not want to
distinguish between parallel and anti-parallel vectors. The equilibrium
angle :math:`\theta` should be between 0 and 180 degrees for
multiplicity 1 and between 0 and 90 degrees for multiplicity 2.

.. _dihedralrestraint:

Dihedral restraints
~~~~~~~~~~~~~~~~~~~

These are used to restrain the dihedral angle :math:`\phi` defined by
four particles as in an improper dihedral (sec. :ref:`imp`) but with a
slightly modified potential. Using:

.. math:: \phi' = \left(\phi-\phi_0\right) ~{\rm MOD}~ 2\pi
          :label: eqndphi

where :math:`\phi_0` is the reference angle, the potential is defined
as:

.. math:: V_{dihr}(\phi') ~=~ \left\{
          \begin{array}{lcllll}
          {\frac{1}{2}}k_{dihr}(\phi'-\Delta\phi)^2      
                          &\mbox{for}&     \|\phi'\| & >   & \Delta\phi       \\[1.5ex]
          0               &\mbox{for}&     \|\phi'\| & \le & \Delta\phi       \\[1.5ex]
          \end{array}\right.
          :label: eqndihre

where :math:`\Delta\phi` is a user defined angle and :math:`k_{dihr}`
is the force constant. **Note** that in the input in topology files,
angles are given in degrees and force constants in
kJ/mol/rad\ :math:`^2`.

.. _distancerestraint:

Distance restraints
~~~~~~~~~~~~~~~~~~~

Distance restraints add a penalty to the potential when the distance
between specified pairs of atoms exceeds a threshold value. They are
normally used to impose experimental restraints from, for instance,
experiments in nuclear magnetic resonance (NMR), on the motion of the
system. Thus, MD can be used for structure refinement using NMR data. In
|Gromacs| there are three ways to impose restraints on pairs of atoms:

-  Simple harmonic restraints: use ``[ bonds ]`` type 6 (see sec. :ref:`excl`).

-  Piecewise linear/harmonic restraints: ``[ bonds ]`` type
   10.

-  Complex NMR distance restraints, optionally with pair, time and/or
   ensemble averaging.

The last two options will be detailed now.

The potential form for distance restraints is quadratic below a
specified lower bound and between two specified upper bounds, and linear
beyond the largest bound (see :numref:`Fig. %s <fig-dist>`).

.. math:: V_{dr}(r_{ij}) ~=~ \left\{
          \begin{array}{lcllllll}
          {\frac{1}{2}}k_{dr}(r_{ij}-r_0)^2      
                          &\mbox{for}&     &     & r_{ij} & < & r_0       \\[1.5ex]
          0               &\mbox{for}& r_0 & \le & r_{ij} & < & r_1       \\[1.5ex]
          {\frac{1}{2}}k_{dr}(r_{ij}-r_1)^2      
                          &\mbox{for}& r_1 & \le & r_{ij} & < & r_2       \\[1.5ex]
          {\frac{1}{2}}k_{dr}(r_2-r_1)(2r_{ij}-r_2-r_1)  
                          &\mbox{for}& r_2 & \le & r_{ij} &   &
          \end{array}\right.
          :label: eqndisre

.. _fig-dist:

.. figure:: plots/f-dr.*
   :width: 8.00000cm

   Distance Restraint potential.

The forces are

.. math:: \mathbf{F}_i~=~ \left\{
          \begin{array}{lcllllll}
          -k_{dr}(r_{ij}-r_0)\frac{\mathbf{r}_ij}{r_{ij}} 
                          &\mbox{for}&     &     & r_{ij} & < & r_0       \\[1.5ex]
          0               &\mbox{for}& r_0 & \le & r_{ij} & < & r_1       \\[1.5ex]
          -k_{dr}(r_{ij}-r_1)\frac{\mathbf{r}_ij}{r_{ij}} 
                          &\mbox{for}& r_1 & \le & r_{ij} & < & r_2       \\[1.5ex]
          -k_{dr}(r_2-r_1)\frac{\mathbf{r}_ij}{r_{ij}}    
                          &\mbox{for}& r_2 & \le & r_{ij} &   &
          \end{array} \right.
          :label: eqndisreforce

For restraints not derived from NMR data, this functionality will
usually suffice and a section of ``[ bonds ]`` type 10 can be used to apply individual
restraints between pairs of atoms, see :ref:`topfile`. For applying
restraints derived from NMR measurements, more complex functionality
might be required, which is provided through the ``[ distance_restraints ]`` section and is
described below.

Time averaging
^^^^^^^^^^^^^^

Distance restraints based on instantaneous distances can potentially
reduce the fluctuations in a molecule significantly. This problem can be
overcome by restraining to a *time averaged*
distance \ :ref:`91 <refTorda89>`. The forces with time averaging are:

.. math:: \mathbf{F}_i~=~ \left\{
          \begin{array}{lcllllll}
          -k^a_{dr}(\bar{r}_{ij}-r_0)\frac{\mathbf{r}_ij}{r_{ij}}   
                          &\mbox{for}&     &     & \bar{r}_{ij} & < & r_0 \\[1.5ex]
          0               &\mbox{for}& r_0 & \le & \bar{r}_{ij} & < & r_1 \\[1.5ex]
          -k^a_{dr}(\bar{r}_{ij}-r_1)\frac{\mathbf{r}_ij}{r_{ij}}   
                          &\mbox{for}& r_1 & \le & \bar{r}_{ij} & < & r_2 \\[1.5ex]
          -k^a_{dr}(r_2-r_1)\frac{\mathbf{r}_ij}{r_{ij}}    
                          &\mbox{for}& r_2 & \le & \bar{r}_{ij} &   &
          \end{array} \right.
          :label: eqntimeaveragerest

where :math:`\bar{r}_{ij}` is given by an exponential running average
with decay time :math:`\tau`:

.. math:: \bar{r}_{ij} ~=~ < r_{ij}^{-3} >^{-1/3}
          :label: eqnrav

The force constant :math:`k^a_{dr}` is switched on slowly to compensate
for the lack of history at the beginning of the simulation:

.. math:: k^a_{dr} = k_{dr} \left(1-\exp\left(-\frac{t}{\tau}\right)\right)
          :label: eqnforceconstantswitch

Because of the time averaging, we can no longer speak of a distance
restraint potential.

This way an atom can satisfy two incompatible distance restraints *on
average* by moving between two positions. An example would be an amino
acid side-chain that is rotating around its :math:`\chi` dihedral angle,
thereby coming close to various other groups. Such a mobile side chain
can give rise to multiple NOEs that can not be fulfilled by a single
structure.

The computation of the time averaged distance in the
:ref:`mdrun <gmx mdrun>` program is done in the following fashion:

.. math:: \begin{array}{rcl}
          \overline{r^{-3}}_{ij}(0)       &=& r_{ij}(0)^{-3}      \\
          \overline{r^{-3}}_{ij}(t)       &=& \overline{r^{-3}}_{ij}(t-\Delta t)~\exp{\left(-\frac{\Delta t}{\tau}\right)} + r_{ij}(t)^{-3}\left[1-\exp{\left(-\frac{\Delta t}{\tau}\right)}\right]
          \end{array}
          :label: eqnravdisre

When a pair is within the bounds, it can still feel a force because the
time averaged distance can still be beyond a bound. To prevent the
protons from being pulled too close together, a mixed approach can be
used. In this approach, the penalty is zero when the instantaneous
distance is within the bounds, otherwise the violation is the square
root of the product of the instantaneous violation and the time averaged
violation:

.. math:: \mathbf{F}_i~=~ \left\{
          \begin{array}{lclll}
          k^a_{dr}\sqrt{(r_{ij}-r_0)(\bar{r}_{ij}-r_0)}\frac{\mathbf{r}_ij}{r_{ij}}   
              & \mbox{for} & r_{ij} < r_0 & \mbox{and} & \bar{r}_{ij} < r_0 \\[1.5ex]
          -k^a _{dr} \,
            \mbox{min}\left(\sqrt{(r_{ij}-r_1)(\bar{r}_{ij}-r_1)},r_2-r_1\right)
            \frac{\mathbf{r}_ij}{r_{ij}}   
              & \mbox{for} & r_{ij} > r_1 & \mbox{and} & \bar{r}_{ij} > r_1 \\[1.5ex]
          0               &\mbox{otherwise}
          \end{array} \right.
          :label: eqntimeaverageviolation

Averaging over multiple pairs
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Sometimes it is unclear from experimental data which atom pair gives
rise to a single NOE, in other occasions it can be obvious that more
than one pair contributes due to the symmetry of the system, *e.g.* a
methyl group with three protons. For such a group, it is not possible to
distinguish between the protons, therefore they should all be taken into
account when calculating the distance between this methyl group and
another proton (or group of protons). Due to the physical nature of
magnetic resonance, the intensity of the NOE signal is inversely
proportional to the sixth power of the inter-atomic distance. Thus, when
combining atom pairs, a fixed list of :math:`N` restraints may be taken
together, where the apparent “distance” is given by:

.. math:: r_N(t) = \left [\sum_{n=1}^{N} \bar{r}_{n}(t)^{-6} \right]^{-1/6}
          :label: eqnrsix

where we use :math:`r_{ij}` or :eq:`eqn. %s <eqnrav>` for the
:math:`\bar{r}_{n}`. The :math:`r_N` of the instantaneous and
time-averaged distances can be combined to do a mixed restraining, as
indicated above. As more pairs of protons contribute to the same NOE
signal, the intensity will increase, and the summed “distance” will be
shorter than any of its components due to the reciprocal summation.

There are two options for distributing the forces over the atom pairs.
In the conservative option, the force is defined as the derivative of
the restraint potential with respect to the coordinates. This results in
a conservative potential when time averaging is not used. The force
distribution over the pairs is proportional to :math:`r^{-6}`. This
means that a close pair feels a much larger force than a distant pair,
which might lead to a molecule that is “too rigid.” The other option is
an equal force distribution. In this case each pair feels :math:`1/N` of
the derivative of the restraint potential with respect to :math:`r_N`.
The advantage of this method is that more conformations might be
sampled, but the non-conservative nature of the forces can lead to local
heating of the protons.

It is also possible to use *ensemble averaging* using multiple (protein)
molecules. In this case the bounds should be lowered as in:

.. math:: \begin{array}{rcl}
          r_1     &~=~&   r_1 * M^{-1/6}  \\
          r_2     &~=~&   r_2 * M^{-1/6}
          \end{array}
          :label: eqnrestforceensembleaverage

where :math:`M` is the number of molecules. The |Gromacs| preprocessor
:ref:`grompp <gmx grompp>` can do this automatically when the appropriate
option is given. The resulting “distance” is then used to calculate the
scalar force according to:

.. math:: \mathbf{F}_i~=~\left\{
          \begin{array}{rcl}
          ~& 0 \hspace{4cm}  & r_{N} < r_1         \\
           & k_{dr}(r_{N}-r_1)\frac{\mathbf{r}_ij}{r_{ij}} & r_1 \le r_{N} < r_2 \\
           & k_{dr}(r_2-r_1)\frac{\mathbf{r}_ij}{r_{ij}}    & r_{N} \ge r_2 
          \end{array} \right.
          :label: eqnrestscalarforce

where :math:`i` and :math:`j` denote the atoms of all the pairs that
contribute to the NOE signal.

Using distance restraints
^^^^^^^^^^^^^^^^^^^^^^^^^

A list of distance restrains based on NOE data can be added to a
molecule definition in your topology file, like in the following
example:

::

    [ distance_restraints ]
    ; ai   aj   type   index   type'      low     up1     up2     fac
    10     16      1       0       1      0.0     0.3     0.4     1.0
    10     28      1       1       1      0.0     0.3     0.4     1.0
    10     46      1       1       1      0.0     0.3     0.4     1.0
    16     22      1       2       1      0.0     0.3     0.4     2.5
    16     34      1       3       1      0.0     0.5     0.6     1.0

In this example a number of features can be found. In columns ai and aj
you find the atom numbers of the particles to be restrained. The type
column should always be 1. As explained in  :ref:`distancerestraint`,
multiple distances can contribute to a single NOE signal. In the
topology this can be set using the index column. In our example, the
restraints 10-28 and 10-46 both have index 1, therefore they are treated
simultaneously. An extra requirement for treating restraints together is
that the restraints must be on successive lines, without any other
intervening restraint. The type’ column will usually be 1, but can be
set to 2 to obtain a distance restraint that will never be time- and
ensemble-averaged; this can be useful for restraining hydrogen bonds.
The columns ``low``, ``up1``, and
``up2`` hold the values of :math:`r_0`, :math:`r_1`, and
:math:`r_2` from  :eq:`eqn. %s <eqndisre>`. In some cases it
can be useful to have different force constants for some restraints;
this is controlled by the column ``fac``. The force constant
in the parameter file is multiplied by the value in the column
``fac`` for each restraint. Information for each restraint
is stored in the energy file and can be processed and plotted with
:ref:`gmx nmr`.

Orientation restraints
~~~~~~~~~~~~~~~~~~~~~~

This section describes how orientations between vectors, as measured in
certain NMR experiments, can be calculated and restrained in MD
simulations. The presented refinement methodology and a comparison of
results with and without time and ensemble averaging have been
published \ :ref:`92 <refHess2003>`.

Theory
^^^^^^

In an NMR experiment, orientations of vectors can be measured when a
molecule does not tumble completely isotropically in the solvent. Two
examples of such orientation measurements are residual dipolar couplings
(between two nuclei) or chemical shift anisotropies. An observable for a
vector :math:`\mathbf{r}_i` can be written as follows:

.. math:: \delta_i = \frac{2}{3} \mbox{tr}({{\mathbf S}}{{\mathbf D}}_i)
          :label: eqnorrestvector

where :math:`{{\mathbf S}}` is the dimensionless order tensor of the
molecule. The tensor :math:`{{\mathbf D}}_i` is given by:

.. math:: {{\mathbf D}}_i = \frac{c_i}{\|\mathbf{r}_i\|^\alpha} \left(
          \begin{array}{lll}
          3 x x - 1 & 3 x y     & 3 x z     \\
          3 x y     & 3 y y - 1 & 3 y z     \\
          3 x z     & 3 y z     & 3 z z - 1 \\
          \end{array} \right)
          :label: eqnorientdef

.. math:: \mbox{with:} \quad 
          x=\frac{r_{i,x}}{\|\mathbf{r}_i\|}, \quad
          y=\frac{r_{i,y}}{\|\mathbf{r}_i\|}, \quad 
          z=\frac{r_{i,z}}{\|\mathbf{r}_i\|}
          :label: eqnorientdef2

For a dipolar coupling :math:`\mathbf{r}_i` is the vector
connecting the two nuclei, :math:`\alpha=3` and the constant :math:`c_i`
is given by:

.. math:: c_i = \frac{\mu_0}{4\pi} \gamma_1^i \gamma_2^i \frac{\hbar}{4\pi}
          :label: eqnorrestconstant

where :math:`\gamma_1^i` and :math:`\gamma_2^i` are the gyromagnetic
ratios of the two nuclei.

The order tensor is symmetric and has trace zero. Using a rotation
matrix :math:`{\mathbf T}` it can be transformed into the following
form:

.. math:: {\mathbf T}^T {{\mathbf S}}{\mathbf T} = s \left( \begin{array}{ccc}
          -\frac{1}{2}(1-\eta) & 0                    & 0 \\
          0                    & -\frac{1}{2}(1+\eta) & 0 \\
          0                    & 0                    & 1
          \end{array} \right)
          :label: eqnorresttensor

where :math:`-1 \leq s \leq 1` and :math:`0 \leq \eta \leq 1`.
:math:`s` is called the order parameter and :math:`\eta` the asymmetry
of the order tensor :math:`{{\mathbf S}}`. When the molecule tumbles
isotropically in the solvent, :math:`s` is zero, and no orientational
effects can be observed because all :math:`\delta_i` are zero.

Calculating orientations in a simulation
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

For reasons which are explained below, the :math:`{{\mathbf D}}`
matrices are calculated which respect to a reference orientation of the
molecule. The orientation is defined by a rotation matrix
:math:`{{\mathbf R}}`, which is needed to least-squares fit the current
coordinates of a selected set of atoms onto a reference conformation.
The reference conformation is the starting conformation of the
simulation. In case of ensemble averaging, which will be treated later,
the structure is taken from the first subsystem. The calculated
:math:`{{\mathbf D}}_i^c` matrix is given by:

.. math:: {{\mathbf D}}_i^c(t) = {{\mathbf R}}(t) {{\mathbf D}}_i(t) {{\mathbf R}}^T(t)
          :label: eqnDrot

The calculated orientation for vector :math:`i` is given by:

.. math:: \delta^c_i(t) = \frac{2}{3} \mbox{tr}({{\mathbf S}}(t){{\mathbf D}}_i^c(t))
          :label: eqnDrotvector

The order tensor :math:`{{\mathbf S}}(t)` is usually unknown. A
reasonable choice for the order tensor is the tensor which minimizes the
(weighted) mean square difference between the calculated and the
observed orientations:

.. math:: MSD(t) = \left(\sum_{i=1}^N w_i\right)^{-1} \sum_{i=1}^N w_i (\delta_i^c (t) -\delta_i^{exp})^2
          :label: eqnSmsd

To properly combine different types of measurements, the unit of
:math:`w_i` should be such that all terms are dimensionless. This means
the unit of :math:`w_i` is the unit of :math:`\delta_i` to the power
:math:`-2`. **Note** that scaling all :math:`w_i` with a constant factor
does not influence the order tensor.

Time averaging
^^^^^^^^^^^^^^

Since the tensors :math:`{{\mathbf D}}_i` fluctuate rapidly in time,
much faster than can be observed in an experiment, they should be
averaged over time in the simulation. However, in a simulation the time
and the number of copies of a molecule are limited. Usually one can not
obtain a converged average of the :math:`{{\mathbf D}}_i` tensors over
all orientations of the molecule. If one assumes that the average
orientations of the :math:`\mathbf{r}_i` vectors within
the molecule converge much faster than the tumbling time of the
molecule, the tensor can be averaged in an axis system that rotates with
the molecule, as expressed by :eq:`equation %s <eqnDrot>`). The time-averaged
tensors are calculated using an exponentially decaying memory function:

.. math:: {{\mathbf D}}^a_i(t) = \frac{\displaystyle
          \int_{u=t_0}^t {{\mathbf D}}^c_i(u) \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
          }{\displaystyle
          \int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
          }
          :label: eqnorresttimeaverage

Assuming that the order tensor :math:`{{\mathbf S}}` fluctuates slower
than the :math:`{{\mathbf D}}_i`, the time-averaged orientation can be
calculated as:

.. math:: \delta_i^a(t) = \frac{2}{3} \mbox{tr}({{\mathbf S}}(t) {{\mathbf D}}_i^a(t))
          :label: eqnorresttimeaveorient

where the order tensor :math:`{{\mathbf S}}(t)` is calculated using
expression :eq:`%s <eqnSmsd>` with :math:`\delta_i^c(t)` replaced by
:math:`\delta_i^a(t)`.

Restraining
^^^^^^^^^^^

The simulated structure can be restrained by applying a force
proportional to the difference between the calculated and the
experimental orientations. When no time averaging is applied, a proper
potential can be defined as:

.. math:: V = \frac{1}{2} k \sum_{i=1}^N w_i (\delta_i^c (t) -\delta_i^{exp})^2
          :label: eqnorrestsimrest

where the unit of :math:`k` is the unit of energy. Thus the effective
force constant for restraint :math:`i` is :math:`k w_i`. The forces are
given by minus the gradient of :math:`V`. The force
:math:`\mathbf{F}\!_i` working on vector
:math:`\mathbf{r}_i` is:

.. math:: \begin{aligned}
          \mathbf{F}\!_i(t) 
          & = & - \frac{\mbox{d} V}{\mbox{d}\mathbf{r}_i} \\
          & = & -k w_i (\delta_i^c (t) -\delta_i^{exp}) \frac{\mbox{d} \delta_i (t)}{\mbox{d}\mathbf{r}_i} \\
          & = & -k w_i (\delta_i^c (t) -\delta_i^{exp})
          \frac{2 c_i}{\|\mathbf{r}\|^{2+\alpha}} \left(2 {{\mathbf R}}^T {{\mathbf S}}{{\mathbf R}}\mathbf{r}_i - \frac{2+\alpha}{\|\mathbf{r}\|^2} \mbox{tr}({{\mathbf R}}^T {{\mathbf S}}{{\mathbf R}}\mathbf{r}_i \mathbf{r}_i^T) \mathbf{r}_i \right)\end{aligned}
          :label: eqnorrestsimrestforce

Ensemble averaging
^^^^^^^^^^^^^^^^^^

Ensemble averaging can be applied by simulating a system of :math:`M`
subsystems that each contain an identical set of orientation restraints.
The systems only interact via the orientation restraint potential which
is defined as:

.. math:: V = M \frac{1}{2} k \sum_{i=1}^N w_i 
          \langle \delta_i^c (t) -\delta_i^{exp} \rangle^2
          :label: eqnorrestensembleave

The force on vector :math:`\mathbf{r}_{i,m}` in subsystem
:math:`m` is given by:

.. math:: \mathbf{F}\!_{i,m}(t) = - \frac{\mbox{d} V}{\mbox{d}\mathbf{r}_{i,m}} =
          -k w_i \langle \delta_i^c (t) -\delta_i^{exp} \rangle \frac{\mbox{d} \delta_{i,m}^c (t)}{\mbox{d}\mathbf{r}_{i,m}}
          :label: eqnorrestensaveforce 

Time averaging
^^^^^^^^^^^^^^

When using time averaging it is not possible to define a potential. We
can still define a quantity that gives a rough idea of the energy stored
in the restraints:

.. math:: V = M \frac{1}{2} k^a \sum_{i=1}^N w_i 
          \langle \delta_i^a (t) -\delta_i^{exp} \rangle^2
          :label: eqntimeavepot

The force constant :math:`k_a` is switched on slowly to compensate for
the lack of history at times close to :math:`t_0`. It is exactly
proportional to the amount of average that has been accumulated:

.. math:: k^a =
          k \, \frac{1}{\tau}\int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
          :label: eqntimeaveforceswitch

What really matters is the definition of the force. It is chosen to be
proportional to the square root of the product of the time-averaged and
the instantaneous deviation. Using only the time-averaged deviation
induces large oscillations. The force is given by:

.. math:: \mathbf{F}\!_{i,m}(t) =
          \left\{ \begin{array}{ll}
          0 & \quad \mbox{for} \quad a\, b \leq 0 \\
          \displaystyle
          k^a w_i \frac{a}{|a|} \sqrt{a\, b} \, \frac{\mbox{d} \delta_{i,m}^c (t)}{\mbox{d}\mathbf{r}_{i,m}}
          & \quad \mbox{for} \quad a\, b > 0 
          \end{array}
          \right.
          :label: eqntimeaveforce

.. math:: \begin{aligned}
          a &=& \langle \delta_i^a (t) -\delta_i^{exp} \rangle \\
          b &=& \langle \delta_i^c (t) -\delta_i^{exp} \rangle\end{aligned}
          :label: eqntimeaveforce2

Using orientation restraints
^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Orientation restraints can be added to a molecule definition in the
topology file in the section ``[ orientation_restraints ]``.
Here we give an example section containing five N-H residual dipolar
coupling restraints:

::

    [ orientation_restraints ]
    ; ai   aj  type  exp.  label  alpha    const.     obs.   weight
    ;                                Hz      nm^3       Hz    Hz^-2
      31   32     1     1      3      3     6.083    -6.73      1.0
      43   44     1     1      4      3     6.083    -7.87      1.0
      55   56     1     1      5      3     6.083    -7.13      1.0
      65   66     1     1      6      3     6.083    -2.57      1.0
      73   74     1     1      7      3     6.083    -2.10      1.0

The unit of the observable is Hz, but one can choose any other unit. In
columns ``ai`` and ``aj`` you find the atom numbers of the particles to be
restrained. The ``type`` column should always be 1. The ``exp.`` column denotes
the experiment number, starting at 1. For each experiment a separate
order tensor :math:`{{\mathbf S}}` is optimized. The label should be a
unique number larger than zero for each restraint. The ``alpha`` column
contains the power :math:`\alpha` that is used in
:eq:`equation %s <eqnorientdef>`) to calculate the orientation. The ``const.`` column
contains the constant :math:`c_i` used in the same equation. The
constant should have the unit of the observable times
nm\ :math:`^\alpha`. The column ``obs.`` contains the observable, in any
unit you like. The last column contains the weights :math:`w_i`; the
unit should be the inverse of the square of the unit of the observable.

Some parameters for orientation restraints can be specified in the
:ref:`grompp <gmx grompp>` :ref:`mdp` file, for a study of the effect of different
force constants and averaging times and ensemble averaging see \ :ref:`92 <refHess2003>`.
Information for each restraint is stored in the energy
file and can be processed and plotted with :ref:`gmx nmr`.

.. raw:: latex

    \clearpage