Bonded interactions

Bonded interactions are based on a fixed list of atoms. They are not exclusively pair interactions, but include 3- and 4-body interactions as well. There are bond stretching (2-body), bond angle (3-body), and dihedral angle (4-body) interactions. A special type of dihedral interaction (called improper dihedral) is used to force atoms to remain in a plane or to prevent transition to a configuration of opposite chirality (a mirror image).

Bond stretching

Harmonic potential

The bond stretching between two covalently bonded atoms \(i\) and \(j\) is represented by a harmonic potential:

Fig. 19 Principle of bond stretching (left), and the bond stretching potential (right).

(1)\[V_b~({r_{ij}}) = {\frac{1}{2}}k^b_{ij}({r_{ij}}-b_{ij})^2\]

See also Fig. 19, with the force given by:

(2)\[\mathbf{F}_i(\mathbf{r}_{ij}) = k^b_{ij}({r_{ij}}-b_{ij}) {\frac{{\mathbf{r}_{ij}}}{{r_{ij}}}}\]

Fourth power potential

In the GROMOS-96 force field 77, the covalent bond potential is, for reasons of computational efficiency, written as:

(3)\[V_b~({r_{ij}}) = \frac{1}{4}k^b_{ij}\left({r_{ij}}^2-b_{ij}^2\right)^2\]

The corresponding force is:

(4)\[\mathbf{F}_i(\mathbf{r}_{ij}) = k^b_{ij}({r_{ij}}^2-b_{ij}^2)~\mathbf{r}_ij\]

The force constants for this form of the potential are related to the usual harmonic force constant \(k^{b,\mathrm{harm}}\) (sec. Bond stretching) as

(5)\[2 k^b b_{ij}^2 = k^{b,\mathrm{harm}}\]

The force constants are mostly derived from the harmonic ones used in GROMOS-87 78. Although this form is computationally more efficient (because no square root has to be evaluated), it is conceptually more complex. One particular disadvantage is that since the form is not harmonic, the average energy of a single bond is not equal to \({\frac{1}{2}}kT\) as it is for the normal harmonic potential.

Morse potential bond stretching

For some systems that require an anharmonic bond stretching potential, the Morse potential 79 between two atoms i and j is available in GROMACS. This potential differs from the harmonic potential in that it has an asymmetric potential well and a zero force at infinite distance. The functional form is:

(6)\[\displaystyle V_{morse} (r_{ij}) = D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2,\]

See also Fig. 20, and the corresponding force is:

(7)\[\begin{split}\begin{array}{rcl} \displaystyle {\bf F}_{morse} ({\bf r}_{ij})&=&2 D_{ij} \beta_{ij} \exp(-\beta_{ij}(r_{ij}-b_{ij})) * \\ \displaystyle \: & \: &[1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))] \frac{\displaystyle {\bf r}_{ij}}{\displaystyle r_{ij}}, \end{array}\end{split}\]

where \(\displaystyle D_{ij}\) is the depth of the well in kJ/mol, \(\displaystyle \beta_{ij}\) defines the steepness of the well (in nm\(^{-1}\)), and \(\displaystyle b_{ij}\) is the equilibrium distance in nm. The steepness parameter \(\displaystyle \beta_{ij}\) can be expressed in terms of the reduced mass of the atoms i and j, the fundamental vibration frequency \(\displaystyle\omega_{ij}\) and the well depth \(\displaystyle D_{ij}\):

(8)\[\displaystyle \beta_{ij}= \omega_{ij} \sqrt{\frac{\mu_{ij}}{2 D_{ij}}}\]

and because \(\displaystyle \omega = \sqrt{k/\mu}\), one can rewrite \(\displaystyle \beta_{ij}\) in terms of the harmonic force constant \(\displaystyle k_{ij}\):

(9)\[\displaystyle \beta_{ij}= \sqrt{\frac{k_{ij}}{2 D_{ij}}}\]

For small deviations \(\displaystyle (r_{ij}-b_{ij})\), one can approximate the \(\displaystyle \exp\)-term to first-order using a Taylor expansion:

(10)\[\displaystyle \exp(-x) \approx 1-x\]

and substituting (9) and (10) in the functional form:

(11)\[\begin{split}\begin{array}{rcl} \displaystyle V_{morse} (r_{ij})&=&D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2\\ \displaystyle \:&=&D_{ij} [1 - (1 -\sqrt{\frac{k_{ij}}{2 D_{ij}}}(r_{ij}-b_{ij}))]^2\\ \displaystyle \:&=&\frac{1}{2} k_{ij} (r_{ij}-b_{ij}))^2 \end{array}\end{split}\]

