Bonded interactions#
Bonded interactions are based on a fixed list of atoms. They are not exclusively pair interactions, but include 3- and 4-body interactions as well. There are bond stretching (2-body), bond angle (3-body), and dihedral angle (4-body) interactions. A special type of dihedral interaction (called improper dihedral) is used to force atoms to remain in a plane or to prevent transition to a configuration of opposite chirality (a mirror image).
Bond stretching#
Harmonic potential#
The bond stretching between two covalently bonded atoms
See also Fig. 19, with the force given by:
Fourth power potential#
In the GROMOS-96 force field 77, the covalent bond potential is, for reasons of computational efficiency, written as:
The corresponding force is:
The force constants for this form of the potential are related to the
usual harmonic force constant
The force constants are mostly derived from the harmonic ones used in
GROMOS-87 78. Although this form is
computationally more efficient (because no square root has to be
evaluated), it is conceptually more complex. One particular disadvantage
is that since the form is not harmonic, the average energy of a single
bond is not equal to
Morse potential bond stretching#
For some systems that require an anharmonic bond stretching potential, the Morse potential 79 between two atoms i and j is available in GROMACS. This potential differs from the harmonic potential in that it has an asymmetric potential well and a zero force at infinite distance. The functional form is:
See also Fig. 20, and the corresponding force is:
where
and because
For small deviations
and substituting (177) and (178) in the functional form:
we recover the harmonic bond stretching potential.
Cubic bond stretching potential#
Another anharmonic bond stretching potential that is slightly simpler than the Morse potential adds a cubic term in the distance to the simple harmonic form:
A flexible water model (based on the SPC water model 80)
including a cubic bond stretching potential for the O-H bond was
developed by Ferguson 81. This model was found to yield a
reasonable infrared spectrum. The Ferguson water model is available in
the GROMACS library (flexwat-ferguson.itp
). It should be noted that the
potential is asymmetric: overstretching leads to infinitely low
energies. The integration timestep is therefore limited to 1 fs.
The force corresponding to this potential is:
FENE bond stretching potential#
In coarse-grained polymer simulations the beads are often connected by a FENE (finitely extensible nonlinear elastic) potential 82:
The potential looks complicated, but the expression for the force is simpler:
At short distances the potential asymptotically goes to a harmonic
potential with force constant
Harmonic angle potential#
The bond-angle vibration between a triplet of atoms
As the bond-angle vibration is represented by a harmonic potential, the form is the same as the bond stretching (Fig. 19).
The force equations are given by the chain rule:
The numbering
Cosine based angle potential#
In the GROMOS-96 force field a simplified function is used to represent angle vibrations:
where
The corresponding force can be derived by partial differentiation with
respect to the atomic positions. The force constants in this function
are related to the force constants in the harmonic form
In the GROMOS-96 manual there is a much more complicated conversion formula which is temperature dependent. The formulas are equivalent at 0 K and the differences at 300 K are on the order of 0.1 to 0.2%. Note that in the input in topology files, angles are given in degrees and force constants in kJ/mol.
Restricted bending potential#
The restricted bending (ReB) potential 83 prevents the
bending angle
To systematically hinder the bending angles from reaching the
Figure 22 shows the comparison between the ReB potential, (189), and the standard one (186).
The wall of the ReB potential is very repulsive in the region close to
Due to its construction, the restricted bending potential cannot be
used for equilibrium
Urey-Bradley potential#
The Urey-Bradley bond-angle vibration between a triplet of atoms
The force equations can be deduced from sections Harmonic potential and Harmonic angle potential.
Linear Angle potential#
The linear angle potential was designed especially for linear compounds such as nitriles and for carbon dioxide 190. It avoids the calculation of the angle per se, since the angle force is not well-defined if the angle is 180 degrees. Rather, it computes the deviation of a central atom in a triplet i,j,k from a reference position
where a is defined by the bond-length i-j and j-k, in a symmetric
molecule such as carbon dioxide a = 1/2. If the compound has different
bond lengths
If the order of atoms is changed to k,j,i, a needs to be replaced by 1-a. The energy is now given by
with
Bond-Bond cross term#
The bond-bond cross term for three particles
where
The force on atom
Bond-Angle cross term#
The bond-angle cross term for three particles
where
Quartic angle potential#
For special purposes there is an angle potential that uses a fourth order polynomial:
Improper dihedrals#
Improper dihedrals are meant to keep planar groups (e.g. aromatic rings) planar, or to prevent molecules from flipping over to their mirror images, see Fig. 23.
Improper dihedrals: harmonic type#
The simplest improper dihedral potential is a harmonic potential; it is plotted in Fig. 25.
Since the potential is harmonic it is discontinuous, but since the
discontinuity is chosen at 180
Improper dihedrals: periodic type#
This potential is identical to the periodic proper dihedral (see below). There is a separate dihedral type for this (type 4) only to be able to distinguish improper from proper dihedrals in the parameter section and the output.
Proper dihedrals#
For the normal dihedral interaction there is a choice of either the
GROMOS periodic function or a function based on expansion in powers of
Proper dihedrals: periodic type#
Proper dihedral angles are defined according to the IUPAC/IUB
convention, where [ dihedral ]
section when multiple parameters are
defined for the same atomtypes in the [ dihedraltypes ]
section.
Proper dihedrals: Ryckaert-Bellemans function#
(199)#
An example of constants for
9.28 |
-13.12 |
26.24 |
|||
12.16 |
-3.06 |
-31.5 |
(Note: The use of this potential implies exclusion of LJ
interactions between the first and the last atom of the dihedral, and
(200)#
(201)#
Proper dihedrals: Fourier function#
(202)#
Proper dihedrals: Restricted torsion potential#
In a manner very similar to the restricted bending potential (see Restricted bending potential), a restricted torsion/dihedral potential is introduced:
with the advantages of being a function of
Proper dihedrals: Combined bending-torsion potential#
When the four particles forming the dihedral angle become collinear
(this situation will never happen in atomistic simulations, but it can
occur in coarse-grained simulations) the calculation of the torsion
angle and potential leads to numerical instabilities. One way to avoid
this is to use the restricted bending potential (see Restricted bending potential) that
prevents the dihedral from reaching the
Another way is to disregard any effects of the dihedral becoming ill-defined, keeping the dihedral force and potential calculation continuous in entire angle range by coupling the torsion potential (in a cosine form) with the bending potentials of the adjacent bending angles in a unique expression:
This combined bending-torsion (CBT) potential has been proposed by 88 for polymer melt simulations and is extensively described in 83.
This potential has two main advantages:
it does not only depend on the dihedral angle
(between the , , and beads) but also on the bending angles and defined from three adjacent beads ( , and , and , and , respectively). The two pre-factors, tentatively suggested by 89 and theoretically discussed by 90, cancel the torsion potential and force when either of the two bending angles approaches the value of .its dependence on
is expressed through a polynomial in that avoids the singularities in or in calculating the torsional force.
These two properties make the CBT potential well-behaved for MD
simulations with weak constraints on the bending angles or even for
steered / non-equilibrium MD in which the bending and torsion angles
suffer major modifications. When using the CBT potential, the bending
potentials for the adjacent
Additionally, the derivative of
The CBT is based on a cosine form without multiplicity, so it can only
be symmetrical around
Bonded pair and 1-4 interactions#
Most force fields do not use normal Lennard-Jones and Coulomb interactions for atoms separated by three bonds, the so-called 1-4 interactions. These interactions are still affected by the modified electronic distributions due to the chemical bonds and they are modified in the force field by the dihedral terms. For this reason the Lennard-Jones and Coulomb 1-4 interactions are often scaled down, by a fixed factor given by the force field. These factors can be supplied in the topology and the parameters can also be overriden per 1-4 interaction or atom type pair. The pair interactions can be used for any atom pair in a molecule, not only 1-4 pairs. The non-bonded interactions between such pairs should be excluded to avoid double interactions. Plain Lennard-Jones and Coulomb interactions are used which are not affected by the non-bonded interaction treatment and potential modifiers.
Tabulated bonded interaction functions#
(206)#