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
34 :label: eqnposrestform
36 The potential is plotted in :numref:`Fig. %s <fig-positionrestraint>`.
38 .. _fig-positionrestraint:
40 .. figure:: plots/f-pr.*
43 Position restraint potential.
45 The potential form can be rewritten without loss of generality as:
47 .. 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]
48 :label: eqnposrestgeneral
52 .. math:: \begin{array}{rcl}
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)
57 :label: eqnposrestforce
59 Using three different force constants the position restraints can be
60 turned on or off in each spatial dimension; this means that atoms can be
61 harmonically restrained to a plane or a line. Position restraints are
62 applied to a special fixed list of atoms. Such a list is usually
63 generated by the :ref:`pdb2gmx <gmx pdb2gmx>` program.
65 .. _reference-manual-position-restraints:
67 Note that position restraints make the potential dependent on absolute
68 coordinates in space. Therefore, in general the pressure (and virial)
69 is not well defined, as the pressure is the derivative of the free-energy
70 of the system with respect to the volume. When the reference coordinates
71 are scaled along with the system, which can be selected with the mdp option
72 :mdp-value:`refcoord-scaling=all`, the pressure and virial are well defined.
74 Flat-bottomed position restraints
75 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
77 Flat-bottomed position restraints can be used to restrain particles to
78 part of the simulation volume. No force acts on the restrained particle
79 within the flat-bottomed region of the potential, however a harmonic
80 force acts to move the particle to the flat-bottomed region if it is
81 outside it. It is possible to apply normal and flat-bottomed position
82 restraints on the same particle (however, only with the same reference
83 position :math:`\mathbf{R}_i`). The following general
84 potential is used (:numref:`Figure %s <fig-fbposres>` A):
86 .. 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}],
87 :label: eqnflatbottomposrest
89 where :math:`\mathbf{R}_i` is the reference position,
90 :math:`r_\mathrm{fb}` is the distance from the center with a flat
91 potential, :math:`k_\mathrm{fb}` the force constant, and :math:`H` is
92 the Heaviside step function. The distance
93 :math:`d_g(\mathbf{r}_i;\mathbf{R}_i)` from
94 the reference position depends on the geometry :math:`g` of the
95 flat-bottomed potential.
99 .. figure:: plots/fbposres.*
102 Flat-bottomed position restraint potential. (A) Not inverted, (B)
105 | The following geometries for the flat-bottomed potential are
108 | **Sphere** (:math:`g =1`): The
109 particle is kept in a sphere of given radius. The force acts towards
110 the center of the sphere. The following distance calculation is used:
112 .. math:: d_g(\mathbf{r}_i;\mathbf{R}_i) = | \mathbf{r}_i-\mathbf{R}_i |
113 :label: eqnfbsphereposrest
115 | **Cylinder** (:math:`g=6,7,8`): The particle is kept in a cylinder of
116 given radius parallel to the :math:`x` (:math:`g=6`), :math:`y`
117 (:math:`g=7`), or :math:`z`-axis (:math:`g=8`). For backwards
118 compatibility, setting :math:`g=2` is mapped to :math:`g=8` in the
119 code so that old :ref:`tpr` files and topologies work. The
120 force from the flat-bottomed potential acts towards the axis of the
121 cylinder. The component of the force parallel to the cylinder axis is
122 zero. For a cylinder aligned along the :math:`z`-axis:
124 .. math:: d_g(\mathbf{r}_i;\mathbf{R}_i) = \sqrt{ (x_i-X_i)^2 + (y_i - Y_i)^2 }
125 :label: eqnfbcylinderposrest
127 | **Layer** (:math:`g=3,4,5`): The particle is kept in a layer defined
128 by the thickness and the normal of the layer. The layer normal can be
129 parallel to the :math:`x`, :math:`y`, or :math:`z`-axis. The force
130 acts parallel to the layer normal.
132 .. math:: d_g(\mathbf{r}_i;\mathbf{R}_i) = |x_i-X_i|, \;\;\;\mbox{or}\;\;\;
133 d_g(\mathbf{r}_i;\mathbf{R}_i) = |y_i-Y_i|, \;\;\;\mbox{or}\;\;\;
134 d_g(\mathbf{r}_i;\mathbf{R}_i) = |z_i-Z_i|.
135 :label: eqnfblayerposrest
137 It is possible to apply multiple independent flat-bottomed position
138 restraints of different geometry on one particle. For example, applying
139 a cylinder and a layer in :math:`z` keeps a particle within a disk.
140 Applying three layers in :math:`x`, :math:`y`, and :math:`z` keeps the
141 particle within a cuboid.
143 In addition, it is possible to invert the restrained region with the
144 unrestrained region, leading to a potential that acts to keep the
145 particle *outside* of the volume defined by
146 :math:`\mathbf{R}_i`, :math:`g`, and
147 :math:`r_\mathrm{fb}`. That feature is switched on by defining a
148 negative :math:`r_\mathrm{fb}` in the topology. The following potential
149 is used (:numref:`Figure %s <fig-fbposres>` B):
151 .. math:: V_\mathrm{fb}^{\mathrm{inv}}(\mathbf{r}_i) = \frac{1}{2}k_\mathrm{fb}
152 [d_g(\mathbf{r}_i;\mathbf{R}_i) - | r_\mathrm{fb} | ]^2\,
153 H[ -(d_g(\mathbf{r}_i;\mathbf{R}_i) - | r_\mathrm{fb} | )].
154 :label: eqninvertrest
159 These are used to restrain the angle between two pairs of particles or
160 between one pair of particles and the :math:`z`-axis. The functional
161 form is similar to that of a proper dihedral. For two pairs of atoms:
163 .. math:: V_{ar}(\mathbf{r}_i,\mathbf{r}_j,\mathbf{r}_k,\mathbf{r}_l)
164 = k_{ar}(1 - \cos(n (\theta - \theta_0))
167 \theta = \arccos\left(\frac{\mathbf{r}_j -\mathbf{r}_i}{\|\mathbf{r}_j -\mathbf{r}_i\|}
168 \cdot \frac{\mathbf{r}_l -\mathbf{r}_k}{\|\mathbf{r}_l -\mathbf{r}_k\|} \right)
171 For one pair of atoms and the :math:`z`-axis:
173 .. math:: V_{ar}(\mathbf{r}_i,\mathbf{r}_j) = k_{ar}(1 - \cos(n (\theta - \theta_0))
176 \theta = \arccos\left(\frac{\mathbf{r}_j -\mathbf{r}_i}{\|\mathbf{r}_j -\mathbf{r}_i\|}
177 \cdot \left( \begin{array}{c} 0 \\ 0 \\ 1 \\ \end{array} \right) \right)
178 :label: eqnanglerestzaxis
180 A multiplicity (:math:`n`) of 2 is useful when you do not want to
181 distinguish between parallel and anti-parallel vectors. The equilibrium
182 angle :math:`\theta` should be between 0 and 180 degrees for
183 multiplicity 1 and between 0 and 90 degrees for multiplicity 2.
185 .. _dihedralrestraint:
190 These are used to restrain the dihedral angle :math:`\phi` defined by
191 four particles as in an improper dihedral (sec. :ref:`imp`) but with a
192 slightly modified potential. Using:
194 .. math:: \phi' = \left(\phi-\phi_0\right) ~{\rm MOD}~ 2\pi
197 where :math:`\phi_0` is the reference angle, the potential is defined
200 .. math:: V_{dihr}(\phi') ~=~ \left\{
201 \begin{array}{lcllll}
202 {\frac{1}{2}}k_{dihr}(\phi'-\Delta\phi)^2
203 &\mbox{for}& \|\phi'\| & > & \Delta\phi \\[1.5ex]
204 0 &\mbox{for}& \|\phi'\| & \le & \Delta\phi \\[1.5ex]
208 where :math:`\Delta\phi` is a user defined angle and :math:`k_{dihr}`
209 is the force constant. **Note** that in the input in topology files,
210 angles are given in degrees and force constants in
211 kJ/mol/rad\ :math:`^2`.
213 .. _distancerestraint:
218 Distance restraints add a penalty to the potential when the distance
219 between specified pairs of atoms exceeds a threshold value. They are
220 normally used to impose experimental restraints from, for instance,
221 experiments in nuclear magnetic resonance (NMR), on the motion of the
222 system. Thus, MD can be used for structure refinement using NMR data. In
223 |Gromacs| there are three ways to impose restraints on pairs of atoms:
225 - Simple harmonic restraints: use ``[ bonds ]`` type 6 (see sec. :ref:`excl`).
227 - Piecewise linear/harmonic restraints: ``[ bonds ]`` type
230 - Complex NMR distance restraints, optionally with pair, time and/or
233 The last two options will be detailed now.
235 The potential form for distance restraints is quadratic below a
236 specified lower bound and between two specified upper bounds, and linear
237 beyond the largest bound (see :numref:`Fig. %s <fig-dist>`).
239 .. math:: V_{dr}(r_{ij}) ~=~ \left\{
240 \begin{array}{lcllllll}
241 {\frac{1}{2}}k_{dr}(r_{ij}-r_0)^2
242 &\mbox{for}& & & r_{ij} & < & r_0 \\[1.5ex]
243 0 &\mbox{for}& r_0 & \le & r_{ij} & < & r_1 \\[1.5ex]
244 {\frac{1}{2}}k_{dr}(r_{ij}-r_1)^2
245 &\mbox{for}& r_1 & \le & r_{ij} & < & r_2 \\[1.5ex]
246 {\frac{1}{2}}k_{dr}(r_2-r_1)(2r_{ij}-r_2-r_1)
247 &\mbox{for}& r_2 & \le & r_{ij} & &
253 .. figure:: plots/f-dr.*
256 Distance Restraint potential.
260 .. math:: \mathbf{F}_i~=~ \left\{
261 \begin{array}{lcllllll}
262 -k_{dr}(r_{ij}-r_0)\frac{\mathbf{r}_{ij}}{r_{ij}}
263 &\mbox{for}& & & r_{ij} & < & r_0 \\[1.5ex]
264 0 &\mbox{for}& r_0 & \le & r_{ij} & < & r_1 \\[1.5ex]
265 -k_{dr}(r_{ij}-r_1)\frac{\mathbf{r}_{ij}}{r_{ij}}
266 &\mbox{for}& r_1 & \le & r_{ij} & < & r_2 \\[1.5ex]
267 -k_{dr}(r_2-r_1)\frac{\mathbf{r}_{ij}}{r_{ij}}
268 &\mbox{for}& r_2 & \le & r_{ij} & &
270 :label: eqndisreforce
272 For restraints not derived from NMR data, this functionality will
273 usually suffice and a section of ``[ bonds ]`` type 10 can be used to apply individual
274 restraints between pairs of atoms, see :ref:`topfile`. For applying
275 restraints derived from NMR measurements, more complex functionality
276 might be required, which is provided through the ``[ distance_restraints ]`` section and is
282 Distance restraints based on instantaneous distances can potentially
283 reduce the fluctuations in a molecule significantly. This problem can be
284 overcome by restraining to a *time averaged*
285 distance \ :ref:`91 <refTorda89>`. The forces with time averaging are:
287 .. math:: \mathbf{F}_i~=~ \left\{
288 \begin{array}{lcllllll}
289 -k^a_{dr}(\bar{r}_{ij}-r_0)\frac{\mathbf{r}_{ij}}{r_{ij}}
290 &\mbox{for}& & & \bar{r}_{ij} & < & r_0 \\[1.5ex]
291 0 &\mbox{for}& r_0 & \le & \bar{r}_{ij} & < & r_1 \\[1.5ex]
292 -k^a_{dr}(\bar{r}_{ij}-r_1)\frac{\mathbf{r}_{ij}}{r_{ij}}
293 &\mbox{for}& r_1 & \le & \bar{r}_{ij} & < & r_2 \\[1.5ex]
294 -k^a_{dr}(r_2-r_1)\frac{\mathbf{r}_{ij}}{r_{ij}}
295 &\mbox{for}& r_2 & \le & \bar{r}_{ij} & &
297 :label: eqntimeaveragerest
299 where :math:`\bar{r}_{ij}` is given by an exponential running average
300 with decay time :math:`\tau`:
302 .. math:: \bar{r}_{ij} ~=~ < r_{ij}^{-3} >^{-1/3}
305 The force constant :math:`k^a_{dr}` is switched on slowly to compensate
306 for the lack of history at the beginning of the simulation:
308 .. math:: k^a_{dr} = k_{dr} \left(1-\exp\left(-\frac{t}{\tau}\right)\right)
309 :label: eqnforceconstantswitch
311 Because of the time averaging, we can no longer speak of a distance
314 This way an atom can satisfy two incompatible distance restraints *on
315 average* by moving between two positions. An example would be an amino
316 acid side-chain that is rotating around its :math:`\chi` dihedral angle,
317 thereby coming close to various other groups. Such a mobile side chain
318 can give rise to multiple NOEs that can not be fulfilled by a single
321 The computation of the time averaged distance in the
322 :ref:`mdrun <gmx mdrun>` program is done in the following fashion:
324 .. math:: \begin{array}{rcl}
325 \overline{r^{-3}}_{ij}(0) &=& r_{ij}(0)^{-3} \\
326 \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]
330 When a pair is within the bounds, it can still feel a force because the
331 time averaged distance can still be beyond a bound. To prevent the
332 protons from being pulled too close together, a mixed approach can be
333 used. In this approach, the penalty is zero when the instantaneous
334 distance is within the bounds, otherwise the violation is the square
335 root of the product of the instantaneous violation and the time averaged
338 .. math:: \mathbf{F}_i~=~ \left\{
340 k^a_{dr}\sqrt{(r_{ij}-r_0)(\bar{r}_{ij}-r_0)}\frac{\mathbf{r}_{ij}}{r_{ij}}
341 & \mbox{for} & r_{ij} < r_0 & \mbox{and} & \bar{r}_{ij} < r_0 \\[1.5ex]
343 \mbox{min}\left(\sqrt{(r_{ij}-r_1)(\bar{r}_{ij}-r_1)},r_2-r_1\right)
344 \frac{\mathbf{r}_{ij}}{r_{ij}}
345 & \mbox{for} & r_{ij} > r_1 & \mbox{and} & \bar{r}_{ij} > r_1 \\[1.5ex]
348 :label: eqntimeaverageviolation
350 Averaging over multiple pairs
351 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
353 Sometimes it is unclear from experimental data which atom pair gives
354 rise to a single NOE, in other occasions it can be obvious that more
355 than one pair contributes due to the symmetry of the system, *e.g.* a
356 methyl group with three protons. For such a group, it is not possible to
357 distinguish between the protons, therefore they should all be taken into
358 account when calculating the distance between this methyl group and
359 another proton (or group of protons). Due to the physical nature of
360 magnetic resonance, the intensity of the NOE signal is inversely
361 proportional to the sixth power of the inter-atomic distance. Thus, when
362 combining atom pairs, a fixed list of :math:`N` restraints may be taken
363 together, where the apparent “distance” is given by:
365 .. math:: r_N(t) = \left [\sum_{n=1}^{N} \bar{r}_{n}(t)^{-6} \right]^{-1/6}
368 where we use :math:`r_{ij}` or :eq:`eqn. %s <eqnrav>` for the
369 :math:`\bar{r}_{n}`. The :math:`r_N` of the instantaneous and
370 time-averaged distances can be combined to do a mixed restraining, as
371 indicated above. As more pairs of protons contribute to the same NOE
372 signal, the intensity will increase, and the summed “distance” will be
373 shorter than any of its components due to the reciprocal summation.
375 There are two options for distributing the forces over the atom pairs.
376 In the conservative option, the force is defined as the derivative of
377 the restraint potential with respect to the coordinates. This results in
378 a conservative potential when time averaging is not used. The force
379 distribution over the pairs is proportional to :math:`r^{-6}`. This
380 means that a close pair feels a much larger force than a distant pair,
381 which might lead to a molecule that is “too rigid.” The other option is
382 an equal force distribution. In this case each pair feels :math:`1/N` of
383 the derivative of the restraint potential with respect to :math:`r_N`.
384 The advantage of this method is that more conformations might be
385 sampled, but the non-conservative nature of the forces can lead to local
386 heating of the protons.
388 It is also possible to use *ensemble averaging* using multiple (protein)
389 molecules. In this case the bounds should be lowered as in:
391 .. math:: \begin{array}{rcl}
392 r_1 &~=~& r_1 * M^{-1/6} \\
393 r_2 &~=~& r_2 * M^{-1/6}
395 :label: eqnrestforceensembleaverage
397 where :math:`M` is the number of molecules. The |Gromacs| preprocessor
398 :ref:`grompp <gmx grompp>` can do this automatically when the appropriate
399 option is given. The resulting “distance” is then used to calculate the
400 scalar force according to:
402 .. math:: \mathbf{F}_i~=~\left\{
404 ~& 0 \hspace{4cm} & r_{N} < r_1 \\
405 & k_{dr}(r_{N}-r_1)\frac{\mathbf{r}_{ij}}{r_{ij}} & r_1 \le r_{N} < r_2 \\
406 & k_{dr}(r_2-r_1)\frac{\mathbf{r}_{ij}}{r_{ij}} & r_{N} \ge r_2
408 :label: eqnrestscalarforce
410 where :math:`i` and :math:`j` denote the atoms of all the pairs that
411 contribute to the NOE signal.
413 Using distance restraints
414 ^^^^^^^^^^^^^^^^^^^^^^^^^
416 A list of distance restrains based on NOE data can be added to a
417 molecule definition in your topology file, like in the following
422 [ distance_restraints ]
423 ; ai aj type index type' low up1 up2 fac
424 10 16 1 0 1 0.0 0.3 0.4 1.0
425 10 28 1 1 1 0.0 0.3 0.4 1.0
426 10 46 1 1 1 0.0 0.3 0.4 1.0
427 16 22 1 2 1 0.0 0.3 0.4 2.5
428 16 34 1 3 1 0.0 0.5 0.6 1.0
430 In this example a number of features can be found. In columns ai and aj
431 you find the atom numbers of the particles to be restrained. The type
432 column should always be 1. As explained in :ref:`distancerestraint`,
433 multiple distances can contribute to a single NOE signal. In the
434 topology this can be set using the index column. In our example, the
435 restraints 10-28 and 10-46 both have index 1, therefore they are treated
436 simultaneously. An extra requirement for treating restraints together is
437 that the restraints must be on successive lines, without any other
438 intervening restraint. The type’ column will usually be 1, but can be
439 set to 2 to obtain a distance restraint that will never be time- and
440 ensemble-averaged; this can be useful for restraining hydrogen bonds.
441 The columns ``low``, ``up1``, and
442 ``up2`` hold the values of :math:`r_0`, :math:`r_1`, and
443 :math:`r_2` from :eq:`eqn. %s <eqndisre>`. In some cases it
444 can be useful to have different force constants for some restraints;
445 this is controlled by the column ``fac``. The force constant
446 in the parameter file is multiplied by the value in the column
447 ``fac`` for each restraint. Information for each restraint
448 is stored in the energy file and can be processed and plotted with
451 Orientation restraints
452 ~~~~~~~~~~~~~~~~~~~~~~
454 This section describes how orientations between vectors, as measured in
455 certain NMR experiments, can be calculated and restrained in MD
456 simulations. The presented refinement methodology and a comparison of
457 results with and without time and ensemble averaging have been
458 published \ :ref:`92 <refHess2003>`.
463 In an NMR experiment, orientations of vectors can be measured when a
464 molecule does not tumble completely isotropically in the solvent. Two
465 examples of such orientation measurements are residual dipolar couplings
466 (between two nuclei) or chemical shift anisotropies. An observable for a
467 vector :math:`\mathbf{r}_i` can be written as follows:
469 .. math:: \delta_i = \frac{2}{3} \mbox{tr}({{\mathbf S}}{{\mathbf D}}_i)
470 :label: eqnorrestvector
472 where :math:`{{\mathbf S}}` is the dimensionless order tensor of the
473 molecule. The tensor :math:`{{\mathbf D}}_i` is given by:
475 .. math:: {{\mathbf D}}_i = \frac{c_i}{\|\mathbf{r}_i\|^\alpha} \left(
477 3 x x - 1 & 3 x y & 3 x z \\
478 3 x y & 3 y y - 1 & 3 y z \\
479 3 x z & 3 y z & 3 z z - 1 \\
483 .. math:: \mbox{with:} \quad
484 x=\frac{r_{i,x}}{\|\mathbf{r}_i\|}, \quad
485 y=\frac{r_{i,y}}{\|\mathbf{r}_i\|}, \quad
486 z=\frac{r_{i,z}}{\|\mathbf{r}_i\|}
487 :label: eqnorientdef2
489 For a dipolar coupling :math:`\mathbf{r}_i` is the vector
490 connecting the two nuclei, :math:`\alpha=3` and the constant :math:`c_i`
493 .. math:: c_i = \frac{\mu_0}{4\pi} \gamma_1^i \gamma_2^i \frac{\hbar}{4\pi}
494 :label: eqnorrestconstant
496 where :math:`\gamma_1^i` and :math:`\gamma_2^i` are the gyromagnetic
497 ratios of the two nuclei.
499 The order tensor is symmetric and has trace zero. Using a rotation
500 matrix :math:`{\mathbf T}` it can be transformed into the following
503 .. math:: {\mathbf T}^T {{\mathbf S}}{\mathbf T} = s \left( \begin{array}{ccc}
504 -\frac{1}{2}(1-\eta) & 0 & 0 \\
505 0 & -\frac{1}{2}(1+\eta) & 0 \\
508 :label: eqnorresttensor
510 where :math:`-1 \leq s \leq 1` and :math:`0 \leq \eta \leq 1`.
511 :math:`s` is called the order parameter and :math:`\eta` the asymmetry
512 of the order tensor :math:`{{\mathbf S}}`. When the molecule tumbles
513 isotropically in the solvent, :math:`s` is zero, and no orientational
514 effects can be observed because all :math:`\delta_i` are zero.
516 Calculating orientations in a simulation
517 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
519 For reasons which are explained below, the :math:`{{\mathbf D}}`
520 matrices are calculated which respect to a reference orientation of the
521 molecule. The orientation is defined by a rotation matrix
522 :math:`{{\mathbf R}}`, which is needed to least-squares fit the current
523 coordinates of a selected set of atoms onto a reference conformation.
524 The reference conformation is the starting conformation of the
525 simulation. In case of ensemble averaging, which will be treated later,
526 the structure is taken from the first subsystem. The calculated
527 :math:`{{\mathbf D}}_i^c` matrix is given by:
529 .. math:: {{\mathbf D}}_i^c(t) = {{\mathbf R}}(t) {{\mathbf D}}_i(t) {{\mathbf R}}^T(t)
532 The calculated orientation for vector :math:`i` is given by:
534 .. math:: \delta^c_i(t) = \frac{2}{3} \mbox{tr}({{\mathbf S}}(t){{\mathbf D}}_i^c(t))
535 :label: eqnDrotvector
537 The order tensor :math:`{{\mathbf S}}(t)` is usually unknown. A
538 reasonable choice for the order tensor is the tensor which minimizes the
539 (weighted) mean square difference between the calculated and the
540 observed orientations:
542 .. 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
545 To properly combine different types of measurements, the unit of
546 :math:`w_i` should be such that all terms are dimensionless. This means
547 the unit of :math:`w_i` is the unit of :math:`\delta_i` to the power
548 :math:`-2`. **Note** that scaling all :math:`w_i` with a constant factor
549 does not influence the order tensor.
554 Since the tensors :math:`{{\mathbf D}}_i` fluctuate rapidly in time,
555 much faster than can be observed in an experiment, they should be
556 averaged over time in the simulation. However, in a simulation the time
557 and the number of copies of a molecule are limited. Usually one can not
558 obtain a converged average of the :math:`{{\mathbf D}}_i` tensors over
559 all orientations of the molecule. If one assumes that the average
560 orientations of the :math:`\mathbf{r}_i` vectors within
561 the molecule converge much faster than the tumbling time of the
562 molecule, the tensor can be averaged in an axis system that rotates with
563 the molecule, as expressed by :eq:`equation %s <eqnDrot>`). The time-averaged
564 tensors are calculated using an exponentially decaying memory function:
566 .. math:: {{\mathbf D}}^a_i(t) = \frac{\displaystyle
567 \int_{u=t_0}^t {{\mathbf D}}^c_i(u) \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
569 \int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
571 :label: eqnorresttimeaverage
573 Assuming that the order tensor :math:`{{\mathbf S}}` fluctuates slower
574 than the :math:`{{\mathbf D}}_i`, the time-averaged orientation can be
577 .. math:: \delta_i^a(t) = \frac{2}{3} \mbox{tr}({{\mathbf S}}(t) {{\mathbf D}}_i^a(t))
578 :label: eqnorresttimeaveorient
580 where the order tensor :math:`{{\mathbf S}}(t)` is calculated using
581 expression :eq:`%s <eqnSmsd>` with :math:`\delta_i^c(t)` replaced by
582 :math:`\delta_i^a(t)`.
587 The simulated structure can be restrained by applying a force
588 proportional to the difference between the calculated and the
589 experimental orientations. When no time averaging is applied, a proper
590 potential can be defined as:
592 .. math:: V = \frac{1}{2} k \sum_{i=1}^N w_i (\delta_i^c (t) -\delta_i^{exp})^2
593 :label: eqnorrestsimrest
595 where the unit of :math:`k` is the unit of energy. Thus the effective
596 force constant for restraint :math:`i` is :math:`k w_i`. The forces are
597 given by minus the gradient of :math:`V`. The force
598 :math:`\mathbf{F}\!_i` working on vector
599 :math:`\mathbf{r}_i` is:
601 .. math:: \begin{aligned}
603 & = & - \frac{\mbox{d} V}{\mbox{d}\mathbf{r}_i} \\
604 & = & -k w_i (\delta_i^c (t) -\delta_i^{exp}) \frac{\mbox{d} \delta_i (t)}{\mbox{d}\mathbf{r}_i} \\
605 & = & -k w_i (\delta_i^c (t) -\delta_i^{exp})
606 \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}
607 :label: eqnorrestsimrestforce
612 Ensemble averaging can be applied by simulating a system of :math:`M`
613 subsystems that each contain an identical set of orientation restraints.
614 The systems only interact via the orientation restraint potential which
617 .. math:: V = M \frac{1}{2} k \sum_{i=1}^N w_i
618 \langle \delta_i^c (t) -\delta_i^{exp} \rangle^2
619 :label: eqnorrestensembleave
621 The force on vector :math:`\mathbf{r}_{i,m}` in subsystem
622 :math:`m` is given by:
624 .. math:: \mathbf{F}\!_{i,m}(t) = - \frac{\mbox{d} V}{\mbox{d}\mathbf{r}_{i,m}} =
625 -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}}
626 :label: eqnorrestensaveforce
631 When using time averaging it is not possible to define a potential. We
632 can still define a quantity that gives a rough idea of the energy stored
635 .. math:: V = M \frac{1}{2} k^a \sum_{i=1}^N w_i
636 \langle \delta_i^a (t) -\delta_i^{exp} \rangle^2
637 :label: eqntimeavepot
639 The force constant :math:`k_a` is switched on slowly to compensate for
640 the lack of history at times close to :math:`t_0`. It is exactly
641 proportional to the amount of average that has been accumulated:
644 k \, \frac{1}{\tau}\int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
645 :label: eqntimeaveforceswitch
647 What really matters is the definition of the force. It is chosen to be
648 proportional to the square root of the product of the time-averaged and
649 the instantaneous deviation. Using only the time-averaged deviation
650 induces large oscillations. The force is given by:
652 .. math:: \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
660 :label: eqntimeaveforce
662 .. math:: \begin{aligned}
663 a &=& \langle \delta_i^a (t) -\delta_i^{exp} \rangle \\
664 b &=& \langle \delta_i^c (t) -\delta_i^{exp} \rangle\end{aligned}
665 :label: eqntimeaveforce2
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`.