4 Bonded interactions are based on a fixed list of atoms. They are not
5 exclusively pair interactions, but include 3- and 4-body interactions as
6 well. There are *bond stretching* (2-body), *bond angle* (3-body), and
7 *dihedral angle* (4-body) interactions. A special type of dihedral
8 interaction (called *improper dihedral*) is used to force atoms to
9 remain in a plane or to prevent transition to a configuration of
10 opposite chirality (a mirror image).
22 The bond stretching between two covalently bonded atoms :math:`i` and
23 :math:`j` is represented by a harmonic potential:
27 .. figure:: plots/bstretch.*
30 Principle of bond stretching (left), and the bond stretching
33 .. math:: V_b~({r_{ij}}) = {\frac{1}{2}}k^b_{ij}({r_{ij}}-b_{ij})^2
34 :label: eqnharmbondstretch
36 See also :numref:`Fig. %s <fig-bstretch1>`, with the force given by:
38 .. math:: \mathbf{F}_i(\mathbf{r}_{ij}) = k^b_{ij}({r_{ij}}-b_{ij}) {\frac{{\mathbf{r}_{ij}}}{{r_{ij}}}}
39 :label: eqnharmbondstretchforce
43 Fourth power potential
44 ^^^^^^^^^^^^^^^^^^^^^^
46 In the GROMOS-96 force field \ :ref:`77 <refgromos96>`, the covalent bond
47 potential is, for reasons of computational efficiency, written as:
49 .. math:: V_b~({r_{ij}}) = \frac{1}{4}k^b_{ij}\left({r_{ij}}^2-b_{ij}^2\right)^2
52 The corresponding force is:
54 .. math:: \mathbf{F}_i(\mathbf{r}_{ij}) = k^b_{ij}({r_{ij}}^2-b_{ij}^2)~\mathbf{r}_ij
55 :label: eqng96bondforce
57 The force constants for this form of the potential are related to the
58 usual harmonic force constant :math:`k^{b,\mathrm{harm}}`
59 (sec. :ref:`bondpot`) as
61 .. math:: 2 k^b b_{ij}^2 = k^{b,\mathrm{harm}}
62 :label: eqn96harmrelation
64 The force constants are mostly derived from the harmonic ones used in
65 GROMOS-87 :ref:`78 <refbiomos>`. Although this form is
66 computationally more efficient (because no square root has to be
67 evaluated), it is conceptually more complex. One particular disadvantage
68 is that since the form is not harmonic, the average energy of a single
69 bond is not equal to :math:`{\frac{1}{2}}kT` as it is for the normal
74 Morse potential bond stretching
75 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
77 For some systems that require an anharmonic bond stretching potential,
78 the Morse potential \ :ref:`79 <refMorse29>` between two atoms *i* and *j* is
79 available in |Gromacs|. This potential differs from the harmonic potential
80 in that it has an asymmetric potential well and a zero force at infinite
81 distance. The functional form is:
83 .. math:: \displaystyle V_{morse} (r_{ij}) = D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2,
86 See also :numref:`Fig. %s <fig-morse>`, and the corresponding force is:
88 .. math:: \begin{array}{rcl}
89 \displaystyle {\bf F}_{morse} ({\bf r}_{ij})&=&2 D_{ij} \beta_{ij} \exp(-\beta_{ij}(r_{ij}-b_{ij})) * \\
90 \displaystyle \: & \: &[1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))] \frac{\displaystyle {\bf r}_{ij}}{\displaystyle r_{ij}},
92 :label: eqnmorsebondforce
94 where :math:`\displaystyle D_{ij}` is the depth of the well in
95 kJ/mol, :math:`\displaystyle \beta_{ij}` defines the steepness of the
96 well (in nm\ :math:`^{-1}`), and :math:`\displaystyle b_{ij}` is the
97 equilibrium distance in nm. The steepness parameter
98 :math:`\displaystyle \beta_{ij}` can be expressed in terms of the reduced mass of the atoms *i* and
99 *j*, the fundamental vibration frequency :math:`\displaystyle\omega_{ij}` and the well depth :math:`\displaystyle D_{ij}`:
101 .. math:: \displaystyle \beta_{ij}= \omega_{ij} \sqrt{\frac{\mu_{ij}}{2 D_{ij}}}
104 and because :math:`\displaystyle \omega = \sqrt{k/\mu}`, one can
105 rewrite :math:`\displaystyle \beta_{ij}` in terms of the harmonic
106 force constant :math:`\displaystyle k_{ij}`:
108 .. math:: \displaystyle \beta_{ij}= \sqrt{\frac{k_{ij}}{2 D_{ij}}}
111 For small deviations :math:`\displaystyle (r_{ij}-b_{ij})`, one can
112 approximate the :math:`\displaystyle \exp`-term to first-order using a
115 .. math:: \displaystyle \exp(-x) \approx 1-x
118 and substituting :eq:`eqn. %s <eqnbetaij>` and :eq:`eqn. %s <eqnexpminx>` in the
121 .. math:: \begin{array}{rcl}
122 \displaystyle V_{morse} (r_{ij})&=&D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2\\
123 \displaystyle \:&=&D_{ij} [1 - (1 -\sqrt{\frac{k_{ij}}{2 D_{ij}}}(r_{ij}-b_{ij}))]^2\\
124 \displaystyle \:&=&\frac{1}{2} k_{ij} (r_{ij}-b_{ij}))^2
126 :label: eqnharmfrommorse
128 we recover the harmonic bond stretching potential.
132 .. figure:: plots/f-morse.*
135 The Morse potential well, with bond length 0.15 nm.
137 Cubic bond stretching potential
138 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
140 Another anharmonic bond stretching potential that is slightly simpler
141 than the Morse potential adds a cubic term in the distance to the simple
144 .. math:: V_b~({r_{ij}}) = k^b_{ij}({r_{ij}}-b_{ij})^2 + k^b_{ij}k^{cub}_{ij}({r_{ij}}-b_{ij})^3
147 A flexible water model (based on the SPC water model \ :ref:`80 <refBerendsen81>`)
148 including a cubic bond stretching potential for the O-H bond was
149 developed by Ferguson \ :ref:`81 <refFerguson95>`. This model was found to yield a
150 reasonable infrared spectrum. The Ferguson water model is available in
151 the |Gromacs| library (``flexwat-ferguson.itp``). It should be noted that the
152 potential is asymmetric: overstretching leads to infinitely low
153 energies. The integration timestep is therefore limited to 1 fs.
155 The force corresponding to this potential is:
157 .. math:: \mathbf{F}_i(\mathbf{r}_{ij}) = 2k^b_{ij}({r_{ij}}-b_{ij})~{\frac{{\mathbf{r}_{ij}}}{{r_{ij}}}}+ 3k^b_{ij}k^{cub}_{ij}({r_{ij}}-b_{ij})^2~{\frac{{\mathbf{r}_{ij}}}{{r_{ij}}}}
158 :label: eqncubicbondforce
160 FENE bond stretching potential
161 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
163 In coarse-grained polymer simulations the beads are often connected by a
164 FENE (finitely extensible nonlinear elastic) potential \ :ref:`82 <refWarner72>`:
166 .. math:: V_{\mbox{FENE}}({r_{ij}}) =
167 -{\frac{1}{2}}k^b_{ij} b^2_{ij} \log\left(1 - \frac{{r_{ij}}^2}{b^2_{ij}}\right)
170 The potential looks complicated, but the expression for the force is
173 .. math:: F_{\mbox{FENE}}(\mathbf{r}_{ij}) =
174 -k^b_{ij} \left(1 - \frac{{r_{ij}}^2}{b^2_{ij}}\right)^{-1} \mathbf{r}_{ij}
175 :label: eqnfenebondforce
177 At short distances the potential asymptotically goes to a harmonic
178 potential with force constant :math:`k^b`, while it diverges at distance
183 Harmonic angle potential
184 ~~~~~~~~~~~~~~~~~~~~~~~~
186 The bond-angle vibration between a triplet of atoms :math:`i` -
187 :math:`j` - :math:`k` is also represented by a harmonic potential on the
188 angle :math:`{\theta_{ijk}}`
192 .. figure:: plots/angle.*
195 Principle of angle vibration (left) and the bond angle potential.
197 .. math:: V_a({\theta_{ijk}}) = {\frac{1}{2}}k^{\theta}_{ijk}({\theta_{ijk}}-{\theta_{ijk}}^0)^2
200 As the bond-angle vibration is represented by a harmonic potential, the
201 form is the same as the bond stretching
202 (:numref:`Fig. %s <fig-bstretch1>`).
204 The force equations are given by the chain rule:
206 .. math:: \begin{array}{l}
207 \mathbf{F}_i ~=~ -\displaystyle\frac{d V_a({\theta_{ijk}})}{d \mathbf{r}_i} \\
208 \mathbf{F}_k ~=~ -\displaystyle\frac{d V_a({\theta_{ijk}})}{d \mathbf{r}_k} \\
209 \mathbf{F}_j ~=~ -\mathbf{F}_i-\mathbf{F}_k
211 ~ \mbox{ ~ where ~ } ~
212 {\theta_{ijk}}= \arccos \frac{(\mathbf{r}_{ij} \cdot \mathbf{r}_{kj})}{r_{ij}r_{kj}}
213 :label: eqnharmangleforce
215 The numbering :math:`i,j,k` is in sequence of covalently bonded atoms.
216 Atom :math:`j` is in the middle; atoms :math:`i` and :math:`k` are at
217 the ends (see :numref:`Fig. %s <fig-angle>`). **Note** that in the input in topology
218 files, angles are given in degrees and force constants in
219 kJ/mol/rad\ :math:`^2`.
223 Cosine based angle potential
224 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
226 In the GROMOS-96 force field a simplified function is used to represent
229 .. math:: V_a({\theta_{ijk}}) = {\frac{1}{2}}k^{\theta}_{ijk}\left(\cos({\theta_{ijk}}) - \cos({\theta_{ijk}}^0)\right)^2
234 .. math:: \cos({\theta_{ijk}}) = \frac{\mathbf{r}_{ij}\cdot\mathbf{r}_{kj}}{{r_{ij}}r_{kj}}
235 :label: eqnG96anglecos
237 The corresponding force can be derived by partial differentiation with
238 respect to the atomic positions. The force constants in this function
239 are related to the force constants in the harmonic form
240 :math:`k^{\theta,\mathrm{harm}}` (:ref:`harmonicangle`) by:
242 .. math:: k^{\theta} \sin^2({\theta_{ijk}}^0) = k^{\theta,\mathrm{harm}}
243 :label: eqnG96angleFC
245 In the GROMOS-96 manual there is a much more complicated conversion
246 formula which is temperature dependent. The formulas are equivalent at 0
247 K and the differences at 300 K are on the order of 0.1 to 0.2%. **Note**
248 that in the input in topology files, angles are given in degrees and
249 force constants in kJ/mol.
253 Restricted bending potential
254 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
256 The restricted bending (ReB) potential \ :ref:`83 <refMonicaGoga2013>` prevents the
257 bending angle :math:`\theta` from reaching the :math:`180^{\circ}`
258 value. In this way, the numerical instabilities due to the calculation
259 of the torsion angle and potential are eliminated when performing
260 coarse-grained molecular dynamics simulations.
262 To systematically hinder the bending angles from reaching the
263 :math:`180^{\circ}` value, the bending potential :eq:`eqn %s <eqnG96angle>` is
264 divided by a :math:`\sin^2\theta` factor:
266 .. math:: V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} \frac{(\cos\theta_i - \cos\theta_0)^2}{\sin^2\theta_i}.
269 :numref:`Figure %s <fig-ReB>` shows the comparison between the ReB potential,
270 :eq:`%s <eqnReB>`, and the standard one :eq:`%s <eqnG96angle>`.
274 .. figure:: plots/fig-02.*
277 Bending angle potentials: cosine harmonic (solid black line), angle
278 harmonic (dashed black line) and restricted bending (red) with the
279 same bending constant :math:`k_{\theta}=85` kJ mol\ :math:`^{-1}` and
280 equilibrium angle :math:`\theta_0=130^{\circ}`. The orange line
281 represents the sum of a cosine harmonic (:math:`k =50` kJ
282 mol\ :math:`^{-1}`) with a restricted bending (:math:`k =25` kJ
283 mol\ :math:`^{-1}`) potential, both with
284 :math:`\theta_0=130^{\circ}`.
286 The wall of the ReB potential is very repulsive in the region close to
287 :math:`180^{\circ}` and, as a result, the bending angles are kept within
288 a safe interval, far from instabilities. The power :math:`2` of
289 :math:`\sin\theta_i` in the denominator has been chosen to guarantee
290 this behavior and allows an elegant differentiation:
292 .. math:: F_{\rm ReB}(\theta_i) = \frac{2k_{\theta}}{\sin^4\theta_i}(\cos\theta_i - \cos\theta_0) (1 - \cos\theta_i\cos\theta_0) \frac{\partial \cos\theta_i}{\partial \vec r_{k}}.
295 Due to its construction, the restricted bending potential cannot be
296 used for equilibrium :math:`\theta_0` values too close to
297 :math:`0^{\circ}` or :math:`180^{\circ}` (from experience, at least
298 :math:`10^{\circ}` difference is recommended). It is very important
299 that, in the starting configuration, all the bending angles have to be
300 in the safe interval to avoid initial instabilities. This bending
301 potential can be used in combination with any form of torsion potential.
302 It will always prevent three consecutive particles from becoming
303 collinear and, as a result, any torsion potential will remain free of
304 singularities. It can be also added to a standard bending potential to
305 affect the angle around :math:`180^{\circ}`, but to keep its original
306 form around the minimum (see the orange curve in :numref:`Fig. %s <fig-ReB>`).
308 Urey-Bradley potential
309 ~~~~~~~~~~~~~~~~~~~~~~
311 The Urey-Bradley bond-angle vibration between a triplet of atoms
312 :math:`i` - :math:`j` - :math:`k` is represented by a harmonic potential
313 on the angle :math:`{\theta_{ijk}}` and a harmonic correction term on
314 the distance between the atoms :math:`i` and :math:`k`. Although this
315 can be easily written as a simple sum of two terms, it is convenient to
316 have it as a single entry in the topology file and in the output as a
317 separate energy term. It is used mainly in the CHARMm force
318 field \ :ref:`84 <refBBrooks83>`. The energy is given by:
320 .. math:: V_a({\theta_{ijk}}) = {\frac{1}{2}}k^{\theta}_{ijk}({\theta_{ijk}}-{\theta_{ijk}}^0)^2 + {\frac{1}{2}}k^{UB}_{ijk}(r_{ik}-r_{ik}^0)^2
323 The force equations can be deduced from sections :ref:`harmonicbond`
324 and :ref:`harmonicangle`.
329 The bond-bond cross term for three particles :math:`i, j, k` forming
330 bonds :math:`i-j` and :math:`k-j` is given
331 by \ :ref:`85 <refLawrence2003b>`:
333 .. math:: V_{rr'} ~=~ k_{rr'} \left(\left|\mathbf{r}_{i}-\mathbf{r}_j\right|-r_{1e}\right) \left(\left|\mathbf{r}_{k}-\mathbf{r}_j\right|-r_{2e}\right)
336 where :math:`k_{rr'}` is the force constant, and :math:`r_{1e}` and
337 :math:`r_{2e}` are the equilibrium bond lengths of the :math:`i-j` and
338 :math:`k-j` bonds respectively. The force associated with this potential
339 on particle :math:`i` is:
341 .. math:: \mathbf{F}_{i} = -k_{rr'}\left(\left|\mathbf{r}_{k}-\mathbf{r}_j\right|-r_{2e}\right)\frac{\mathbf{r}_i-\mathbf{r}_j}{\left|\mathbf{r}_{i}-\mathbf{r}_j\right|}
342 :label: eqncrossbbforce
344 The force on atom :math:`k` can be obtained by swapping :math:`i` and
345 :math:`k` in the above equation. Finally, the force on atom :math:`j`
346 follows from the fact that the sum of internal forces should be zero:
347 :math:`\mathbf{F}_j = -\mathbf{F}_i-\mathbf{F}_k`.
349 Bond-Angle cross term
350 ~~~~~~~~~~~~~~~~~~~~~
352 The bond-angle cross term for three particles :math:`i, j, k` forming
353 bonds :math:`i-j` and :math:`k-j` is given
354 by \ :ref:`85 <refLawrence2003b>`:
356 .. math:: V_{r\theta} ~=~ k_{r\theta} \left(\left|\mathbf{r}_{i}-\mathbf{r}_k\right|-r_{3e} \right) \left(\left|\mathbf{r}_{i}-\mathbf{r}_j\right|-r_{1e} + \left|\mathbf{r}_{k}-\mathbf{r}_j\right|-r_{2e}\right)
359 where :math:`k_{r\theta}` is the force constant, :math:`r_{3e}` is the
360 :math:`i-k` distance, and the other constants are the same as in
361 :eq:`Equation %s <eqncrossbb>`. The force associated with the potential on atom
364 .. math:: \mathbf{F}_{i} ~=~ -k_{r\theta}
367 \left| \mathbf{r}_{i} - \mathbf{r}_{k}\right|
370 \mathbf{r}_{i}-\mathbf{r}_j}
371 { \left| \mathbf{r}_{i}-\mathbf{r}_{j}\right|
374 \left| \mathbf{r}_{i}-\mathbf{r}_{j}\right|
376 + \left| \mathbf{r}_{k}-\mathbf{r}_{j}\right|
379 \mathbf{r}_{i}-\mathbf{r}_{k}}
380 {\left| \mathbf{r}_{i}-\mathbf{r}_{k}\right|
383 :label: eqncrossbaforce
385 Quartic angle potential
386 ~~~~~~~~~~~~~~~~~~~~~~~
388 For special purposes there is an angle potential that uses a fourth
391 .. math:: V_q({\theta_{ijk}}) ~=~ \sum_{n=0}^5 C_n ({\theta_{ijk}}-{\theta_{ijk}}^0)^n
392 :label: eqnquarticangle
399 Improper dihedrals are meant to keep planar groups (*e.g.* aromatic
400 rings) planar, or to prevent molecules from flipping over to their
401 mirror images, see :numref:`Fig. %s <fig-imp>`.
405 .. figure:: plots/ring-imp.*
408 Principle of improper dihedral angles. Out of plane bending for rings.
409 The improper dihedral angle :math:`\xi` is defined as the angle between
410 planes (i,j,k) and (j,k,l).
412 .. figure:: plots/subst-im.*
415 .. figure:: plots/tetra-im.*
418 Principle of improper dihedral angles. Out of tetrahedral angle.
419 The improper dihedral angle :math:`\xi` is defined
420 as the angle between planes (i,j,k) and (j,k,l).
422 Improper dihedrals: harmonic type
423 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
425 The simplest improper dihedral potential is a harmonic potential; it is
426 plotted in :numref:`Fig. %s <fig-imps>`.
428 .. math:: V_{id}(\xi_{ijkl}) = {\frac{1}{2}}k_{\xi}(\xi_{ijkl}-\xi_0)^2
429 :label: eqnharmimpdihedral
431 Since the potential is harmonic it is discontinuous, but since the
432 discontinuity is chosen at 180\ :math:`^\circ` distance from
433 :math:`\xi_0` this will never cause problems. **Note** that in the input
434 in topology files, angles are given in degrees and force constants in
435 kJ/mol/rad\ :math:`^2`.
439 .. figure:: plots/f-imps.*
442 Improper dihedral potential.
444 Improper dihedrals: periodic type
445 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
447 This potential is identical to the periodic proper dihedral (see below).
448 There is a separate dihedral type for this (type 4) only to be able to
449 distinguish improper from proper dihedrals in the parameter section and
455 For the normal dihedral interaction there is a choice of either the
456 GROMOS periodic function or a function based on expansion in powers of
457 :math:`\cos \phi` (the so-called Ryckaert-Bellemans potential). This
458 choice has consequences for the inclusion of special interactions
459 between the first and the fourth atom of the dihedral quadruple. With
460 the periodic GROMOS potential a special 1-4 LJ-interaction must be
461 included; with the Ryckaert-Bellemans potential *for alkanes* the 1-4
462 interactions must be excluded from the non-bonded list. **Note:**
463 Ryckaert-Bellemans potentials are also used in *e.g.* the OPLS force
464 field in combination with 1-4 interactions. You should therefore not
465 modify topologies generated by :ref:`pdb2gmx <gmx pdb2gmx>` in this case.
467 Proper dihedrals: periodic type
468 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
470 Proper dihedral angles are defined according to the IUPAC/IUB
471 convention, where :math:`\phi` is the angle between the :math:`ijk` and
472 the :math:`jkl` planes, with **zero** corresponding to the *cis*
473 configuration (:math:`i` and :math:`l` on the same side). There are two
474 dihedral function types in |Gromacs| topology files. There is the standard
475 type 1 which behaves like any other bonded interactions. For certain
476 force fields, type 9 is useful. Type 9 allows multiple potential
477 functions to be applied automatically to a single dihedral in the
478 ``[ dihedral ]`` section when multiple parameters are
479 defined for the same atomtypes in the ``[ dihedraltypes ]``
484 .. figure:: plots/f-dih.*
487 Principle of proper dihedral angle (left, in *trans* form) and the
488 dihedral angle potential (right).
490 .. math:: V_d(\phi_{ijkl}) = k_{\phi}(1 + \cos(n \phi - \phi_s))
491 :label: eqnperiodicpropdihedral
493 Proper dihedrals: Ryckaert-Bellemans function
494 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
496 | For alkanes, the following proper dihedral potential is often used
497 (see :numref:`Fig. %s <fig-rbdih>`):
499 .. math:: V_{rb}(\phi_{ijkl}) = \sum_{n=0}^5 C_n( \cos(\psi ))^n,
500 :label: eqnRBproperdihedral
502 | where :math:`\psi = \phi - 180^\circ`.
503 | **Note:** A conversion from one convention to another can be achieved
504 by multiplying every coefficient :math:`\displaystyle C_n` by
505 :math:`\displaystyle (-1)^n`.
507 An example of constants for :math:`C` is given in :numref:`Table %s <tab-crb>`.
512 Constants for Ryckaert-Bellemans potential (\ :math:`\mathrm{kJ mol}^{-1}`).
516 +-------------+-------+-------------+--------+-------------+-------+
517 | :math:`C_0` | 9.28 | :math:`C_2` | -13.12 | :math:`C_4` | 26.24 |
518 +-------------+-------+-------------+--------+-------------+-------+
519 | :math:`C_1` | 12.16 | :math:`C_3` | -3.06 | :math:`C_5` | -31.5 |
520 +-------------+-------+-------------+--------+-------------+-------+
525 .. figure:: plots/f-rbs.*
528 Ryckaert-Bellemans dihedral potential.
530 (**Note:** The use of this potential implies exclusion of LJ
531 interactions between the first and the last atom of the dihedral, and
532 :math:`\psi` is defined according to the “polymer convention”
533 (:math:`\psi_{trans}=0`).)
535 | The RB dihedral function can also be used to include Fourier dihedrals
538 .. math:: V_{rb} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[F_1(1+\cos(\phi)) + F_2(
539 1-\cos(2\phi)) + F_3(1+\cos(3\phi)) + F_4(1-\cos(4\phi))\right]
540 :label: eqnRBproperdihedralFourier
542 | Because of the equalities :math:`\cos(2\phi) = 2\cos^2(\phi) - 1`,
543 :math:`\cos(3\phi) = 4\cos^3(\phi) - 3\cos(\phi)` and
544 :math:`\cos(4\phi) = 8\cos^4(\phi) - 8\cos^2(\phi) + 1` one can
545 translate the OPLS parameters to Ryckaert-Bellemans parameters as
548 .. math:: \displaystyle
550 \displaystyle C_0&=&F_2 + \frac{1}{2} (F_1 + F_3)\\
551 \displaystyle C_1&=&\frac{1}{2} (- F_1 + 3 \, F_3)\\
552 \displaystyle C_2&=& -F_2 + 4 \, F_4\\
553 \displaystyle C_3&=&-2 \, F_3\\
554 \displaystyle C_4&=&-4 \, F_4\\
555 \displaystyle C_5&=&0
557 :label: eqnoplsRBconversion
559 | with OPLS parameters in protein convention and RB parameters in
560 polymer convention (this yields a minus sign for the odd powers of
561 cos\ :math:`(\phi)`).
562 | **Note:** Mind the conversion from **kcal mol**\ :math:`^{-1}` for
563 literature OPLS and RB parameters to **kJ mol**\ :math:`^{-1}` in
566 Proper dihedrals: Fourier function
567 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
569 | The OPLS potential function is given as the first three
570 :ref:`86 <refJorgensen1996>` or four \ :ref:`87 <refRobertson2015a>`
571 cosine terms of a Fourier series. In |Gromacs| the four term function is
574 .. math:: V_{F} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[C_1(1+\cos(\phi)) + C_2(
575 1-\cos(2\phi)) + C_3(1+\cos(3\phi)) + C_4(1-\cos(4\phi))\right],
576 :label: eqnfourierproperdihedral
578 | Internally, |Gromacs| uses the Ryckaert-Bellemans code to compute
579 Fourier dihedrals (see above), because this is more efficient.
580 | **Note:** Mind the conversion from *k*\ cal mol\ :math:`^{-1}` for
581 literature OPLS parameters to **kJ mol**\ :math:`^{-1}` in |Gromacs|.
583 Proper dihedrals: Restricted torsion potential
584 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
586 In a manner very similar to the restricted bending potential (see
587 :ref:`ReB`), a restricted torsion/dihedral potential is introduced:
589 .. math:: V_{\rm ReT}(\phi_i) = \frac{1}{2} k_{\phi} \frac{(\cos\phi_i - \cos\phi_0)^2}{\sin^2\phi_i}
592 with the advantages of being a function of :math:`\cos\phi` (no
593 problems taking the derivative of :math:`\sin\phi`) and of keeping the
594 torsion angle at only one minimum value. In this case, the factor
595 :math:`\sin^2\phi` does not allow the dihedral angle to move from the
596 [:math:`-180^{\circ}`:0] to [0::math:`180^{\circ}`] interval, i.e. it
597 cannot have maxima both at :math:`-\phi_0` and :math:`+\phi_0` maxima,
598 but only one of them. For this reason, all the dihedral angles of the
599 starting configuration should have their values in the desired angles
600 interval and the equilibrium :math:`\phi_0` value should not be too
601 close to the interval limits (as for the restricted bending potential,
602 described in :ref:`ReB`, at least :math:`10^{\circ}` difference is
605 Proper dihedrals: Combined bending-torsion potential
606 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
608 When the four particles forming the dihedral angle become collinear
609 (this situation will never happen in atomistic simulations, but it can
610 occur in coarse-grained simulations) the calculation of the torsion
611 angle and potential leads to numerical instabilities. One way to avoid
612 this is to use the restricted bending potential (see :ref:`ReB`) that
613 prevents the dihedral from reaching the :math:`180^{\circ}` value.
615 Another way is to disregard any effects of the dihedral becoming
616 ill-defined, keeping the dihedral force and potential calculation
617 continuous in entire angle range by coupling the torsion potential (in a
618 cosine form) with the bending potentials of the adjacent bending angles
619 in a unique expression:
621 .. math:: V_{\rm CBT}(\theta_{i-1}, \theta_i, \phi_i) = k_{\phi} \sin^3\theta_{i-1} \sin^3\theta_{i} \sum_{n=0}^4 { a_n \cos^n\phi_i}.
624 This combined bending-torsion (CBT) potential has been proposed
625 by \ :ref:`88 <refBulacuGiessen2005>` for polymer melt simulations and is
626 extensively described in \ :ref:`83 <refMonicaGoga2013>`.
628 This potential has two main advantages:
630 - it does not only depend on the dihedral angle :math:`\phi_i` (between
631 the :math:`i-2`, :math:`i-1`, :math:`i` and :math:`i+1` beads) but
632 also on the bending angles :math:`\theta_{i-1}` and :math:`\theta_i`
633 defined from three adjacent beads (:math:`i-2`, :math:`i-1` and
634 :math:`i`, and :math:`i-1`, :math:`i` and :math:`i+1`, respectively).
635 The two :math:`\sin^3\theta` pre-factors, tentatively suggested
636 by \ :ref:`89 <refScottScheragator1966>` and theoretically discussed by
637 \ :ref:`90 <refPaulingBond>`, cancel the torsion potential and force when either of the two
638 bending angles approaches the value of :math:`180^\circ`.
640 - its dependence on :math:`\phi_i` is expressed through a polynomial in
641 :math:`\cos\phi_i` that avoids the singularities in
642 :math:`\phi=0^\circ` or :math:`180^\circ` in calculating the
645 These two properties make the CBT potential well-behaved for MD
646 simulations with weak constraints on the bending angles or even for
647 steered / non-equilibrium MD in which the bending and torsion angles
648 suffer major modifications. When using the CBT potential, the bending
649 potentials for the adjacent :math:`\theta_{i-1}` and :math:`\theta_i`
650 may have any form. It is also possible to leave out the two angle
651 bending terms (:math:`\theta_{i-1}` and :math:`\theta_{i}`) completely.
652 :numref:`Fig. %s <fig-CBT>` illustrates the difference between a torsion potential
653 with and without the :math:`\sin^{3}\theta` factors (blue and gray
654 curves, respectively).
658 .. figure:: plots/fig-04.*
661 Blue: surface plot of the combined bending-torsion potential
662 (:eq:`%s <eqnCBT>` with :math:`k = 10` kJ mol\ :math:`^{-1}`,
663 :math:`a_0=2.41`, :math:`a_1=-2.95`, :math:`a_2=0.36`,
664 :math:`a_3=1.33`) when, for simplicity, the bending angles behave the
665 same (:math:`\theta_1=\theta_2=\theta`). Gray: the same torsion
666 potential without the :math:`\sin^{3}\theta` terms
667 (Ryckaert-Bellemans type). :math:`\phi` is the dihedral angle.
669 Additionally, the derivative of :math:`V_{CBT}` with respect to the
670 Cartesian variables is straightforward:
672 .. math:: \frac{\partial V_{\rm CBT}(\theta_{i-1},\theta_i,\phi_i)} {\partial \vec r_{l}} = \frac{\partial V_{\rm CBT}}{\partial \theta_{i-1}} \frac{\partial \theta_{i-1}}{\partial \vec r_{l}} +
673 \frac{\partial V_{\rm CBT}}{\partial \theta_{i }} \frac{\partial \theta_{i }}{\partial \vec r_{l}} +
674 \frac{\partial V_{\rm CBT}}{\partial \phi_{i }} \frac{\partial \phi_{i }}{\partial \vec r_{l}}
677 The CBT is based on a cosine form without multiplicity, so it can only
678 be symmetrical around :math:`0^{\circ}`. To obtain an asymmetrical
679 dihedral angle distribution (e.g. only one maximum in
680 [:math:`-180^{\circ}`::math:`180^{\circ}`] interval), a standard torsion
681 potential such as harmonic angle or periodic cosine potentials should be
682 used instead of a CBT potential. However, these two forms have the
683 inconveniences of the force derivation (:math:`1/\sin\phi`) and of the
684 alignment of beads (:math:`\theta_i` or
685 :math:`\theta_{i-1} = 0^{\circ}, 180^{\circ}`). Coupling such
686 non-\ :math:`\cos\phi` potentials with :math:`\sin^3\theta` factors does
687 not improve simulation stability since there are cases in which
688 :math:`\theta` and :math:`\phi` are simultaneously :math:`180^{\circ}`.
689 The integration at this step would be possible (due to the cancelling of
690 the torsion potential) but the next step would be singular
691 (:math:`\theta` is not :math:`180^{\circ}` and :math:`\phi` is very
692 close to :math:`180^{\circ}`).
694 Tabulated bonded interaction functions
695 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
697 | For full flexibility, any functional shape can be used for bonds,
698 angles and dihedrals through user-supplied tabulated functions. The
699 functional shapes are:
701 .. math:: \begin{aligned}
702 V_b(r_{ij}) &=& k \, f^b_n(r_{ij}) \\
703 V_a({\theta_{ijk}}) &=& k \, f^a_n({\theta_{ijk}}) \\
704 V_d(\phi_{ijkl}) &=& k \, f^d_n(\phi_{ijkl})\end{aligned}
705 :label: eqntabuöatedbond
707 | where :math:`k` is a force constant in units of energy and :math:`f`
708 is a cubic spline function; for details see :ref:`cubicspline`. For
709 each interaction, the force constant :math:`k` and the table number
710 :math:`n` are specified in the topology. There are two different types
711 of bonds, one that generates exclusions (type 8) and one that does not
712 (type 9). For details see :numref:`Table %s <tab-topfile2>`. The table files are
713 supplied to the :ref:`mdrun <gmx mdrun>` program. After the table file name an
714 underscore, the letter “b” for bonds, “a” for angles or “d” for
715 dihedrals and the table number must be appended. For example, a
716 tabulated bond with :math:`n=0` can be read from the file
717 table\_b0.xvg. Multiple tables can be supplied simply by adding files
718 with different values of :math:`n`, and are applied to the appropriate
719 bonds, as specified in the topology (:numref:`Table %s <tab-topfile2>`). The format
720 for the table files is three fixed-format columns of any suitable
721 width. These columns must contain :math:`x`, :math:`f(x)`,
722 :math:`-f'(x)`, and the values of :math:`x` should be uniformly
723 spaced. Requirements for entries in the topology are given
724 in :numref:`Table %s <tab-topfile2>`. The setup of the tables is as follows:
725 | **bonds**: :math:`x` is the distance in nm. For distances beyond the
726 table length, :ref:`mdrun <gmx mdrun>` will quit with an error message.
727 | **angles**: :math:`x` is the angle in degrees. The table should go
728 from 0 up to and including 180 degrees; the derivative is taken in
730 | **dihedrals**: :math:`x` is the dihedral angle in degrees. The table
731 should go from -180 up to and including 180 degrees; the IUPAC/IUB
732 convention is used, *i.e.* zero is cis, the derivative is taken in