Update instructions in containers.rst
[gromacs.git] / docs / reference-manual / functions / bonded-interactions.rst
blobe8b9cffd1eb76136b0bed5ddc97dc51683f981fd
1 Bonded interactions
2 -------------------
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).
12 .. _bondpot:
14 Bond stretching
15 ~~~~~~~~~~~~~~~
17 .. _harmonicbond:
19 Harmonic potential
20 ^^^^^^^^^^^^^^^^^^
22 The bond stretching between two covalently bonded atoms :math:`i` and
23 :math:`j` is represented by a harmonic potential:
25 .. _fig-bstretch1:
27 .. figure:: plots/bstretch.*
28    :width: 7.00000cm
30    Principle of bond stretching (left), and the bond stretching
31    potential (right).
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
41 .. _g96bond:
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
50           :label: eqng96bond
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
70 harmonic potential.
72 .. _morsebond:
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,
84           :label: eqnmorsebond
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}},
91           \end{array}
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}}}
102           :label: eqnmorsefreq
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}}}
109           :label: eqnbetaij
111 For small deviations :math:`\displaystyle (r_{ij}-b_{ij})`, one can
112 approximate the :math:`\displaystyle \exp`-term to first-order using a
113 Taylor expansion:
115 .. math:: \displaystyle \exp(-x) \approx 1-x
116           :label: eqnexpminx
118 and substituting :eq:`eqn. %s <eqnbetaij>` and :eq:`eqn. %s <eqnexpminx>` in the
119 functional form:
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
125           \end{array}
126           :label: eqnharmfrommorse
128 we recover the harmonic bond stretching potential.
130 .. _fig-morse:
132 .. figure:: plots/f-morse.*
133    :width: 7.00000cm
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
142 harmonic form:
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
145           :label: eqncubicbond
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)
168           :label: eqnfenebond
170 The potential looks complicated, but the expression for the force is
171 simpler:
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
179 :math:`b`.
181 .. _harmonicangle:
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}}`
190 .. _fig-angle:
192 .. figure:: plots/angle.*
193    :width: 7.00000cm
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
198           :label: eqnharmangle
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
210           \end{array}
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`.
221 .. _g96angle:
223 Cosine based angle potential
224 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
226 In the GROMOS-96 force field a simplified function is used to represent
227 angle vibrations:
229 .. math:: V_a({\theta_{ijk}}) = {\frac{1}{2}}k^{\theta}_{ijk}\left(\cos({\theta_{ijk}}) - \cos({\theta_{ijk}}^0)\right)^2
230           :label: eqnG96angle
232 where
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.
251 .. _reb:
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}.
267           :label: eqnReB
269 :numref:`Figure %s <fig-ReB>` shows the comparison between the ReB potential,
270 :eq:`%s <eqnReB>`, and the standard one :eq:`%s <eqnG96angle>`.
272 .. _fig-ReB:
274 .. figure:: plots/fig-02.*
275    :width: 10.00000cm
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}}.
293           :label: eqdiffReB
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
321           :label: eqnUBAngle
323 The force equations can be deduced from sections :ref:`harmonicbond`
324 and :ref:`harmonicangle`.
326 Bond-Bond cross term
327 ~~~~~~~~~~~~~~~~~~~~
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)
334           :label: eqncrossbb
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)
357           :label: eqncrossba
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
362 :math:`i` is:
364 .. math:: \mathbf{F}_{i} ~=~ -k_{r\theta}
365           \left[
366           \left(
367           \left| \mathbf{r}_{i} - \mathbf{r}_{k}\right|
368           -r_{3e}\right)
369           \frac{
370                 \mathbf{r}_{i}-\mathbf{r}_j}
371                 { \left| \mathbf{r}_{i}-\mathbf{r}_{j}\right| 
372                 }
373           + \left(
374             \left| \mathbf{r}_{i}-\mathbf{r}_{j}\right|
375           -r_{1e}
376           + \left| \mathbf{r}_{k}-\mathbf{r}_{j}\right|
377           -r_{2e}\right)
378           \frac{
379                 \mathbf{r}_{i}-\mathbf{r}_{k}}
380                 {\left| \mathbf{r}_{i}-\mathbf{r}_{k}\right|
381                 }
382           \right]
383           :label: eqncrossbaforce
385 Quartic angle potential
386 ~~~~~~~~~~~~~~~~~~~~~~~
388 For special purposes there is an angle potential that uses a fourth
389 order polynomial:
391 .. math:: V_q({\theta_{ijk}}) ~=~ \sum_{n=0}^5 C_n ({\theta_{ijk}}-{\theta_{ijk}}^0)^n
392           :label: eqnquarticangle
394 .. _imp:
396 Improper dihedrals
397 ~~~~~~~~~~~~~~~~~~
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>`.
403 .. _fig-imp:
405 .. figure:: plots/ring-imp.*
406         :width: 4.00000cm
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.*
413         :width: 3.00000cm
415 .. figure:: plots/tetra-im.*
416         :width: 3.00000cm
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`.
437 .. _fig-imps:
439 .. figure:: plots/f-imps.*
440    :width: 10.00000cm
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
450 the output.
452 Proper dihedrals
453 ~~~~~~~~~~~~~~~~
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 ]``
480 section.
482 .. _fig-pdihf:
484 .. figure:: plots/f-dih.*
485    :width: 7.00000cm
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>`.
509 .. _tab-crb:
511 .. table:: 
512     Constants for Ryckaert-Bellemans potential (\ :math:`\mathrm{kJ mol}^{-1}`).
513     :widths: auto
514     :align: center
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     +-------------+-------+-------------+--------+-------------+-------+
523 .. _fig-rbdih:
525 .. figure:: plots/f-rbs.*
526    :width: 8.00000cm
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
536   (see below):
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
546   follows:
548   .. math:: \displaystyle
549             \begin{array}{rcl}
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
556             \end{array}
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
564   |Gromacs|.
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
572   implemented:
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}
590           :label: eqnReT
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
603 recommended).
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}.
622           :label: eqnCBT
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
643    torsional force.
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).
656 .. _fig-CBT:
658 .. figure:: plots/fig-04.*
659    :width: 10.00000cm
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}}
675           :label: eqnforcecbt
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
729   degrees.
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
733   degrees.