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
35 See also :numref:`Fig. %s <fig-bstretch1>`, with the force given by:
37 .. math:: \mathbf{F}_i(\mathbf{r}_ij) = k^b_{ij}({r_{ij}}-b_{ij}) {\frac{{\mathbf{r}_{ij}}}{{r_{ij}}}}
41 Fourth power potential
42 ^^^^^^^^^^^^^^^^^^^^^^
44 In the GROMOS-96 force field \ :ref:`77 <refgromos96>`, the covalent bond
45 potential is, for reasons of computational efficiency, written as:
47 .. math:: V_b~({r_{ij}}) = \frac{1}{4}k^b_{ij}\left({r_{ij}}^2-b_{ij}^2\right)^2
49 The corresponding force is:
51 .. math:: \mathbf{F}_i(\mathbf{r}_ij) = k^b_{ij}({r_{ij}}^2-b_{ij}^2)~\mathbf{r}_ij
53 The force constants for this form of the potential are related to the
54 usual harmonic force constant :math:`k^{b,\mathrm{harm}}`
55 (sec. :ref:`bondpot`) as
57 .. math:: 2 k^b b_{ij}^2 = k^{b,\mathrm{harm}}
59 The force constants are mostly derived from the harmonic ones used in
60 GROMOS-87 :ref:`78 <refbiomos>`. Although this form is
61 computationally more efficient (because no square root has to be
62 evaluated), it is conceptually more complex. One particular disadvantage
63 is that since the form is not harmonic, the average energy of a single
64 bond is not equal to :math:`{\frac{1}{2}}kT` as it is for the normal
69 Morse potential bond stretching
70 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
72 For some systems that require an anharmonic bond stretching potential,
73 the Morse potential \ :ref:`79 <refMorse29>` between two atoms *i* and *j* is
74 available in |Gromacs|. This potential differs from the harmonic potential
75 in that it has an asymmetric potential well and a zero force at infinite
76 distance. The functional form is:
78 .. math:: \displaystyle V_{morse} (r_{ij}) = D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2,
80 See also :numref:`Fig. %s <fig-morse>`, and the corresponding force is:
85 \displaystyle {\bf F}_{morse} ({\bf r}_{ij})&=&2 D_{ij} \beta_{ij} \exp(-\beta_{ij}(r_{ij}-b_{ij})) * \\
86 \displaystyle \: & \: &[1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))] \frac{\displaystyle {\bf r}_{ij}}{\displaystyle r_{ij}},
89 where :math:`\displaystyle D_{ij}` is the depth of the well in
90 kJ/mol, :math:`\displaystyle \beta_{ij}` defines the steepness of the
91 well (in nm\ :math:`^{-1}`), and :math:`\displaystyle b_{ij}` is the
92 equilibrium distance in nm. The steepness parameter
93 :math:`\displaystyle \beta_{ij}` can be expressed in terms of the reduced mass of the atoms *i* and
94 *j*, the fundamental vibration frequency :math:`\displaystyle\omega_{ij}` and the well depth :math:`\displaystyle D_{ij}`:
96 .. math:: \displaystyle \beta_{ij}= \omega_{ij} \sqrt{\frac{\mu_{ij}}{2 D_{ij}}}
98 and because :math:`\displaystyle \omega = \sqrt{k/\mu}`, one can
99 rewrite :math:`\displaystyle \beta_{ij}` in terms of the harmonic
100 force constant :math:`\displaystyle k_{ij}`:
102 .. math:: \displaystyle \beta_{ij}= \sqrt{\frac{k_{ij}}{2 D_{ij}}}
105 For small deviations :math:`\displaystyle (r_{ij}-b_{ij})`, one can
106 approximate the :math:`\displaystyle \exp`-term to first-order using a
109 .. math:: \displaystyle \exp(-x) \approx 1-x
112 and substituting :eq:`eqn. %s <eqnbetaij>` and :eq:`eqn. %s <eqnexpminx>` in the
118 \displaystyle V_{morse} (r_{ij})&=&D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2\\
119 \displaystyle \:&=&D_{ij} [1 - (1 -\sqrt{\frac{k_{ij}}{2 D_{ij}}}(r_{ij}-b_{ij}))]^2\\
120 \displaystyle \:&=&\frac{1}{2} k_{ij} (r_{ij}-b_{ij}))^2
123 we recover the harmonic bond stretching potential.
127 .. figure:: plots/f-morse.*
130 The Morse potential well, with bond length 0.15 nm.
132 Cubic bond stretching potential
133 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
135 Another anharmonic bond stretching potential that is slightly simpler
136 than the Morse potential adds a cubic term in the distance to the simple
139 .. math:: V_b~({r_{ij}}) = k^b_{ij}({r_{ij}}-b_{ij})^2 + k^b_{ij}k^{cub}_{ij}({r_{ij}}-b_{ij})^3
141 A flexible water model (based on the SPC water model \ :ref:`80 <refBerendsen81>`)
142 including a cubic bond stretching potential for the O-H bond was
143 developed by Ferguson \ :ref:`81 <refFerguson95>`. This model was found to yield a
144 reasonable infrared spectrum. The Ferguson water model is available in
145 the |Gromacs| library (``flexwat-ferguson.itp``). It should be noted that the
146 potential is asymmetric: overstretching leads to infinitely low
147 energies. The integration timestep is therefore limited to 1 fs.
149 The force corresponding to this potential is:
151 .. 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}}}}
153 FENE bond stretching potential
154 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
156 In coarse-grained polymer simulations the beads are often connected by a
157 FENE (finitely extensible nonlinear elastic) potential \ :ref:`82 <refWarner72>`:
161 V_{\mbox{\small FENE}}({r_{ij}}) =
162 -{\frac{1}{2}}k^b_{ij} b^2_{ij} \log\left(1 - \frac{{r_{ij}}^2}{b^2_{ij}}\right)
164 The potential looks complicated, but the expression for the force is
169 F_{\mbox{\small FENE}}(\mathbf{r}_ij) =
170 -k^b_{ij} \left(1 - \frac{{r_{ij}}^2}{b^2_{ij}}\right)^{-1} \mathbf{r}_ij
172 At short distances the potential asymptotically goes to a harmonic
173 potential with force constant :math:`k^b`, while it diverges at distance
178 Harmonic angle potential
179 ~~~~~~~~~~~~~~~~~~~~~~~~
181 The bond-angle vibration between a triplet of atoms :math:`i` -
182 :math:`j` - :math:`k` is also represented by a harmonic potential on the
183 angle :math:`{\theta_{ijk}}`
187 .. figure:: plots/angle.*
190 Principle of angle vibration (left) and the bond angle potential
193 .. math:: V_a({\theta_{ijk}}) = {\frac{1}{2}}k^{\theta}_{ijk}({\theta_{ijk}}-{\theta_{ijk}}^0)^2
195 As the bond-angle vibration is represented by a harmonic potential, the
196 form is the same as the bond stretching
197 (:numref:`Fig. %s <fig-bstretch1>`).
199 The force equations are given by the chain rule:
204 \mathbf{F}_i ~=~ -\displaystyle\frac{d V_a({\theta_{ijk}})}{d \mathbf{r}_i} \\
205 \mathbf{F}_k ~=~ -\displaystyle\frac{d V_a({\theta_{ijk}})}{d \mathbf{r}_k} \\
206 \mathbf{F}_j ~=~ -\mathbf{F}_i-\mathbf{F}_k
208 ~ \mbox{ ~ where ~ } ~
209 {\theta_{ijk}}= \arccos \frac{(\mathbf{r}_ij \cdot \mathbf{r}_{kj})}{r_{ij}r_{kj}}
211 The numbering :math:`i,j,k` is in sequence of covalently bonded atoms.
212 Atom :math:`j` is in the middle; atoms :math:`i` and :math:`k` are at
213 the ends (see :numref:`Fig. %s <fig-angle>`). **Note** that in the input in topology
214 files, angles are given in degrees and force constants in
215 kJ/mol/rad\ :math:`^2`.
219 Cosine based angle potential
220 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
222 In the GROMOS-96 force field a simplified function is used to represent
225 .. math:: V_a({\theta_{ijk}}) = {\frac{1}{2}}k^{\theta}_{ijk}\left(\cos({\theta_{ijk}}) - \cos({\theta_{ijk}}^0)\right)^2
230 .. math:: \cos({\theta_{ijk}}) = \frac{\mathbf{r}_ij\cdot\mathbf{r}_{kj}}{{r_{ij}}r_{kj}}
232 The corresponding force can be derived by partial differentiation with
233 respect to the atomic positions. The force constants in this function
234 are related to the force constants in the harmonic form
235 :math:`k^{\theta,\mathrm{harm}}` (:ref:`harmonicangle`) by:
237 .. math:: k^{\theta} \sin^2({\theta_{ijk}}^0) = k^{\theta,\mathrm{harm}}
239 In the GROMOS-96 manual there is a much more complicated conversion
240 formula which is temperature dependent. The formulas are equivalent at 0
241 K and the differences at 300 K are on the order of 0.1 to 0.2%. **Note**
242 that in the input in topology files, angles are given in degrees and
243 force constants in kJ/mol.
247 Restricted bending potential
248 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
250 The restricted bending (ReB) potential \ :ref:`83 <refMonicaGoga2013>` prevents the
251 bending angle :math:`\theta` from reaching the :math:`180^{\circ}`
252 value. In this way, the numerical instabilities due to the calculation
253 of the torsion angle and potential are eliminated when performing
254 coarse-grained molecular dynamics simulations.
256 To systematically hinder the bending angles from reaching the
257 :math:`180^{\circ}` value, the bending potential :eq:`eqn %s <eqG96angle>` is
258 divided by a :math:`\sin^2\theta` factor:
260 .. math:: V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} \frac{(\cos\theta_i - \cos\theta_0)^2}{\sin^2\theta_i}.
263 :numref:`Figure %s <fig-ReB>` shows the comparison between the ReB potential,
264 :eq:`%s <eqReB>`, and the standard one :eq:`%s <eqG96angle>`.
268 .. figure:: plots/fig-02.*
271 Bending angle potentials: cosine harmonic (solid black line), angle
272 harmonic (dashed black line) and restricted bending (red) with the
273 same bending constant :math:`k_{\theta}=85` kJ mol\ :math:`^{-1}` and
274 equilibrium angle :math:`\theta_0=130^{\circ}`. The orange line
275 represents the sum of a cosine harmonic (:math:`k =50` kJ
276 mol\ :math:`^{-1}`) with a restricted bending (:math:`k =25` kJ
277 mol\ :math:`^{-1}`) potential, both with
278 :math:`\theta_0=130^{\circ}`.
280 The wall of the ReB potential is very repulsive in the region close to
281 :math:`180^{\circ}` and, as a result, the bending angles are kept within
282 a safe interval, far from instabilities. The power :math:`2` of
283 :math:`\sin\theta_i` in the denominator has been chosen to guarantee
284 this behavior and allows an elegant differentiation:
286 .. 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}}.
289 Due to its construction, the restricted bending potential cannot be
290 used for equilibrium :math:`\theta_0` values too close to
291 :math:`0^{\circ}` or :math:`180^{\circ}` (from experience, at least
292 :math:`10^{\circ}` difference is recommended). It is very important
293 that, in the starting configuration, all the bending angles have to be
294 in the safe interval to avoid initial instabilities. This bending
295 potential can be used in combination with any form of torsion potential.
296 It will always prevent three consecutive particles from becoming
297 collinear and, as a result, any torsion potential will remain free of
298 singularities. It can be also added to a standard bending potential to
299 affect the angle around :math:`180^{\circ}`, but to keep its original
300 form around the minimum (see the orange curve in :numref:`Fig. %s <fig-ReB>`).
302 Urey-Bradley potential
303 ~~~~~~~~~~~~~~~~~~~~~~
305 The Urey-Bradley bond-angle vibration between a triplet of atoms
306 :math:`i` - :math:`j` - :math:`k` is represented by a harmonic potential
307 on the angle :math:`{\theta_{ijk}}` and a harmonic correction term on
308 the distance between the atoms :math:`i` and :math:`k`. Although this
309 can be easily written as a simple sum of two terms, it is convenient to
310 have it as a single entry in the topology file and in the output as a
311 separate energy term. It is used mainly in the CHARMm force
312 field \ :ref:`84 <refBBrooks83>`. The energy is given by:
314 .. 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
316 The force equations can be deduced from sections :ref:`harmonicbond`
317 and :ref:`harmonicangle`.
322 The bond-bond cross term for three particles :math:`i, j, k` forming
323 bonds :math:`i-j` and :math:`k-j` is given
324 by \ :ref:`85 <refLawrence2003b>`:
326 .. 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)
329 where :math:`k_{rr'}` is the force constant, and :math:`r_{1e}` and
330 :math:`r_{2e}` are the equilibrium bond lengths of the :math:`i-j` and
331 :math:`k-j` bonds respectively. The force associated with this potential
332 on particle :math:`i` is:
334 .. 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|}
336 The force on atom :math:`k` can be obtained by swapping :math:`i` and
337 :math:`k` in the above equation. Finally, the force on atom :math:`j`
338 follows from the fact that the sum of internal forces should be zero:
339 :math:`\mathbf{F}_j = -\mathbf{F}_i-\mathbf{F}_k`.
341 Bond-Angle cross term
342 ~~~~~~~~~~~~~~~~~~~~~
344 The bond-angle cross term for three particles :math:`i, j, k` forming
345 bonds :math:`i-j` and :math:`k-j` is given
346 by \ :ref:`85 <refLawrence2003b>`:
348 .. 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)
350 where :math:`k_{r\theta}` is the force constant, :math:`r_{3e}` is the
351 :math:`i-k` distance, and the other constants are the same as in
352 :eq:`Equation %s <eqncrossbb>`. The force associated with the potential on atom
357 \mathbf{F}_{i} ~=~ -k_{r\theta}
360 \left| \mathbf{r}_{i} - \mathbf{r}_{k}\right|
363 \mathbf{r}_{i}-\mathbf{r}_j}
364 { \left| \mathbf{r}_{i}-\mathbf{r}_{j}\right|
367 \left| \mathbf{r}_{i}-\mathbf{r}_{j}\right|
369 + \left| \mathbf{r}_{k}-\mathbf{r}_{j}\right|
372 \mathbf{r}_{i}-\mathbf{r}_{k}}
373 {\left| \mathbf{r}_{i}-\mathbf{r}_{k}\right|
377 Quartic angle potential
378 ~~~~~~~~~~~~~~~~~~~~~~~
380 For special purposes there is an angle potential that uses a fourth
383 .. math:: V_q({\theta_{ijk}}) ~=~ \sum_{n=0}^5 C_n ({\theta_{ijk}}-{\theta_{ijk}}^0)^n
390 Improper dihedrals are meant to keep planar groups (*e.g.* aromatic
391 rings) planar, or to prevent molecules from flipping over to their
392 mirror images, see :numref:`Fig. %s <fig-imp>`.
396 .. figure:: plots/ring-imp.*
399 Principle of improper dihedral angles. Out of plane bending for rings.
400 The improper dihedral angle :math:`\xi` is defined as the angle between
401 planes (i,j,k) and (j,k,l).
403 .. figure:: plots/subst-im.*
406 .. figure:: plots/tetra-im.*
409 Principle of improper dihedral angles. Out of tetrahedral angle.
410 The improper dihedral angle :math:`\xi` is defined
411 as the angle between planes (i,j,k) and (j,k,l).
413 Improper dihedrals: harmonic type
414 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
416 The simplest improper dihedral potential is a harmonic potential; it is
417 plotted in :numref:`Fig. %s <fig-imps>`.
419 .. math:: V_{id}(\xi_{ijkl}) = {\frac{1}{2}}k_{\xi}(\xi_{ijkl}-\xi_0)^2
421 Since the potential is harmonic it is discontinuous, but since the
422 discontinuity is chosen at 180\ :math:`^\circ` distance from
423 :math:`\xi_0` this will never cause problems. **Note** that in the input
424 in topology files, angles are given in degrees and force constants in
425 kJ/mol/rad\ :math:`^2`.
429 .. figure:: plots/f-imps.*
432 Improper dihedral potential.
434 Improper dihedrals: periodic type
435 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
437 This potential is identical to the periodic proper dihedral (see below).
438 There is a separate dihedral type for this (type 4) only to be able to
439 distinguish improper from proper dihedrals in the parameter section and
445 For the normal dihedral interaction there is a choice of either the
446 GROMOS periodic function or a function based on expansion in powers of
447 :math:`\cos \phi` (the so-called Ryckaert-Bellemans potential). This
448 choice has consequences for the inclusion of special interactions
449 between the first and the fourth atom of the dihedral quadruple. With
450 the periodic GROMOS potential a special 1-4 LJ-interaction must be
451 included; with the Ryckaert-Bellemans potential *for alkanes* the 1-4
452 interactions must be excluded from the non-bonded list. **Note:**
453 Ryckaert-Bellemans potentials are also used in *e.g.* the OPLS force
454 field in combination with 1-4 interactions. You should therefore not
455 modify topologies generated by :ref:`pdb2gmx <gmx pdb2gmx>` in this case.
457 Proper dihedrals: periodic type
458 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
460 Proper dihedral angles are defined according to the IUPAC/IUB
461 convention, where :math:`\phi` is the angle between the :math:`ijk` and
462 the :math:`jkl` planes, with **zero** corresponding to the *cis*
463 configuration (:math:`i` and :math:`l` on the same side). There are two
464 dihedral function types in |Gromacs| topology files. There is the standard
465 type 1 which behaves like any other bonded interactions. For certain
466 force fields, type 9 is useful. Type 9 allows multiple potential
467 functions to be applied automatically to a single dihedral in the
468 ``[ dihedral ]`` section when multiple parameters are
469 defined for the same atomtypes in the ``[ dihedraltypes ]``
474 .. figure:: plots/f-dih.*
477 Principle of proper dihedral angle (left, in *trans* form) and the
478 dihedral angle potential (right).
480 .. math:: V_d(\phi_{ijkl}) = k_{\phi}(1 + \cos(n \phi - \phi_s))
482 Proper dihedrals: Ryckaert-Bellemans function
483 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
485 | For alkanes, the following proper dihedral potential is often used
486 (see :numref:`Fig. %s <fig-rbdih>`):
488 .. math:: V_{rb}(\phi_{ijkl}) = \sum_{n=0}^5 C_n( \cos(\psi ))^n,
490 | where :math:`\psi = \phi - 180^\circ`.
491 | **Note:** A conversion from one convention to another can be achieved
492 by multiplying every coefficient :math:`\displaystyle C_n` by
493 :math:`\displaystyle (-1)^n`.
495 An example of constants for :math:`C` is given in :numref:`Table %s <tab-crb>`.
500 Constants for Ryckaert-Bellemans potential (\ :math:`\mathrm{kJ mol}^{-1}`).
504 +-------------+-------+-------------+--------+-------------+-------+
505 | :math:`C_0` | 9.28 | :math:`C_2` | -13.12 | :math:`C_4` | 26.24 |
506 +-------------+-------+-------------+--------+-------------+-------+
507 | :math:`C_1` | 12.16 | :math:`C_3` | -3.06 | :math:`C_5` | -31.5 |
508 +-------------+-------+-------------+--------+-------------+-------+
513 .. figure:: plots/f-rbs.*
516 Ryckaert-Bellemans dihedral potential.
518 (**Note:** The use of this potential implies exclusion of LJ
519 interactions between the first and the last atom of the dihedral, and
520 :math:`\psi` is defined according to the “polymer convention”
521 (:math:`\psi_{trans}=0`).)
523 | The RB dihedral function can also be used to include Fourier dihedrals
528 V_{rb} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[F_1(1+\cos(\phi)) + F_2(
529 1-\cos(2\phi)) + F_3(1+\cos(3\phi)) + F_4(1-\cos(4\phi))\right]
531 | Because of the equalities :math:`\cos(2\phi) = 2\cos^2(\phi) - 1`,
532 :math:`\cos(3\phi) = 4\cos^3(\phi) - 3\cos(\phi)` and
533 :math:`\cos(4\phi) = 8\cos^4(\phi) - 8\cos^2(\phi) + 1` one can
534 translate the OPLS parameters to Ryckaert-Bellemans parameters as
541 \displaystyle C_0&=&F_2 + \frac{1}{2} (F_1 + F_3)\\
542 \displaystyle C_1&=&\frac{1}{2} (- F_1 + 3 \, F_3)\\
543 \displaystyle C_2&=& -F_2 + 4 \, F_4\\
544 \displaystyle C_3&=&-2 \, F_3\\
545 \displaystyle C_4&=&-4 \, F_4\\
546 \displaystyle C_5&=&0
549 | with OPLS parameters in protein convention and RB parameters in
550 polymer convention (this yields a minus sign for the odd powers of
551 cos\ :math:`(\phi)`).
552 | **Note:** Mind the conversion from **kcal mol**\ :math:`^{-1}` for
553 literature OPLS and RB parameters to **kJ mol**\ :math:`^{-1}` in
556 Proper dihedrals: Fourier function
557 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
559 | The OPLS potential function is given as the first three
560 :ref:`86 <refJorgensen1996>` or four \ :ref:`87 <refRobertson2015a>`
561 cosine terms of a Fourier series. In |Gromacs| the four term function is
566 V_{F} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[C_1(1+\cos(\phi)) + C_2(
567 1-\cos(2\phi)) + C_3(1+\cos(3\phi)) + C_4(1-\cos(4\phi))\right],
569 | Internally, |Gromacs| uses the Ryckaert-Bellemans code to compute
570 Fourier dihedrals (see above), because this is more efficient.
571 | **Note:** Mind the conversion from *k*\ cal mol\ :math:`^{-1}` for
572 literature OPLS parameters to **kJ mol**\ :math:`^{-1}` in |Gromacs|.
574 Proper dihedrals: Restricted torsion potential
575 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
577 In a manner very similar to the restricted bending potential (see
578 :ref:`ReB`), a restricted torsion/dihedral potential is introduced:
580 .. math:: V_{\rm ReT}(\phi_i) = \frac{1}{2} k_{\phi} \frac{(\cos\phi_i - \cos\phi_0)^2}{\sin^2\phi_i}
583 with the advantages of being a function of :math:`\cos\phi` (no
584 problems taking the derivative of :math:`\sin\phi`) and of keeping the
585 torsion angle at only one minimum value. In this case, the factor
586 :math:`\sin^2\phi` does not allow the dihedral angle to move from the
587 [:math:`-180^{\circ}`:0] to [0::math:`180^{\circ}`] interval, i.e. it
588 cannot have maxima both at :math:`-\phi_0` and :math:`+\phi_0` maxima,
589 but only one of them. For this reason, all the dihedral angles of the
590 starting configuration should have their values in the desired angles
591 interval and the the equilibrium :math:`\phi_0` value should not be too
592 close to the interval limits (as for the restricted bending potential,
593 described in :ref:`ReB`, at least :math:`10^{\circ}` difference is
596 Proper dihedrals: Combined bending-torsion potential
597 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
599 When the four particles forming the dihedral angle become collinear
600 (this situation will never happen in atomistic simulations, but it can
601 occur in coarse-grained simulations) the calculation of the torsion
602 angle and potential leads to numerical instabilities. One way to avoid
603 this is to use the restricted bending potential (see :ref:`ReB`) that
604 prevents the dihedral from reaching the :math:`180^{\circ}` value.
606 Another way is to disregard any effects of the dihedral becoming
607 ill-defined, keeping the dihedral force and potential calculation
608 continuous in entire angle range by coupling the torsion potential (in a
609 cosine form) with the bending potentials of the adjacent bending angles
610 in a unique expression:
612 .. 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}.
615 This combined bending-torsion (CBT) potential has been proposed
616 by \ :ref:`88 <refBulacuGiessen2005>` for polymer melt simulations and is
617 extensively described in \ :ref:`83 <refMonicaGoga2013>`.
619 This potential has two main advantages:
621 - it does not only depend on the dihedral angle :math:`\phi_i` (between
622 the :math:`i-2`, :math:`i-1`, :math:`i` and :math:`i+1` beads) but
623 also on the bending angles :math:`\theta_{i-1}` and :math:`\theta_i`
624 defined from three adjacent beads (:math:`i-2`, :math:`i-1` and
625 :math:`i`, and :math:`i-1`, :math:`i` and :math:`i+1`, respectively).
626 The two :math:`\sin^3\theta` pre-factors, tentatively suggested
627 by \ :ref:`89 <refScottScheragator1966>` and theoretically discussed by
628 \ :ref:`90 <refPaulingBond>`, cancel the torsion potential and force when either of the two
629 bending angles approaches the value of :math:`180^\circ`.
631 - its dependence on :math:`\phi_i` is expressed through a polynomial in
632 :math:`\cos\phi_i` that avoids the singularities in
633 :math:`\phi=0^\circ` or :math:`180^\circ` in calculating the
636 These two properties make the CBT potential well-behaved for MD
637 simulations with weak constraints on the bending angles or even for
638 steered / non-equilibrium MD in which the bending and torsion angles
639 suffer major modifications. When using the CBT potential, the bending
640 potentials for the adjacent :math:`\theta_{i-1}` and :math:`\theta_i`
641 may have any form. It is also possible to leave out the two angle
642 bending terms (:math:`\theta_{i-1}` and :math:`\theta_{i}`) completely.
643 :numref:`Fig. %s <fig-CBT>` illustrates the difference between a torsion potential
644 with and without the :math:`\sin^{3}\theta` factors (blue and gray
645 curves, respectively).
649 .. figure:: plots/fig-04.*
652 Blue: surface plot of the combined bending-torsion potential
653 (:eq:`%s <eqCBT>` with :math:`k = 10` kJ mol\ :math:`^{-1}`,
654 :math:`a_0=2.41`, :math:`a_1=-2.95`, :math:`a_2=0.36`,
655 :math:`a_3=1.33`) when, for simplicity, the bending angles behave the
656 same (:math:`\theta_1=\theta_2=\theta`). Gray: the same torsion
657 potential without the :math:`\sin^{3}\theta` terms
658 (Ryckaert-Bellemans type). :math:`\phi` is the dihedral angle.
660 Additionally, the derivative of :math:`V_{CBT}` with respect to the
661 Cartesian variables is straightforward:
663 .. 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}} +
664 \frac{\partial V_{\rm CBT}}{\partial \theta_{i }} \frac{\partial \theta_{i }}{\partial \vec r_{l}} +
665 \frac{\partial V_{\rm CBT}}{\partial \phi_{i }} \frac{\partial \phi_{i }}{\partial \vec r_{l}}
668 The CBT is based on a cosine form without multiplicity, so it can only
669 be symmetrical around :math:`0^{\circ}`. To obtain an asymmetrical
670 dihedral angle distribution (e.g. only one maximum in
671 [:math:`-180^{\circ}`::math:`180^{\circ}`] interval), a standard torsion
672 potential such as harmonic angle or periodic cosine potentials should be
673 used instead of a CBT potential. However, these two forms have the
674 inconveniences of the force derivation (:math:`1/\sin\phi`) and of the
675 alignment of beads (:math:`\theta_i` or
676 :math:`\theta_{i-1} = 0^{\circ}, 180^{\circ}`). Coupling such
677 non-\ :math:`\cos\phi` potentials with :math:`\sin^3\theta` factors does
678 not improve simulation stability since there are cases in which
679 :math:`\theta` and :math:`\phi` are simultaneously :math:`180^{\circ}`.
680 The integration at this step would be possible (due to the cancelling of
681 the torsion potential) but the next step would be singular
682 (:math:`\theta` is not :math:`180^{\circ}` and :math:`\phi` is very
683 close to :math:`180^{\circ}`).
685 Tabulated bonded interaction functions
686 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
688 | For full flexibility, any functional shape can be used for bonds,
689 angles and dihedrals through user-supplied tabulated functions. The
690 functional shapes are:
695 V_b(r_{ij}) &=& k \, f^b_n(r_{ij}) \\
696 V_a({\theta_{ijk}}) &=& k \, f^a_n({\theta_{ijk}}) \\
697 V_d(\phi_{ijkl}) &=& k \, f^d_n(\phi_{ijkl})\end{aligned}
699 | where :math:`k` is a force constant in units of energy and :math:`f`
700 is a cubic spline function; for details see :ref:`cubicspline`. For
701 each interaction, the force constant :math:`k` and the table number
702 :math:`n` are specified in the topology. There are two different types
703 of bonds, one that generates exclusions (type 8) and one that does not
704 (type 9). For details see :numref:`Table %s <tab-topfile2>`. The table files are
705 supplied to the :ref:`mdrun <gmx mdrun>` program. After the table file name an
706 underscore, the letter “b” for bonds, “a” for angles or “d” for
707 dihedrals and the table number must be appended. For example, a
708 tabulated bond with :math:`n=0` can be read from the file
709 table\_b0.xvg. Multiple tables can be supplied simply by adding files
710 with different values of :math:`n`, and are applied to the appropriate
711 bonds, as specified in the topology (:numref:`Table %s <tab-topfile2>`). The format
712 for the table files is three fixed-format columns of any suitable
713 width. These columns must contain :math:`x`, :math:`f(x)`,
714 :math:`-f'(x)`, and the values of :math:`x` should be uniformly
715 spaced. Requirements for entries in the topology are given
716 in :numref:`Table %s <tab-topfile2>`. The setup of the tables is as follows:
717 | **bonds**: :math:`x` is the distance in nm. For distances beyond the
718 table length, :ref:`mdrun <gmx mdrun>` will quit with an error message.
719 | **angles**: :math:`x` is the angle in degrees. The table should go
720 from 0 up to and including 180 degrees; the derivative is taken in
722 | **dihedrals**: :math:`x` is the dihedral angle in degrees. The table
723 should go from -180 up to and including 180 degrees; the IUPAC/IUB
724 convention is used, *i.e.* zero is cis, the derivative is taken in