4 Special potentials are used for imposing restraints on the motion of the
5 system, either to avoid disastrous deviations, or to include knowledge
6 from experimental data. In either case they are not really part of the
7 force field and the reliability of the parameters is not important. The
8 potential forms, as implemented in |Gromacs|, are mentioned just for the
9 sake of completeness. Restraints and constraints refer to quite
10 different algorithms in |Gromacs|.
12 .. _positionrestraint:
17 These are used to restrain particles to fixed reference positions
18 :math:`\mathbf{R}_i`. They can be used during
19 equilibration in order to avoid drastic rearrangements of critical parts
20 (*e.g.* to restrain motion in a protein that is subjected to large
21 solvent forces when the solvent is not yet equilibrated). Another
22 application is the restraining of particles in a shell around a region
23 that is simulated in detail, while the shell is only approximated
24 because it lacks proper interaction from missing particles outside the
25 shell. Restraining will then maintain the integrity of the inner part.
26 For spherical shells, it is a wise procedure to make the force constant
27 depend on the radius, increasing from zero at the inner boundary to a
28 large value at the outer boundary. This feature has not, however, been
29 implemented in |Gromacs|.
31 The following form is used:
33 .. math:: V_{pr}(\mathbf{r}_i) = {\frac{1}{2}}k_{pr}|\mathbf{r}_i-\mathbf{R}_i|^2
35 The potential is plotted in :numref:`Fig. %s <fig-positionrestraint>`.
37 .. _fig-positionrestraint:
39 .. figure:: plots/f-pr.*
42 Position restraint potential.
44 The potential form can be rewritten without loss of generality as:
46 .. 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]
53 F_i^x &=& -k_{pr}^x~(x_i - X_i) \\
54 F_i^y &=& -k_{pr}^y~(y_i - Y_i) \\
55 F_i^z &=& -k_{pr}^z~(z_i - Z_i)
58 Using three different force constants the position restraints can be
59 turned on or off in each spatial dimension; this means that atoms can be
60 harmonically restrained to a plane or a line. Position restraints are
61 applied to a special fixed list of atoms. Such a list is usually
62 generated by the :ref:`pdb2gmx <gmx pdb2gmx>` program.
64 Flat-bottomed position restraints
65 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
67 Flat-bottomed position restraints can be used to restrain particles to
68 part of the simulation volume. No force acts on the restrained particle
69 within the flat-bottomed region of the potential, however a harmonic
70 force acts to move the particle to the flat-bottomed region if it is
71 outside it. It is possible to apply normal and flat-bottomed position
72 restraints on the same particle (however, only with the same reference
73 position :math:`\mathbf{R}_i`). The following general
74 potential is used (:numref:`Figure %s <fig-fbposres>` A):
76 .. 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}],
78 where :math:`\mathbf{R}_i` is the reference position,
79 :math:`r_\mathrm{fb}` is the distance from the center with a flat
80 potential, :math:`k_\mathrm{fb}` the force constant, and :math:`H` is
81 the Heaviside step function. The distance
82 :math:`d_g(\mathbf{r}_i;\mathbf{R}_i)` from
83 the reference position depends on the geometry :math:`g` of the
84 flat-bottomed potential.
88 .. figure:: plots/fbposres.*
91 Flat-bottomed position restraint potential. (A) Not inverted, (B)
94 | The following geometries for the flat-bottomed potential are
97 | **Sphere** (:math:`g =1`): The
98 particle is kept in a sphere of given radius. The force acts towards
99 the center of the sphere. The following distance calculation is used:
101 .. math:: d_g(\mathbf{r}_i;\mathbf{R}_i) = | \mathbf{r}_i-\mathbf{R}_i |
103 | **Cylinder** (:math:`g=6,7,8`): The particle is kept in a cylinder of
104 given radius parallel to the :math:`x` (:math:`g=6`), :math:`y`
105 (:math:`g=7`), or :math:`z`-axis (:math:`g=8`). For backwards
106 compatibility, setting :math:`g=2` is mapped to :math:`g=8` in the
107 code so that old :ref:`tpr` files and topologies work. The
108 force from the flat-bottomed potential acts towards the axis of the
109 cylinder. The component of the force parallel to the cylinder axis is
110 zero. For a cylinder aligned along the :math:`z`-axis:
112 .. math:: d_g(\mathbf{r}_i;\mathbf{R}_i) = \sqrt{ (x_i-X_i)^2 + (y_i - Y_i)^2 }
114 | **Layer** (:math:`g=3,4,5`): The particle is kept in a layer defined
115 by the thickness and the normal of the layer. The layer normal can be
116 parallel to the :math:`x`, :math:`y`, or :math:`z`-axis. The force
117 acts parallel to the layer normal.
121 d_g(\mathbf{r}_i;\mathbf{R}_i) = |x_i-X_i|, \;\;\;\mbox{or}\;\;\;
122 d_g(\mathbf{r}_i;\mathbf{R}_i) = |y_i-Y_i|, \;\;\;\mbox{or}\;\;\;
123 d_g(\mathbf{r}_i;\mathbf{R}_i) = |z_i-Z_i|.
125 It is possible to apply multiple independent flat-bottomed position
126 restraints of different geometry on one particle. For example, applying
127 a cylinder and a layer in :math:`z` keeps a particle within a disk.
128 Applying three layers in :math:`x`, :math:`y`, and :math:`z` keeps the
129 particle within a cuboid.
131 In addition, it is possible to invert the restrained region with the
132 unrestrained region, leading to a potential that acts to keep the
133 particle *outside* of the volume defined by
134 :math:`\mathbf{R}_i`, :math:`g`, and
135 :math:`r_\mathrm{fb}`. That feature is switched on by defining a
136 negative :math:`r_\mathrm{fb}` in the topology. The following potential
137 is used (:numref:`Figure %s <fig-fbposres>` B):
141 V_\mathrm{fb}^{\mathrm{inv}}(\mathbf{r}_i) = \frac{1}{2}k_\mathrm{fb}
142 [d_g(\mathbf{r}_i;\mathbf{R}_i) - | r_\mathrm{fb} | ]^2\,
143 H[ -(d_g(\mathbf{r}_i;\mathbf{R}_i) - | r_\mathrm{fb} | )].
148 These are used to restrain the angle between two pairs of particles or
149 between one pair of particles and the :math:`z`-axis. The functional
150 form is similar to that of a proper dihedral. For two pairs of atoms:
154 V_{ar}(\mathbf{r}_i,\mathbf{r}_j,\mathbf{r}_k,\mathbf{r}_l)
155 = k_{ar}(1 - \cos(n (\theta - \theta_0))
158 \theta = \arccos\left(\frac{\mathbf{r}_j -\mathbf{r}_i}{\|\mathbf{r}_j -\mathbf{r}_i\|}
159 \cdot \frac{\mathbf{r}_l -\mathbf{r}_k}{\|\mathbf{r}_l -\mathbf{r}_k\|} \right)
161 For one pair of atoms and the :math:`z`-axis:
165 V_{ar}(\mathbf{r}_i,\mathbf{r}_j) = k_{ar}(1 - \cos(n (\theta - \theta_0))
168 \theta = \arccos\left(\frac{\mathbf{r}_j -\mathbf{r}_i}{\|\mathbf{r}_j -\mathbf{r}_i\|}
169 \cdot \left( \begin{array}{c} 0 \\ 0 \\ 1 \\ \end{array} \right) \right)
171 A multiplicity (:math:`n`) of 2 is useful when you do not want to
172 distinguish between parallel and anti-parallel vectors. The equilibrium
173 angle :math:`\theta` should be between 0 and 180 degrees for
174 multiplicity 1 and between 0 and 90 degrees for multiplicity 2.
176 .. _dihedralrestraint:
181 These are used to restrain the dihedral angle :math:`\phi` defined by
182 four particles as in an improper dihedral (sec. :ref:`imp`) but with a
183 slightly modified potential. Using:
185 .. math:: \phi' = \left(\phi-\phi_0\right) ~{\rm MOD}~ 2\pi
188 where :math:`\phi_0` is the reference angle, the potential is defined
191 .. math:: V_{dihr}(\phi') ~=~ \left\{
192 \begin{array}{lcllll}
193 {\frac{1}{2}}k_{dihr}(\phi'-\phi_0-\Delta\phi)^2
194 &\mbox{for}& \phi' & > & \Delta\phi \\[1.5ex]
195 0 &\mbox{for}& \phi' & \le & \Delta\phi \\[1.5ex]
199 where :math:`\Delta\phi` is a user defined angle and :math:`k_{dihr}`
200 is the force constant. **Note** that in the input in topology files,
201 angles are given in degrees and force constants in
202 kJ/mol/rad\ :math:`^2`.
204 .. _distancerestraint:
209 Distance restraints add a penalty to the potential when the distance
210 between specified pairs of atoms exceeds a threshold value. They are
211 normally used to impose experimental restraints from, for instance,
212 experiments in nuclear magnetic resonance (NMR), on the motion of the
213 system. Thus, MD can be used for structure refinement using NMR data. In
214 |Gromacs| there are three ways to impose restraints on pairs of atoms:
216 - Simple harmonic restraints: use ``[ bonds ]`` type 6 (see sec. :ref:`excl`).
218 - Piecewise linear/harmonic restraints: ``[ bonds ]`` type
221 - Complex NMR distance restraints, optionally with pair, time and/or
224 The last two options will be detailed now.
226 The potential form for distance restraints is quadratic below a
227 specified lower bound and between two specified upper bounds, and linear
228 beyond the largest bound (see :numref:`Fig. %s <fig-dist>`).
230 .. math:: V_{dr}(r_{ij}) ~=~ \left\{
231 \begin{array}{lcllllll}
232 {\frac{1}{2}}k_{dr}(r_{ij}-r_0)^2
233 &\mbox{for}& & & r_{ij} & < & r_0 \\[1.5ex]
234 0 &\mbox{for}& r_0 & \le & r_{ij} & < & r_1 \\[1.5ex]
235 {\frac{1}{2}}k_{dr}(r_{ij}-r_1)^2
236 &\mbox{for}& r_1 & \le & r_{ij} & < & r_2 \\[1.5ex]
237 {\frac{1}{2}}k_{dr}(r_2-r_1)(2r_{ij}-r_2-r_1)
238 &\mbox{for}& r_2 & \le & r_{ij} & &
244 .. figure:: plots/f-dr.*
247 Distance Restraint potential.
253 \mathbf{F}_i~=~ \left\{
254 \begin{array}{lcllllll}
255 -k_{dr}(r_{ij}-r_0)\frac{\mathbf{r}_ij}{r_{ij}}
256 &\mbox{for}& & & r_{ij} & < & r_0 \\[1.5ex]
257 0 &\mbox{for}& r_0 & \le & r_{ij} & < & r_1 \\[1.5ex]
258 -k_{dr}(r_{ij}-r_1)\frac{\mathbf{r}_ij}{r_{ij}}
259 &\mbox{for}& r_1 & \le & r_{ij} & < & r_2 \\[1.5ex]
260 -k_{dr}(r_2-r_1)\frac{\mathbf{r}_ij}{r_{ij}}
261 &\mbox{for}& r_2 & \le & r_{ij} & &
264 For restraints not derived from NMR data, this functionality will
265 usually suffice and a section of ``[ bonds ]`` type 10 can be used to apply individual
266 restraints between pairs of atoms, see :ref:`topfile`. For applying
267 restraints derived from NMR measurements, more complex functionality
268 might be required, which is provided through the ``[ distance_restraints ]`` section and is
274 Distance restraints based on instantaneous distances can potentially
275 reduce the fluctuations in a molecule significantly. This problem can be
276 overcome by restraining to a *time averaged*
277 distance \ :ref:`91 <refTorda89>`. The forces with time averaging are:
281 \mathbf{F}_i~=~ \left\{
282 \begin{array}{lcllllll}
283 -k^a_{dr}(\bar{r}_{ij}-r_0)\frac{\mathbf{r}_ij}{r_{ij}}
284 &\mbox{for}& & & \bar{r}_{ij} & < & r_0 \\[1.5ex]
285 0 &\mbox{for}& r_0 & \le & \bar{r}_{ij} & < & r_1 \\[1.5ex]
286 -k^a_{dr}(\bar{r}_{ij}-r_1)\frac{\mathbf{r}_ij}{r_{ij}}
287 &\mbox{for}& r_1 & \le & \bar{r}_{ij} & < & r_2 \\[1.5ex]
288 -k^a_{dr}(r_2-r_1)\frac{\mathbf{r}_ij}{r_{ij}}
289 &\mbox{for}& r_2 & \le & \bar{r}_{ij} & &
292 where :math:`\bar{r}_{ij}` is given by an exponential running average
293 with decay time :math:`\tau`:
295 .. math:: \bar{r}_{ij} ~=~ < r_{ij}^{-3} >^{-1/3}
298 The force constant :math:`k^a_{dr}` is switched on slowly to compensate
299 for the lack of history at the beginning of the simulation:
301 .. math:: k^a_{dr} = k_{dr} \left(1-\exp\left(-\frac{t}{\tau}\right)\right)
303 Because of the time averaging, we can no longer speak of a distance
306 This way an atom can satisfy two incompatible distance restraints *on
307 average* by moving between two positions. An example would be an amino
308 acid side-chain that is rotating around its :math:`\chi` dihedral angle,
309 thereby coming close to various other groups. Such a mobile side chain
310 can give rise to multiple NOEs that can not be fulfilled by a single
313 The computation of the time averaged distance in the
314 :ref:`mdrun <gmx mdrun>` program is done in the following fashion:
316 .. math:: \begin{array}{rcl}
317 \overline{r^{-3}}_{ij}(0) &=& r_{ij}(0)^{-3} \\
318 \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]
322 When a pair is within the bounds, it can still feel a force because the
323 time averaged distance can still be beyond a bound. To prevent the
324 protons from being pulled too close together, a mixed approach can be
325 used. In this approach, the penalty is zero when the instantaneous
326 distance is within the bounds, otherwise the violation is the square
327 root of the product of the instantaneous violation and the time averaged
332 \mathbf{F}_i~=~ \left\{
334 k^a_{dr}\sqrt{(r_{ij}-r_0)(\bar{r}_{ij}-r_0)}\frac{\mathbf{r}_ij}{r_{ij}}
335 & \mbox{for} & r_{ij} < r_0 & \mbox{and} & \bar{r}_{ij} < r_0 \\[1.5ex]
337 \mbox{min}\left(\sqrt{(r_{ij}-r_1)(\bar{r}_{ij}-r_1)},r_2-r_1\right)
338 \frac{\mathbf{r}_ij}{r_{ij}}
339 & \mbox{for} & r_{ij} > r_1 & \mbox{and} & \bar{r}_{ij} > r_1 \\[1.5ex]
343 Averaging over multiple pairs
344 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
346 Sometimes it is unclear from experimental data which atom pair gives
347 rise to a single NOE, in other occasions it can be obvious that more
348 than one pair contributes due to the symmetry of the system, *e.g.* a
349 methyl group with three protons. For such a group, it is not possible to
350 distinguish between the protons, therefore they should all be taken into
351 account when calculating the distance between this methyl group and
352 another proton (or group of protons). Due to the physical nature of
353 magnetic resonance, the intensity of the NOE signal is inversely
354 proportional to the sixth power of the inter-atomic distance. Thus, when
355 combining atom pairs, a fixed list of :math:`N` restraints may be taken
356 together, where the apparent “distance” is given by:
358 .. math:: r_N(t) = \left [\sum_{n=1}^{N} \bar{r}_{n}(t)^{-6} \right]^{-1/6}
361 where we use :math:`r_{ij}` or :eq:`eqn. %s <eqnrav>` for the
362 :math:`\bar{r}_{n}`. The :math:`r_N` of the instantaneous and
363 time-averaged distances can be combined to do a mixed restraining, as
364 indicated above. As more pairs of protons contribute to the same NOE
365 signal, the intensity will increase, and the summed “distance” will be
366 shorter than any of its components due to the reciprocal summation.
368 There are two options for distributing the forces over the atom pairs.
369 In the conservative option, the force is defined as the derivative of
370 the restraint potential with respect to the coordinates. This results in
371 a conservative potential when time averaging is not used. The force
372 distribution over the pairs is proportional to :math:`r^{-6}`. This
373 means that a close pair feels a much larger force than a distant pair,
374 which might lead to a molecule that is “too rigid.” The other option is
375 an equal force distribution. In this case each pair feels :math:`1/N` of
376 the derivative of the restraint potential with respect to :math:`r_N`.
377 The advantage of this method is that more conformations might be
378 sampled, but the non-conservative nature of the forces can lead to local
379 heating of the protons.
381 It is also possible to use *ensemble averaging* using multiple (protein)
382 molecules. In this case the bounds should be lowered as in:
387 r_1 &~=~& r_1 * M^{-1/6} \\
388 r_2 &~=~& r_2 * M^{-1/6}
391 where :math:`M` is the number of molecules. The |Gromacs| preprocessor
392 :ref:`grompp <gmx grompp>` can do this automatically when the appropriate
393 option is given. The resulting “distance” is then used to calculate the
394 scalar force according to:
398 \mathbf{F}_i~=~\left\{
400 ~& 0 \hspace{4cm} & r_{N} < r_1 \\
401 & k_{dr}(r_{N}-r_1)\frac{\mathbf{r}_ij}{r_{ij}} & r_1 \le r_{N} < r_2 \\
402 & k_{dr}(r_2-r_1)\frac{\mathbf{r}_ij}{r_{ij}} & r_{N} \ge r_2
405 where :math:`i` and :math:`j` denote the atoms of all the pairs that
406 contribute to the NOE signal.
408 Using distance restraints
409 ^^^^^^^^^^^^^^^^^^^^^^^^^
411 A list of distance restrains based on NOE data can be added to a
412 molecule definition in your topology file, like in the following
417 [ distance_restraints ]
418 ; ai aj type index type' low up1 up2 fac
419 10 16 1 0 1 0.0 0.3 0.4 1.0
420 10 28 1 1 1 0.0 0.3 0.4 1.0
421 10 46 1 1 1 0.0 0.3 0.4 1.0
422 16 22 1 2 1 0.0 0.3 0.4 2.5
423 16 34 1 3 1 0.0 0.5 0.6 1.0
425 In this example a number of features can be found. In columns ai and aj
426 you find the atom numbers of the particles to be restrained. The type
427 column should always be 1. As explained in :ref:`distancerestraint`,
428 multiple distances can contribute to a single NOE signal. In the
429 topology this can be set using the index column. In our example, the
430 restraints 10-28 and 10-46 both have index 1, therefore they are treated
431 simultaneously. An extra requirement for treating restraints together is
432 that the restraints must be on successive lines, without any other
433 intervening restraint. The type’ column will usually be 1, but can be
434 set to 2 to obtain a distance restraint that will never be time- and
435 ensemble-averaged; this can be useful for restraining hydrogen bonds.
436 The columns ``low``, ``up1``, and
437 ``up2`` hold the values of :math:`r_0`, :math:`r_1`, and
438 :math:`r_2` from :eq:`eqn. %s <eqndisre>`. In some cases it
439 can be useful to have different force constants for some restraints;
440 this is controlled by the column ``fac``. The force constant
441 in the parameter file is multiplied by the value in the column
442 ``fac`` for each restraint. Information for each restraint
443 is stored in the energy file and can be processed and plotted with
446 Orientation restraints
447 ~~~~~~~~~~~~~~~~~~~~~~
449 This section describes how orientations between vectors, as measured in
450 certain NMR experiments, can be calculated and restrained in MD
451 simulations. The presented refinement methodology and a comparison of
452 results with and without time and ensemble averaging have been
453 published \ :ref:`92 <refHess2003>`.
458 In an NMR experiment, orientations of vectors can be measured when a
459 molecule does not tumble completely isotropically in the solvent. Two
460 examples of such orientation measurements are residual dipolar couplings
461 (between two nuclei) or chemical shift anisotropies. An observable for a
462 vector :math:`\mathbf{r}_i` can be written as follows:
464 .. math:: \delta_i = \frac{2}{3} \mbox{tr}({{\mathbf S}}{{\mathbf D}}_i)
466 where :math:`{{\mathbf S}}` is the dimensionless order tensor of the
467 molecule. The tensor :math:`{{\mathbf D}}_i` is given by:
469 .. math:: {{\mathbf D}}_i = \frac{c_i}{\|\mathbf{r}_i\|^\alpha} \left(
471 3 x x - 1 & 3 x y & 3 x z \\
472 3 x y & 3 y y - 1 & 3 y z \\
473 3 x z & 3 y z & 3 z z - 1 \\
480 x=\frac{r_{i,x}}{\|\mathbf{r}_i\|}, \quad
481 y=\frac{r_{i,y}}{\|\mathbf{r}_i\|}, \quad
482 z=\frac{r_{i,z}}{\|\mathbf{r}_i\|}
484 For a dipolar coupling :math:`\mathbf{r}_i` is the vector
485 connecting the two nuclei, :math:`\alpha=3` and the constant :math:`c_i`
488 .. math:: c_i = \frac{\mu_0}{4\pi} \gamma_1^i \gamma_2^i \frac{\hbar}{4\pi}
490 where :math:`\gamma_1^i` and :math:`\gamma_2^i` are the gyromagnetic
491 ratios of the two nuclei.
493 The order tensor is symmetric and has trace zero. Using a rotation
494 matrix :math:`{\mathbf T}` it can be transformed into the following
499 {\mathbf T}^T {{\mathbf S}}{\mathbf T} = s \left( \begin{array}{ccc}
500 -\frac{1}{2}(1-\eta) & 0 & 0 \\
501 0 & -\frac{1}{2}(1+\eta) & 0 \\
505 where :math:`-1 \leq s \leq 1` and :math:`0 \leq \eta \leq 1`.
506 :math:`s` is called the order parameter and :math:`\eta` the asymmetry
507 of the order tensor :math:`{{\mathbf S}}`. When the molecule tumbles
508 isotropically in the solvent, :math:`s` is zero, and no orientational
509 effects can be observed because all :math:`\delta_i` are zero.
511 Calculating orientations in a simulation
512 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
514 For reasons which are explained below, the :math:`{{\mathbf D}}`
515 matrices are calculated which respect to a reference orientation of the
516 molecule. The orientation is defined by a rotation matrix
517 :math:`{{\mathbf R}}`, which is needed to least-squares fit the current
518 coordinates of a selected set of atoms onto a reference conformation.
519 The reference conformation is the starting conformation of the
520 simulation. In case of ensemble averaging, which will be treated later,
521 the structure is taken from the first subsystem. The calculated
522 :math:`{{\mathbf D}}_i^c` matrix is given by:
524 .. math:: {{\mathbf D}}_i^c(t) = {{\mathbf R}}(t) {{\mathbf D}}_i(t) {{\mathbf R}}^T(t)
527 The calculated orientation for vector :math:`i` is given by:
529 .. math:: \delta^c_i(t) = \frac{2}{3} \mbox{tr}({{\mathbf S}}(t){{\mathbf D}}_i^c(t))
531 The order tensor :math:`{{\mathbf S}}(t)` is usually unknown. A
532 reasonable choice for the order tensor is the tensor which minimizes the
533 (weighted) mean square difference between the calculated and the
534 observed orientations:
536 .. 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
539 To properly combine different types of measurements, the unit of
540 :math:`w_i` should be such that all terms are dimensionless. This means
541 the unit of :math:`w_i` is the unit of :math:`\delta_i` to the power
542 :math:`-2`. **Note** that scaling all :math:`w_i` with a constant factor
543 does not influence the order tensor.
548 Since the tensors :math:`{{\mathbf D}}_i` fluctuate rapidly in time,
549 much faster than can be observed in an experiment, they should be
550 averaged over time in the simulation. However, in a simulation the time
551 and the number of copies of a molecule are limited. Usually one can not
552 obtain a converged average of the :math:`{{\mathbf D}}_i` tensors over
553 all orientations of the molecule. If one assumes that the average
554 orientations of the :math:`\mathbf{r}_i` vectors within
555 the molecule converge much faster than the tumbling time of the
556 molecule, the tensor can be averaged in an axis system that rotates with
557 the molecule, as expressed by :eq:`equation %s <eqnDrot>`). The time-averaged
558 tensors are calculated using an exponentially decaying memory function:
562 {{\mathbf D}}^a_i(t) = \frac{\displaystyle
563 \int_{u=t_0}^t {{\mathbf D}}^c_i(u) \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
565 \int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
568 Assuming that the order tensor :math:`{{\mathbf S}}` fluctuates slower
569 than the :math:`{{\mathbf D}}_i`, the time-averaged orientation can be
572 .. math:: \delta_i^a(t) = \frac{2}{3} \mbox{tr}({{\mathbf S}}(t) {{\mathbf D}}_i^a(t))
574 where the order tensor :math:`{{\mathbf S}}(t)` is calculated using
575 expression :eq:`%s <eqnSmsd>` with :math:`\delta_i^c(t)` replaced by
576 :math:`\delta_i^a(t)`.
581 The simulated structure can be restrained by applying a force
582 proportional to the difference between the calculated and the
583 experimental orientations. When no time averaging is applied, a proper
584 potential can be defined as:
586 .. math:: V = \frac{1}{2} k \sum_{i=1}^N w_i (\delta_i^c (t) -\delta_i^{exp})^2
588 where the unit of :math:`k` is the unit of energy. Thus the effective
589 force constant for restraint :math:`i` is :math:`k w_i`. The forces are
590 given by minus the gradient of :math:`V`. The force
591 :math:`\mathbf{F}\!_i` working on vector
592 :math:`\mathbf{r}_i` is:
598 & = & - \frac{\mbox{d} V}{\mbox{d}\mathbf{r}_i} \\
599 & = & -k w_i (\delta_i^c (t) -\delta_i^{exp}) \frac{\mbox{d} \delta_i (t)}{\mbox{d}\mathbf{r}_i} \\
600 & = & -k w_i (\delta_i^c (t) -\delta_i^{exp})
601 \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}
606 Ensemble averaging can be applied by simulating a system of :math:`M`
607 subsystems that each contain an identical set of orientation restraints.
608 The systems only interact via the orientation restraint potential which
613 V = M \frac{1}{2} k \sum_{i=1}^N w_i
614 \langle \delta_i^c (t) -\delta_i^{exp} \rangle^2
616 The force on vector :math:`\mathbf{r}_{i,m}` in subsystem
617 :math:`m` is given by:
621 \mathbf{F}\!_{i,m}(t) = - \frac{\mbox{d} V}{\mbox{d}\mathbf{r}_{i,m}} =
622 -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}} \\
627 When using time averaging it is not possible to define a potential. We
628 can still define a quantity that gives a rough idea of the energy stored
633 V = M \frac{1}{2} k^a \sum_{i=1}^N w_i
634 \langle \delta_i^a (t) -\delta_i^{exp} \rangle^2
636 The force constant :math:`k_a` is switched on slowly to compensate for
637 the lack of history at times close to :math:`t_0`. It is exactly
638 proportional to the amount of average that has been accumulated:
643 k \, \frac{1}{\tau}\int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
645 What really matters is the definition of the force. It is chosen to be
646 proportional to the square root of the product of the time-averaged and
647 the instantaneous deviation. Using only the time-averaged deviation
648 induces large oscillations. The force is given by:
652 \mathbf{F}\!_{i,m}(t) =
653 \left\{ \begin{array}{ll}
654 0 & \quad \mbox{for} \quad a\, b \leq 0 \\
656 k^a w_i \frac{a}{|a|} \sqrt{a\, b} \, \frac{\mbox{d} \delta_{i,m}^c (t)}{\mbox{d}\mathbf{r}_{i,m}}
657 & \quad \mbox{for} \quad a\, b > 0
664 a &=& \langle \delta_i^a (t) -\delta_i^{exp} \rangle \\
665 b &=& \langle \delta_i^c (t) -\delta_i^{exp} \rangle\end{aligned}
667 Using orientation restraints
668 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
670 Orientation restraints can be added to a molecule definition in the
671 topology file in the section ``[ orientation_restraints ]``.
672 Here we give an example section containing five N-H residual dipolar
677 [ orientation_restraints ]
678 ; ai aj type exp. label alpha const. obs. weight
680 31 32 1 1 3 3 6.083 -6.73 1.0
681 43 44 1 1 4 3 6.083 -7.87 1.0
682 55 56 1 1 5 3 6.083 -7.13 1.0
683 65 66 1 1 6 3 6.083 -2.57 1.0
684 73 74 1 1 7 3 6.083 -2.10 1.0
686 The unit of the observable is Hz, but one can choose any other unit. In
687 columns ``ai`` and ``aj`` you find the atom numbers of the particles to be
688 restrained. The ``type`` column should always be 1. The ``exp.`` column denotes
689 the experiment number, starting at 1. For each experiment a separate
690 order tensor :math:`{{\mathbf S}}` is optimized. The label should be a
691 unique number larger than zero for each restraint. The ``alpha`` column
692 contains the power :math:`\alpha` that is used in
693 :eq:`equation %s <eqnorientdef>`) to calculate the orientation. The ``const.`` column
694 contains the constant :math:`c_i` used in the same equation. The
695 constant should have the unit of the observable times
696 nm\ :math:`^\alpha`. The column ``obs.`` contains the observable, in any
697 unit you like. The last column contains the weights :math:`w_i`; the
698 unit should be the inverse of the square of the unit of the observable.
700 Some parameters for orientation restraints can be specified in the
701 :ref:`grompp <gmx grompp>` :ref:`mdp` file, for a study of the effect of different
702 force constants and averaging times and ensemble averaging see \ :ref:`92 <refHess2003>`.
703 Information for each restraint is stored in the energy
704 file and can be processed and plotted with :ref:`gmx nmr`.