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.
Position restraints¶
These are used to restrain particles to fixed reference positions
The following form is used:
The potential is plotted in Fig. 29.
The potential form can be rewritten without loss of generality as:
Now the forces are:
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 pdb2gmx program.
Note that position restraints make the potential dependent on absolute
coordinates in space. Therefore, in general the pressure (and virial)
is not well defined, as the pressure is the derivative of the free-energy
of the system with respect to the volume. When the reference coordinates
are scaled along with the system, which can be selected with the mdp option
refcoord-scaling=all
, the pressure and virial are well defined.
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
where
(210)¶
(211)¶
(212)¶
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
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
Angle restraints¶
These are used to restrain the angle between two pairs of particles or
between one pair of particles and the
For one pair of atoms and the
A multiplicity (
Dihedral restraints¶
These are used to restrain the dihedral angle
where
where
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. Exclusions).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 Fig. 31).
The forces are
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 Topology file. 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 91. The forces with time averaging are:
where
The force constant
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
The computation of the time averaged distance in the mdrun program is done in the following fashion:
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:
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
where we use
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
It is also possible to use ensemble averaging using multiple (protein) molecules. In this case the bounds should be lowered as in:
where
where
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 Distance restraints,
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 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
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 92.
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
where
For a dipolar coupling
where
The order tensor is symmetric and has trace zero. Using a rotation
matrix
where
Calculating orientations in a simulation¶
For reasons which are explained below, the
The calculated orientation for vector
The order tensor
To properly combine different types of measurements, the unit of
Time averaging¶
Since the tensors
Assuming that the order tensor
where the order tensor
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:
where the unit of
Ensemble averaging¶
Ensemble averaging can be applied by simulating a system of
The force on vector
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:
The force constant
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:
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 alpha
column
contains the power const.
column
contains the constant obs.
contains the observable, in any
unit you like. The last column contains the weights
Some parameters for orientation restraints can be specified in the grompp mdp file, for a study of the effect of different force constants and averaging times and ensemble averaging see 92. Information for each restraint is stored in the energy file and can be processed and plotted with gmx nmr.