we recover the harmonic bond stretching potential.

Fig. 20 The Morse potential well, with bond length 0.15 nm.

Cubic bond stretching potential

Another anharmonic bond stretching potential that is slightly simpler than the Morse potential adds a cubic term in the distance to the simple harmonic form:

(12)\[V_b~({r_{ij}}) = k^b_{ij}({r_{ij}}-b_{ij})^2 + k^b_{ij}k^{cub}_{ij}({r_{ij}}-b_{ij})^3\]

A flexible water model (based on the SPC water model 80) including a cubic bond stretching potential for the O-H bond was developed by Ferguson 81. This model was found to yield a reasonable infrared spectrum. The Ferguson water model is available in the GROMACS library (flexwat-ferguson.itp). It should be noted that the potential is asymmetric: overstretching leads to infinitely low energies. The integration timestep is therefore limited to 1 fs.

The force corresponding to this potential is:

(13)\[\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}}}}\]

FENE bond stretching potential

In coarse-grained polymer simulations the beads are often connected by a FENE (finitely extensible nonlinear elastic) potential 82:

(14)\[V_{\mbox{FENE}}({r_{ij}}) = -{\frac{1}{2}}k^b_{ij} b^2_{ij} \log\left(1 - \frac{{r_{ij}}^2}{b^2_{ij}}\right)\]

The potential looks complicated, but the expression for the force is simpler:

(15)\[F_{\mbox{FENE}}(\mathbf{r}_{ij}) = -k^b_{ij} \left(1 - \frac{{r_{ij}}^2}{b^2_{ij}}\right)^{-1} \mathbf{r}_{ij}\]

At short distances the potential asymptotically goes to a harmonic potential with force constant \(k^b\), while it diverges at distance \(b\).

Harmonic angle potential

The bond-angle vibration between a triplet of atoms \(i\) - \(j\) - \(k\) is also represented by a harmonic potential on the angle \({\theta_{ijk}}\)

Fig. 21 Principle of angle vibration (left) and the bond angle potential.

(16)\[V_a({\theta_{ijk}}) = {\frac{1}{2}}k^{\theta}_{ijk}({\theta_{ijk}}-{\theta_{ijk}}^0)^2\]

As the bond-angle vibration is represented by a harmonic potential, the form is the same as the bond stretching (Fig. 19).

The force equations are given by the chain rule:

(17)\[\begin{split}\begin{array}{l} \mathbf{F}_i ~=~ -\displaystyle\frac{d V_a({\theta_{ijk}})}{d \mathbf{r}_i} \\ \mathbf{F}_k ~=~ -\displaystyle\frac{d V_a({\theta_{ijk}})}{d \mathbf{r}_k} \\ \mathbf{F}_j ~=~ -\mathbf{F}_i-\mathbf{F}_k \end{array} ~ \mbox{ ~ where ~ } ~ {\theta_{ijk}}= \arccos \frac{(\mathbf{r}_{ij} \cdot \mathbf{r}_{kj})}{r_{ij}r_{kj}}\end{split}\]

The numbering \(i,j,k\) is in sequence of covalently bonded atoms. Atom \(j\) is in the middle; atoms \(i\) and \(k\) are at the ends (see Fig. 21). Note that in the input in topology files, angles are given in degrees and force constants in kJ/mol/rad\(^2\).

Cosine based angle potential

In the GROMOS-96 force field a simplified function is used to represent angle vibrations:

(18)\[V_a({\theta_{ijk}}) = {\frac{1}{2}}k^{\theta}_{ijk}\left(\cos({\theta_{ijk}}) - \cos({\theta_{ijk}}^0)\right)^2\]

where

(19)\[\cos({\theta_{ijk}}) = \frac{\mathbf{r}_{ij}\cdot\mathbf{r}_{kj}}{{r_{ij}}r_{kj}}\]

The corresponding force can be derived by partial differentiation with respect to the atomic positions. The force constants in this function are related to the force constants in the harmonic form \(k^{\theta,\mathrm{harm}}\) (Harmonic angle potential) by:

(20)\[k^{\theta} \sin^2({\theta_{ijk}}^0) = k^{\theta,\mathrm{harm}}\]

In the GROMOS-96 manual there is a much more complicated conversion formula which is temperature dependent. The formulas are equivalent at 0 K and the differences at 300 K are on the order of 0.1 to 0.2%. Note that in the input in topology files, angles are given in degrees and force constants in kJ/mol.

Restricted bending potential

The restricted bending (ReB) potential 83 prevents the bending angle \(\theta\) from reaching the \(180^{\circ}\) value. In this way, the numerical instabilities due to the calculation of the torsion angle and potential are eliminated when performing coarse-grained molecular dynamics simulations.

To systematically hinder the bending angles from reaching the \(180^{\circ}\) value, the bending potential (18) is divided by a \(\sin^2\theta\) factor:

(21)\[V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} \frac{(\cos\theta_i - \cos\theta_0)^2}{\sin^2\theta_i}.\]

Figure 22 shows the comparison between the ReB potential, (21), and the standard one (18).

Fig. 22 Bending angle potentials: cosine harmonic (solid black line), angle harmonic (dashed black line) and restricted bending (red) with the same bending constant \(k_{\theta}=85\) kJ mol\(^{-1}\) and equilibrium angle \(\theta_0=130^{\circ}\). The orange line represents the sum of a cosine harmonic (\(k =50\) kJ mol\(^{-1}\)) with a restricted bending (\(k =25\) kJ mol\(^{-1}\)) potential, both with \(\theta_0=130^{\circ}\).

The wall of the ReB potential is very repulsive in the region close to \(180^{\circ}\) and, as a result, the bending angles are kept within a safe interval, far from instabilities. The power \(2\) of \(\sin\theta_i\) in the denominator has been chosen to guarantee this behavior and allows an elegant differentiation:

(22)\[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}}.\]

Due to its construction, the restricted bending potential cannot be used for equilibrium \(\theta_0\) values too close to \(0^{\circ}\) or \(180^{\circ}\) (from experience, at least \(10^{\circ}\) difference is recommended). It is very important that, in the starting configuration, all the bending angles have to be in the safe interval to avoid initial instabilities. This bending potential can be used in combination with any form of torsion potential. It will always prevent three consecutive particles from becoming collinear and, as a result, any torsion potential will remain free of singularities. It can be also added to a standard bending potential to affect the angle around \(180^{\circ}\), but to keep its original form around the minimum (see the orange curve in Fig. 22).

Urey-Bradley potential

The Urey-Bradley bond-angle vibration between a triplet of atoms \(i\) - \(j\) - \(k\) is represented by a harmonic potential on the angle \({\theta_{ijk}}\) and a harmonic correction term on the distance between the atoms \(i\) and \(k\). Although this can be easily written as a simple sum of two terms, it is convenient to have it as a single entry in the topology file and in the output as a separate energy term. It is used mainly in the CHARMm force field 84. The energy is given by:

(23)\[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\]

The force equations can be deduced from sections Harmonic potential and Harmonic angle potential.

Bond-Bond cross term

The bond-bond cross term for three particles \(i, j, k\) forming bonds \(i-j\) and \(k-j\) is given by 85:

(24)\[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)\]

where \(k_{rr'}\) is the force constant, and \(r_{1e}\) and \(r_{2e}\) are the equilibrium bond lengths of the \(i-j\) and \(k-j\) bonds respectively. The force associated with this potential on particle \(i\) is:

(25)\[\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|}\]

The force on atom \(k\) can be obtained by swapping \(i\) and \(k\) in the above equation. Finally, the force on atom \(j\) follows from the fact that the sum of internal forces should be zero: \(\mathbf{F}_j = -\mathbf{F}_i-\mathbf{F}_k\).

Bond-Angle cross term

The bond-angle cross term for three particles \(i, j, k\) forming bonds \(i-j\) and \(k-j\) is given by 85:

(26)\[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)\]

where \(k_{r\theta}\) is the force constant, \(r_{3e}\) is the \(i-k\) distance, and the other constants are the same as in (24). The force associated with the potential on atom \(i\) is:

(27)\[\mathbf{F}_{i} ~=~ -k_{r\theta} \left[ \left( \left| \mathbf{r}_{i} - \mathbf{r}_{k}\right| -r_{3e}\right) \frac{ \mathbf{r}_{i}-\mathbf{r}_j} { \left| \mathbf{r}_{i}-\mathbf{r}_{j}\right| } + \left( \left| \mathbf{r}_{i}-\mathbf{r}_{j}\right| -r_{1e} + \left| \mathbf{r}_{k}-\mathbf{r}_{j}\right| -r_{2e}\right) \frac{ \mathbf{r}_{i}-\mathbf{r}_{k}} {\left| \mathbf{r}_{i}-\mathbf{r}_{k}\right| } \right]\]

Quartic angle potential

For special purposes there is an angle potential that uses a fourth order polynomial:

(28)\[V_q({\theta_{ijk}}) ~=~ \sum_{n=0}^5 C_n ({\theta_{ijk}}-{\theta_{ijk}}^0)^n\]

Improper dihedrals

Improper dihedrals are meant to keep planar groups (e.g. aromatic rings) planar, or to prevent molecules from flipping over to their mirror images, see Fig. 23.

Fig. 23 Principle of improper dihedral angles. Out of plane bending for rings. The improper dihedral angle \(\xi\) is defined as the angle between planes (i,j,k) and (j,k,l).

Fig. 24 Principle of improper dihedral angles. Out of tetrahedral angle. The improper dihedral angle \(\xi\) is defined as the angle between planes (i,j,k) and (j,k,l).

Improper dihedrals: harmonic type

The simplest improper dihedral potential is a harmonic potential; it is plotted in Fig. 25.

(29)\[V_{id}(\xi_{ijkl}) = {\frac{1}{2}}k_{\xi}(\xi_{ijkl}-\xi_0)^2\]

Since the potential is harmonic it is discontinuous, but since the discontinuity is chosen at 180\(^\circ\) distance from \(\xi_0\) this will never cause problems. Note that in the input in topology files, angles are given in degrees and force constants in kJ/mol/rad\(^2\).

Fig. 25 Improper dihedral potential.

Improper dihedrals: periodic type

This potential is identical to the periodic proper dihedral (see below). There is a separate dihedral type for this (type 4) only to be able to distinguish improper from proper dihedrals in the parameter section and the output.

Proper dihedrals

For the normal dihedral interaction there is a choice of either the GROMOS periodic function or a function based on expansion in powers of \(\cos \phi\) (the so-called Ryckaert-Bellemans potential). This choice has consequences for the inclusion of special interactions between the first and the fourth atom of the dihedral quadruple. With the periodic GROMOS potential a special 1-4 LJ-interaction must be included; with the Ryckaert-Bellemans potential for alkanes the 1-4 interactions must be excluded from the non-bonded list. Note: Ryckaert-Bellemans potentials are also used in e.g. the OPLS force field in combination with 1-4 interactions. You should therefore not modify topologies generated by pdb2gmx in this case.

Proper dihedrals: periodic type

Proper dihedral angles are defined according to the IUPAC/IUB convention, where \(\phi\) is the angle between the \(ijk\) and the \(jkl\) planes, with zero corresponding to the cis configuration (\(i\) and \(l\) on the same side). There are two dihedral function types in GROMACS topology files. There is the standard type 1 which behaves like any other bonded interactions. For certain force fields, type 9 is useful. Type 9 allows multiple potential functions to be applied automatically to a single dihedral in the [ dihedral ] section when multiple parameters are defined for the same atomtypes in the [ dihedraltypes ] section.

Fig. 26 Principle of proper dihedral angle (left, in trans form) and the dihedral angle potential (right).

(30)\[V_d(\phi_{ijkl}) = k_{\phi}(1 + \cos(n \phi - \phi_s))\]

Proper dihedrals: Ryckaert-Bellemans function

For alkanes, the following proper dihedral potential is often used (see Fig. 27):

(31)\[V_{rb}(\phi_{ijkl}) = \sum_{n=0}^5 C_n( \cos(\psi ))^n,\]

where \(\psi = \phi - 180^\circ\).

Note: A conversion from one convention to another can be achieved by multiplying every coefficient \(\displaystyle C_n\) by \(\displaystyle (-1)^n\).

An example of constants for \(C\) is given in Table 8.

Table 8 Constants for Ryckaert-Bellemans potential (\(\mathrm{kJ mol}^{-1}\)).

\(C_0\)	9.28	\(C_2\)	-13.12	\(C_4\)	26.24
\(C_1\)	12.16	\(C_3\)	-3.06	\(C_5\)	-31.5

Fig. 27 Ryckaert-Bellemans dihedral potential.

(Note: The use of this potential implies exclusion of LJ interactions between the first and the last atom of the dihedral, and \(\psi\) is defined according to the “polymer convention” (\(\psi_{trans}=0\)).)

The RB dihedral function can also be used to include Fourier dihedrals (see below):

(32)\[V_{rb} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[F_1(1+\cos(\phi)) + F_2( 1-\cos(2\phi)) + F_3(1+\cos(3\phi)) + F_4(1-\cos(4\phi))\right]\]

Because of the equalities \(\cos(2\phi) = 2\cos^2(\phi) - 1\), \(\cos(3\phi) = 4\cos^3(\phi) - 3\cos(\phi)\) and \(\cos(4\phi) = 8\cos^4(\phi) - 8\cos^2(\phi) + 1\) one can translate the OPLS parameters to Ryckaert-Bellemans parameters as follows:

(33)\[\begin{split}\displaystyle \begin{array}{rcl} \displaystyle C_0&=&F_2 + \frac{1}{2} (F_1 + F_3)\\ \displaystyle C_1&=&\frac{1}{2} (- F_1 + 3 \, F_3)\\ \displaystyle C_2&=& -F_2 + 4 \, F_4\\ \displaystyle C_3&=&-2 \, F_3\\ \displaystyle C_4&=&-4 \, F_4\\ \displaystyle C_5&=&0 \end{array}\end{split}\]

with OPLS parameters in protein convention and RB parameters in polymer convention (this yields a minus sign for the odd powers of cos\((\phi)\)).

Note: Mind the conversion from kcal mol\(^{-1}\) for literature OPLS and RB parameters to kJ mol\(^{-1}\) in GROMACS.

Proper dihedrals: Fourier function

The OPLS potential function is given as the first three  86 or four 87 cosine terms of a Fourier series. In GROMACS the four term function is implemented:

(34)\[V_{F} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[C_1(1+\cos(\phi)) + C_2( 1-\cos(2\phi)) + C_3(1+\cos(3\phi)) + C_4(1-\cos(4\phi))\right],\]

Internally, GROMACS uses the Ryckaert-Bellemans code to compute Fourier dihedrals (see above), because this is more efficient.

Note: Mind the conversion from kcal mol\(^{-1}\) for literature OPLS parameters to kJ mol\(^{-1}\) in GROMACS.

Proper dihedrals: Restricted torsion potential

In a manner very similar to the restricted bending potential (see Restricted bending potential), a restricted torsion/dihedral potential is introduced:

(35)\[V_{\rm ReT}(\phi_i) = \frac{1}{2} k_{\phi} \frac{(\cos\phi_i - \cos\phi_0)^2}{\sin^2\phi_i}\]

with the advantages of being a function of \(\cos\phi\) (no problems taking the derivative of \(\sin\phi\)) and of keeping the torsion angle at only one minimum value. In this case, the factor \(\sin^2\phi\) does not allow the dihedral angle to move from the [\(-180^{\circ}\):0] to [0:\(180^{\circ}\)] interval, i.e. it cannot have maxima both at \(-\phi_0\) and \(+\phi_0\) maxima, but only one of them. For this reason, all the dihedral angles of the starting configuration should have their values in the desired angles interval and the the equilibrium \(\phi_0\) value should not be too close to the interval limits (as for the restricted bending potential, described in Restricted bending potential, at least \(10^{\circ}\) difference is recommended).

Proper dihedrals: Combined bending-torsion potential

When the four particles forming the dihedral angle become collinear (this situation will never happen in atomistic simulations, but it can occur in coarse-grained simulations) the calculation of the torsion angle and potential leads to numerical instabilities. One way to avoid this is to use the restricted bending potential (see Restricted bending potential) that prevents the dihedral from reaching the \(180^{\circ}\) value.

Another way is to disregard any effects of the dihedral becoming ill-defined, keeping the dihedral force and potential calculation continuous in entire angle range by coupling the torsion potential (in a cosine form) with the bending potentials of the adjacent bending angles in a unique expression:

(36)\[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}.\]

This combined bending-torsion (CBT) potential has been proposed by 88 for polymer melt simulations and is extensively described in 83.

This potential has two main advantages:

• it does not only depend on the dihedral angle \(\phi_i\) (between the \(i-2\), \(i-1\), \(i\) and \(i+1\) beads) but also on the bending angles \(\theta_{i-1}\) and \(\theta_i\) defined from three adjacent beads (\(i-2\), \(i-1\) and \(i\), and \(i-1\), \(i\) and \(i+1\), respectively). The two \(\sin^3\theta\) pre-factors, tentatively suggested by 89 and theoretically discussed by 90, cancel the torsion potential and force when either of the two bending angles approaches the value of \(180^\circ\).
• its dependence on \(\phi_i\) is expressed through a polynomial in \(\cos\phi_i\) that avoids the singularities in \(\phi=0^\circ\) or \(180^\circ\) in calculating the torsional force.

These two properties make the CBT potential well-behaved for MD simulations with weak constraints on the bending angles or even for steered / non-equilibrium MD in which the bending and torsion angles suffer major modifications. When using the CBT potential, the bending potentials for the adjacent \(\theta_{i-1}\) and \(\theta_i\) may have any form. It is also possible to leave out the two angle bending terms (\(\theta_{i-1}\) and \(\theta_{i}\)) completely. Fig. 28 illustrates the difference between a torsion potential with and without the \(\sin^{3}\theta\) factors (blue and gray curves, respectively).

Fig. 28 Blue: surface plot of the combined bending-torsion potential ((36) with \(k = 10\) kJ mol\(^{-1}\), \(a_0=2.41\), \(a_1=-2.95\), \(a_2=0.36\), \(a_3=1.33\)) when, for simplicity, the bending angles behave the same (\(\theta_1=\theta_2=\theta\)). Gray: the same torsion potential without the \(\sin^{3}\theta\) terms (Ryckaert-Bellemans type). \(\phi\) is the dihedral angle.

Additionally, the derivative of \(V_{CBT}\) with respect to the Cartesian variables is straightforward:

(37)\[\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}} + \frac{\partial V_{\rm CBT}}{\partial \theta_{i }} \frac{\partial \theta_{i }}{\partial \vec r_{l}} + \frac{\partial V_{\rm CBT}}{\partial \phi_{i }} \frac{\partial \phi_{i }}{\partial \vec r_{l}}\]

The CBT is based on a cosine form without multiplicity, so it can only be symmetrical around \(0^{\circ}\). To obtain an asymmetrical dihedral angle distribution (e.g. only one maximum in [\(-180^{\circ}\):\(180^{\circ}\)] interval), a standard torsion potential such as harmonic angle or periodic cosine potentials should be used instead of a CBT potential. However, these two forms have the inconveniences of the force derivation (\(1/\sin\phi\)) and of the alignment of beads (\(\theta_i\) or \(\theta_{i-1} = 0^{\circ}, 180^{\circ}\)). Coupling such non-\(\cos\phi\) potentials with \(\sin^3\theta\) factors does not improve simulation stability since there are cases in which \(\theta\) and \(\phi\) are simultaneously \(180^{\circ}\). The integration at this step would be possible (due to the cancelling of the torsion potential) but the next step would be singular (\(\theta\) is not \(180^{\circ}\) and \(\phi\) is very close to \(180^{\circ}\)).

Tabulated bonded interaction functions

For full flexibility, any functional shape can be used for bonds, angles and dihedrals through user-supplied tabulated functions. The functional shapes are:

(38)\[\begin{split}\begin{aligned} V_b(r_{ij}) &=& k \, f^b_n(r_{ij}) \\ V_a({\theta_{ijk}}) &=& k \, f^a_n({\theta_{ijk}}) \\ V_d(\phi_{ijkl}) &=& k \, f^d_n(\phi_{ijkl})\end{aligned}\end{split}\]

where \(k\) is a force constant in units of energy and \(f\) is a cubic spline function; for details see Cubic splines for potentials. For each interaction, the force constant \(k\) and the table number \(n\) are specified in the topology. There are two different types of bonds, one that generates exclusions (type 8) and one that does not (type 9). For details see Table 14. The table files are supplied to the mdrun program. After the table file name an underscore, the letter “b” for bonds, “a” for angles or “d” for dihedrals and the table number must be appended. For example, a tabulated bond with \(n=0\) can be read from the file table_b0.xvg. Multiple tables can be supplied simply by adding files with different values of \(n\), and are applied to the appropriate bonds, as specified in the topology (Table 14). The format for the table files is three fixed-format columns of any suitable width. These columns must contain \(x\), \(f(x)\), \(-f'(x)\), and the values of \(x\) should be uniformly spaced. Requirements for entries in the topology are given in Table 14. The setup of the tables is as follows:

bonds: \(x\) is the distance in nm. For distances beyond the table length, mdrun will quit with an error message.

angles: \(x\) is the angle in degrees. The table should go from 0 up to and including 180 degrees; the derivative is taken in degrees.

dihedrals: \(x\) is the dihedral angle in degrees. The table should go from -180 up to and including 180 degrees; the IUPAC/IUB convention is used, i.e. zero is cis, the derivative is taken in degrees.