The Intramolecular Potential Functions

In this section we catalogue and describe the forms of potential function available in . The keywords required to select potential forms are given in brackets () before each definition. The derivations of the atomic forces, virial and stress tensor are also outlined.

Bond Potentials

Fig. 1 The interatomic bond vector

The bond potentials describe explicit chemical bonds between specified atoms. They are all functions of the interatomic distance. Only the coulomb potential makes an exception as it depends on the charges of the specified atoms. The potential functions available are as follows:

1. Harmonic bond: (harm)

   (2)\[U(r_{ij}) = \frac{1}{2} k (r_{ij}-r_{o})^{2}\]

2. Morse potential: (mors)

   (3)\[U(r_{ij}) = E_{o} [\{1-\exp(-k(r_{ij}-r_{o}))\}^{2}-1]\]

3. 12-6 potential bond: (12-6)

   (4)\[U(r_{ij}) = \left(\frac{A}{r_{ij}^{12}}\right)-\left(\frac{B}{r_{ij}^{6}}\right)\]

4. Lennard-Jones potential: (lj)

   (5)\[U(r_{ij}) = 4\epsilon\left[\left (\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right]\]

5. Restrained harmonic: (rhrm)

   (6)\[\begin{split}U(r_{ij}) = \left\{ \begin{array} {l@{\quad:\quad}l} \frac{1}{2}k(r_{ij}-r_{o})^{2} & |r_{ij}-r_{o}|\le r_{c} \\ \frac{1}{2}kr_{c}^{2}+kr_{c}(|r_{ij}-r_{o}|-r_{c}) & |r_{ij}-r_{o}| > r_{c} \end{array} \right.\end{split}\]

6. Quartic potential: (quar)

   (7)\[U(r_{ij}) = \frac{k}{2}(r_{ij}-r_{o})^{2}+\frac{k'}{3}(r_{ij}-r_{o})^{3}+\frac{k''}{4}(r_{ij}-r_{o})^{4}\]

7. Buckingham potential: (buck)

   (8)\[U(r_{ij}) = A~\exp\left(-\frac{r_{ij}}{\rho}\right)-\frac{C}{r_{ij}^{6}}\]

8. Coulomb potential: (coul)

   (9)\[U(r_{ij}) = k \cdot U^{\text{Electrostatics}}(r_{ij}) \; \left(= \frac{k}{4\pi\epsilon_{0}\epsilon}\frac{q_{i}q_{j}}{r_{ij}}\right)~~,\]

   where \(q_{\ell}\) is the charge on an atom labelled \(\ell\). It is worth noting that the Coulomb potential switches to the particular model of Electrostatics opted in CONTROL.

9. Shifted finitely extendible non-linear elastic (FENE) potential ^10,38,122: (fene)

   (10)\[\begin{split}U(r_{ij}) = \left\{ \begin{array} {l@{\quad:\quad}l} -0.5~k~R_{o}^{2}~ln\left[1-\left(\frac{r_{ij}-\Delta}{R_{o}}\right)^{2}\right] & |r_{ij} - \Delta| < R_{o} \\ \infty & |r_{ij} - \Delta| \ge R_{o} \end{array} \right. \label{FENE}\end{split}\]

   The FENE potential is used to maintain the distance between connected beads and to prevent chains from crossing each other. It is used in combination with the WCA, equation (60), potential to create a potential well for the flexible bonds of a molecule, that maintains the topology of the molecule. This implementation allows for a radius shift of up to half a \(R_{o}\) (\(|\Delta| \le0.5~R_{o}\)) with a default of zero (\(\Delta_{default} = 0\)).

10. MM3 bond stretch potential ⁵: (mmst)

    (11)\[U(r_{ij}) = k~(r_{ij}-r_{o})^{2}\left[1-2.55~(r_{ij}-r_{o})+(7/12)~2.55^{2}~(r_{ij}-r_{o})^{2}\right]\]

11. Tabulated potential: (tab). The potential is defined numerically in TABBND (see Section Setting up Tabulated Intramolecular Force-Field Files and Section  The TABBND, TABANG, TABDIH & TABINV Files).

In these formulae \(r_{ij}\) is the distance between atoms labelled \(i\) and \(j\):

\[r_{ij} = |\underline{r}_{j}-\underline{r}_{i}|,\]

where \(\underline{r}_{\ell}\) is the position vector of an atom labelled \(\ell\).

Note

some DL_POLY_4 routines may use the convention that \(\underline{r_{ij}}=\underline{r}_{i}-\underline{r}_{j}\)

The force on the atom \(j\) arising from a bond potential is obtained using the general formula:

(12)\[\underline{f}_{j} = -\frac{1}{{r}_{ij}} \left[ \frac{\partial }{\partial r_{ij}}U(r_{ij})\right] \underline{r}_{ij}~~. \label{bondf}\]

The force \(\underline{f}_{i}\) acting on atom \(i\) is the negative of this.

The contribution to be added to the atomic virial is given by

\[{\cal W} = -\underline{r}_{ij} \cdot \underline{f}_{j}~~,\]

with only one such contribution from each bond.

The contribution to be added to the atomic stress tensor is given by

(13)\[\sigma^{\alpha \beta} = r_{ij}^{\alpha} f_{j}^{\beta}~~, \label{bonds}\]

where \(\alpha\) and \(\beta\) indicate the \(x,y,z\) components. The atomic stress tensor derived in this way is symmetric.

In DL_POLY_4 bond forces are handled by the routine bonds_forces (and intra_coul called within).

Distance Restraints

In DL_POLY_4 distance restraints, in which the separation between two atoms, is maintained around some preset value \(r_0\) is handled as a special case of bond potentials. As a consequence, distance restraints may be applied only between atoms in the same molecule. Unlike with application of the “pure” bond potentials, the electrostatic and van der Waals interactions between the pair of atoms are still evaluated when distance restraints are applied. All the potential forms of the previous section are available as distance restraints, although they have different key words:

1. Harmonic potential: (-hrm)

2. Morse potential: (-mrs)

3. 12-6 potential bond: (-126)

4. Lennard-Jones potential: (-lj)

5. Restrained harmonic: (-rhm)

6. Quartic potential: (-qur)

7. Buckingham potential: (-bck)

8. Coulomb potential: (-cul)

9. FENE potential: (-fne)

10. MM3 bond stretch potential ⁵: (-m3s)

11. Tabulated potential: (-tab). The potential is defined numerically in TABBND (see Section Setting up Tabulated Intramolecular Force-Field Files and Section  The TABBND, TABANG, TABDIH & TABINV Files).

In DL_POLY_4 distance restraints are handled by the routine bonds_forces (and intra_coul called within).

Valence Angle Potentials

Fig. 2 The valence angle and associated vectors

The valence angle potentials describe the bond bending terms between the specified atoms. They should not be confused with the three-body potentials described later, which are defined by atom types rather than indices.

1. Harmonic: (harm)

   (14)\[U(\theta_{jik}) = {k \over 2} (\theta_{jik} - \theta_0)^{2}\]

2. Quartic: (quar)

   (15)\[U(\theta_{jik}) = {k \over 2}(\theta_{jik} - \theta_0)^{2} + {k' \over 3}(\theta_{jik} - \theta_0)^{3} + {k'' \over 4}(\theta_{jik} - \theta_0)^{4}\]

3. Truncated harmonic: ( thrm)

   (16)\[U(\theta_{jik}) = {k\over 2} (\theta_{jik} - \theta_0)^{2} \exp[-(r_{ij}^{8} + r_{ik}^{8}) / \rho^{8}]\]

4. Screened harmonic: ( shrm)

   (17)\[U(\theta_{jik}) = {k\over 2} (\theta_{jik} - \theta_0)^{2} \exp[-(r_{ij} / \rho_1 + r_{ik} / \rho_2)]\]

5. Screened Vessal ¹²⁰: ( bvs1)

   (18)\[\begin{split}\begin{aligned} U(\theta_{jik})=&{k \over 8(\theta_{0}-\pi)^{2}} \left[ (\theta_{0} -\pi)^{2} -(\theta_{jik}-\pi)^{2}\right]^{2} \times \nonumber \\ & \exp[-(r_{ij} / \rho_{1} + r_{ik} / \rho_{2})]\end{aligned}\end{split}\]

6. Truncated Vessal ¹⁰⁴: ( bvs2 )

   (19)\[\begin{split}\begin{aligned} U(\theta_{jik})=& k~(\theta_{jik}-\theta_{0})^{2}~ \big[~\theta_{jik}^a (\theta_{jik}+\theta_{0}-2\pi)^{2} + \nonumber \\ & {a \over 2} \pi^{a-1} (\theta_{0}-\pi)^{3}~\big]~~\exp[-(r_{ij}^{8} + r_{ik}^{8}) / \rho^{8}]\end{aligned}\end{split}\]

7. Harmonic cosine: ( hcos)

   (20)\[U(\theta_{jik}) = {k\over 2}(\cos(\theta_{jik}) -\cos(\theta_{0}))^{2}\]

8. Cosine: ( cos)

   (21)\[U(\theta_{jik}) = A~[1 + \cos(m~\theta_{jik}-\delta)]\]

9. MM3 stretch-bend ⁵: ( mmsb)

   (22)\[U(\theta_{jik}) = A~(\theta_{jik}-\theta_{0})~(r_{ij} - r_{ij}^{o})~(r_{ik} - r_{ik}^{o})\]

10. Compass stretch-stretch ¹⁰⁷: ( stst)

    (23)\[U(\theta_{jik}) = A~(r_{ij} - r_{ij}^{o})~(r_{ik} - r_{ik}^{o})\]

11. Compass stretch-bend ¹⁰⁷: ( stbe)

    (24)\[U(\theta_{jik}) = A~(\theta_{jik}-\theta_{0})~(r_{ij} - r_{ij}^{o})\]

12. Compass all terms ¹⁰⁷: ( cmps)

    (25)\[U(\theta_{jik}) = A~(r_{ij} - r_{ij}^{o})~(r_{ik} - r_{ik}^{o}) + (\theta_{jik}-\theta_{0})~[B~(r_{ij} - r_{ij}^{o})+C~(r_{ik} - r_{ik}^{o})]\]

13. MM3 angle bend term ⁵: ( m3ab)

    (26)\[\begin{split}\begin{aligned} U(\theta_{jik}) = k~(\theta_{jik}-\theta_{0})^{2}[1 -& 1.4 \cdot 10^{-2}(\theta_{jik}-\theta_{0})^{~}+5.6 \cdot 10^{-5}(\theta_{jik}-\theta_{0})^{2} \nonumber \\ -& 7.0 \cdot 10^{-7}(\theta_{jik}-\theta_{0})^{3}+2.2 \cdot 10^{-8}(\theta_{jik}-\theta_{0})^{4}]\end{aligned}\end{split}\]

14. KKY ⁵⁰: ( kky)

    (27)\[\begin{split}\begin{aligned} U(\theta_{jik}) =& 2~f_{k}~\sqrt{K_{ij} \cdot K_{ik}}~\sin^{2}\left[(\theta_{jik}-\theta_{0})\right] \nonumber \\ K_{ij} =& \frac{1}{\exp\left[g_{r}(r_{ij}-r_{o})\right]+1}\end{aligned}\end{split}\]

15. Tabulated potential: ( tab). The potential is defined numerically in TABANG (see Section Setting up Tabulated Intramolecular Force-Field Files and Section  The TABBND, TABANG, TABDIH & TABINV Files).

In these formulae \(\theta_{jik}\) is the angle between bond vectors \(\underline{r}_{ij}\) and \(\underline{r}_{ik}\):

\[\theta_{jik}=cos^{-1}\left\{\frac{\underline{r}_{ij}\cdot\underline{r}_{ik}} {r_{ij}r_{ik}}\right\}~~.\]

In DL_POLY_4 the most general form for the valence angle potentials can be written as:

\[U(\theta_{jik},r_{ij},r_{ik}) = A(\theta_{jik})~S(r_{ij})~S(r_{ik})~S(r_{ik})~~,\]

where \(A(\theta)\) is a purely angular function and \(S(r)\) is a screening or truncation function. All the function arguments are scalars. With this reduction the force on an atom derived from the valence angle potential is given by:

\[f_{\ell}^{\alpha} = -\frac{\partial}{\partial r_{\ell}^{\alpha}}U(\theta_{jik},r_{ij},r_{ik},r_{jk})~~,\]

with atomic label \(\ell\) being one of \(i,j,k\) and \(\alpha\) indicating the \(x,y,z\) component. The derivative is

\[\begin{split}\begin{aligned} -\frac{\partial}{\partial r_{\ell}^{\alpha}}U(\theta_{jik},r_{ij},r_{ik},r_{jk})=& -S(r_{ij})S(r_{ik})S(r_{jk})\frac{\partial}{\partial r_{\ell}^{\alpha}}A(\theta_{jik}) \nonumber \\ & - A(\theta_{jik})S(r_{ik})S(r_{jk})(\delta_{\ell j}-\delta_{\ell i}) \frac{r_{ij}^{\alpha}}{r_{ij}} \frac{\partial}{\partial r_{ij}}S(r_{ij}) \nonumber \\ & - A(\theta_{jik})S(r_{ij})S(r_{jk})(\delta_{\ell k}-\delta_{\ell i}) \frac{r_{ik}^{\alpha}}{r_{ik}} \frac{\partial}{\partial r_{ik}}S(r_{ik}) \nonumber \\ & - A(\theta_{jik})S(r_{ij})S(r_{ik})(\delta_{\ell k}-\delta_{\ell j}) \frac{r_{jk}^{\alpha}}{r_{jk}} \frac{\partial}{\partial r_{jk}}S(r_{jk})~~,\end{aligned}\end{split}\]

with \(\delta_{ab}=1\) if \(a=b\) and \(\delta_{ab}=0\) if \(a\ne b\) . In the absence of screening terms \(S(r)\), this formula reduces to:

\[-\frac{\partial}{\partial r_{\ell}^{\alpha}}U(\theta_{jik},r_{ij},r_{ik},r_{jk}) = -\frac{\partial}{\partial r_{\ell}^{\alpha}}A(\theta_{jik})~~.\]

The derivative of the angular function is

\[-\frac{\partial}{\partial r_{\ell}^{\alpha}}A(\theta_{jik}) = \left\{\frac{1}{\sin(\theta_{jik})}\right\} \frac{\partial}{\partial \theta_{jik}}A(\theta_{jik}) \frac{\partial}{\partial r_{\ell}^{\alpha}}\left\{ \frac{\underline{r}_{ij}\cdot\underline{r}_{ik}}{r_{ij}r_{ik}}\right\}~~,\]

with

\[\begin{split}\begin{aligned} \frac{\partial}{\partial r_{\ell}^{\alpha}}\left\{ \frac{\underline{r}_{ij}\cdot\underline{r}_{ik}}{r_{ij}r_{ik}}\right\}=& (\delta_{\ell j}-\delta_{\ell i})\frac{r_{ik}^{\alpha}}{r_{ij}r_{ik}}+ (\delta_{\ell k}- \delta_{\ell i})\frac{r_{ij}^{\alpha}}{r_{ij}r_{ik}}-\nonumber \\ & \cos(\theta_{jik}) \left\{(\delta_{\ell j}-\delta_{\ell i})\frac{r_{ij}^{\alpha}}{r_{ij}^{2}}+ (\delta_{\ell k}-\delta_{\ell i})\frac{r_{ik}^{\alpha}}{r_{ik}^{2}}\right\}~~.\end{aligned}\end{split}\]

The atomic forces are then completely specified by the derivatives of the particular functions \(A(\theta)\) and \(S(r)\) .

The contribution to be added to the atomic virial is given by

\[{\cal W} = -(\underline{r}_{ij} \cdot \underline{f}_{j} + \underline{r}_{ik} \cdot \underline{f}_{k})~~.\]

It is worth noting that in the absence of screening terms S(r), the virial is zero ⁹⁷.

The contribution to be added to the atomic stress tensor is given by

(28)\[\sigma^{\alpha \beta} = r_{ij}^{\alpha} f_{j}^{\beta} + r_{ik}^{\alpha} f_{k}^{\beta}\]

and the stress tensor is symmetric.

In DL_POLY_4 valence forces are handled by the routine angles_forces.

Angular Restraints

In DL_POLY_4 angle restraints, in which the angle subtended by a triplet of atoms, is maintained around some preset value \(\theta_{0}\) is handled as a special case of angle potentials. As a consequence angle restraints may be applied only between atoms in the same molecule. Unlike with application of the “pure” angle potentials, the electrostatic and van der Waals interactions between the pair of atoms are still evaluated when distance restraints are applied. All the potential forms of the previous section are available as angular restraints, although they have different key words:

1. Harmonic: (-hrm)

2. Quartic: (-qur)

3. Truncated harmonic: (-thm)

4. Screened harmonic: (-shm)

5. Screened Vessal ¹²⁰: (-bv1)

6. Truncated Vessal ¹⁰⁴: (-bv2)

7. Harmonic cosine: (-hcs)

8. Cosine: (-cos)

9. MM3 stretch-bend ⁵: (-msb)

10. Compass stretch-stretch ¹⁰⁷: (-sts)

11. Compass stretch-bend ¹⁰⁷: (-stb)

12. Compass all terms ¹⁰⁷: (-cmp)

13. MM3 angle bend ⁵: (-m3a)

14. KKY ⁵⁰: (-kky)

15. Tabulated potential: ( -tab). The potential is defined numerically in TABANG (see Section Setting up Tabulated Intramolecular Force-Field Files and Section  The TABBND, TABANG, TABDIH & TABINV Files).

In DL_POLY_4 angular restraints are handled by the routine angles_forces.

Dihedral Angle Potentials

Fig. 3 The dihedral angle and associated vectors

The dihedral angle potentials describe the interaction arising from torsional forces in molecules. (They are sometimes referred to as torsion potentials.) They require the specification of four atomic positions. The potential functions available in DL_POLY_4 are as follows:

1. Cosine potential: (cos)

   (29)\[U(\phi_{ijkn}) = A~\left [ 1 + \cos(m\phi_{ijkn} - \delta)\right]\]

2. Harmonic: (harm)

   (30)\[U(\phi_{ijkn}) = {k \over 2}~(\phi_{ijkn} - \phi_{0})^{2}\]

3. Harmonic cosine: (hcos)

   (31)\[U(\phi_{ijkn}) = {k \over 2}~(\cos(\phi_{ijkn}) - \cos(\phi_{0}))^{2}\]

4. Triple cosine: (cos3)

   (32)\[U(\phi) = {1 \over 2}~\{A_{1}~(1+\cos(\phi)) + A_{2}~(1-\cos(2\phi)) + A_{3}~(1+\cos(3\phi))\}\]

5. Ryckaert-Bellemans ⁸⁴ with fixed constants a-f: (ryck)

   (33)\[\begin{split}\begin{aligned} U(\phi) = A~\{~a~+&~b~\cos(\phi)+c~\cos^{2}(\phi)+d~\cos^{3}(\phi) \nonumber \\ +&~e~\cos^{4}(\phi)+f~\cos^{5}(\phi)~\} \end{aligned}\end{split}\]

6. Fluorinated Ryckaert-Bellemans ⁸⁸ with fixed constants a-h: (rbf)

   (34)\[\begin{split}\begin{aligned} U(\phi) = A~\{~a~+&~b~\cos(\phi)+c~\cos^{2}(\phi)+d~\cos^{3}(\phi)+e~\cos^{4}(\phi) \nonumber \\ +&~f~\cos^{5}(\phi)+~g~\exp(-h(\phi-\pi)^{2}))~\} \end{aligned}\end{split}\]

7. OPLS torsion potential: (opls)

   (35)\[U(\phi)=A_{0}+\frac{1}{2} \left\{ A_{1}~(1+\cos(\phi))+A_{2}~(1-\cos(2\phi))+A_{3}~(1+\cos(3\phi)) \right\}\]

8. Tabulated potential: (tab). The potential is defined numerically in TABDIH (see Section Setting up Tabulated Intramolecular Force-Field Files and Section  The TABBND, TABANG, TABDIH & TABINV Files).

In these formulae \(\phi_{ijkn}\) is the dihedral angle defined by

\[\phi_{ijkn}=\cos^{-1}\{B(\underline{r}_{ij},\underline{r}_{jk},\underline{r}_{kn})\}~~,\]

with

\[B(\underline{r}_{ij},\underline{r}_{jk},\underline{r}_{kn})= \left\{\frac{ (\underline{r}_{ij}\times\underline{r}_{jk})\cdot(\underline{r}_{jk}\times\underline{r}_{kn})} {|\underline{r}_{ij}\times\underline{r}_{jk}||\underline{r}_{jk}\times\underline{r}_{kn}|} \right\}~~.\]

With this definition, the sign of the dihedral angle is positive if the vector product

\[(\underline{r}_{ij} \times \underline{r}_{jk}) \times (\underline{r}_{jk} \times \underline{r}_{kn})\]

is in the same direction as the bond vector \(\underline{r}_{jk}\) and negative if in the opposite direction.

The force on an atom arising from the dihedral potential is given by

(36)\[f_{\ell}^{\alpha}=-\frac{\partial}{\partial r_{\ell}^{\alpha}}U(\phi_{ijkn})~~,\]

with \(\ell\) being one of \(i,j,k,n\) and \(\alpha\) one of \(x,y,z\). This may be expanded into

\[-\frac{\partial}{\partial r_{\ell}^{\alpha}}U(\phi_{ijkn})= \left\{\frac{1}{\sin(\phi_{ijkn})}\right\} \frac{\partial}{\partial \phi_{ijkn}}U(\phi_{ijkn}) \frac{\partial}{\partial r_{\ell}^{\alpha}} B(\underline{r}_{ij},\underline{r}_{jk},\underline{r}_{kn})~~.\]

The derivative of the function \(B(\underline{r}_{ij},\underline{r}_{jk},\underline{r}_{kn})\) is

\[\begin{split}\begin{aligned} &\frac{\partial}{\partial r_{\ell}^{\alpha}} B(\underline{r}_{ij},\underline{r}_{jk},\underline{r}_{kn})~=~ \frac{1}{|\underline{r}_{ij}\times\underline{r}_{jk}||\underline{r}_{jk}\times\underline{r}_{kn}|} \frac{\partial}{\partial r_{\ell}^{\alpha}} \{(\underline{r}_{ij}\times\underline{r}_{jk})\cdot(\underline{r}_{jk}\times\underline{r}_{kn})\}~-~\nonumber \\ & \phantom{x} \frac{\cos(\phi_{ijkn})}{2}\left\{ \frac{1}{|\underline{r}_{ij}\times\underline{r}_{jk}|^{2}} \frac{\partial}{\partial r_{\ell}^{\alpha}} |\underline{r}_{ij}\times\underline{r}_{jk}|^{2}+ \frac{1}{|\underline{r}_{jk}\times\underline{r}_{kn}|^{2}} \frac{\partial}{\partial r_{\ell}^{\alpha}} |\underline{r}_{jk}\times\underline{r}_{kn}|^{2} \right\}~, \end{aligned}\end{split}\]

with

\[\begin{split}\begin{aligned} \frac{\partial}{\partial r_{\ell}^{\alpha}} \{(\underline{r}_{ij}\times\underline{r}_{jk})\cdot(\underline{r}_{jk}\times\underline{r}_{kn})\}= & r_{ij}^{\alpha}([\underline{r}_{jk}\underline{r}_{jk}]_{\alpha}(\delta_{\ell k}-\delta_{\ell n})+ [\underline{r}_{jk}\underline{r}_{kn}]_{\alpha}(\delta_{\ell k}-\delta_{\ell j}))+ \nonumber \\ \phantom{\frac{\partial}{\partial r_{\ell}^{\alpha}}}& r_{jk}^{\alpha}([\underline{r}_{ij}\underline{r}_{jk}]_{\alpha} (\delta_{\ell n}-\delta_{\ell k})+ [\underline{r}_{jk}\underline{r}_{kn}]_{\alpha} (\delta_{\ell j}-\delta_{\ell i}))+ \nonumber \\ \phantom{\frac{\partial}{\partial r_{\ell}^{\alpha}}} & r_{kn}^{\alpha}([\underline{r}_{ij}\underline{r}_{jk}]_{\alpha}(\delta_{\ell k}-\delta_{\ell j})+ [\underline{r}_{jk}\underline{r}_{jk}]_{\alpha}(\delta_{\ell i}-\delta_{\ell j}))+ \nonumber \\ \phantom{\frac{\partial}{\partial r_{\ell}^{\alpha}}} & 2r_{jk}^{\alpha}[\underline{r}_{ij}\underline{r}_{kn}]_{\alpha}(\delta_{\ell j}-\delta_{\ell k})~~,\end{aligned}\end{split}\]

\[\begin{split}\begin{aligned} \frac{\partial}{\partial r_{\ell}^{\alpha}} |\underline{r}_{ij}\times\underline{r}_{jk}|^{2}= & 2r_{ij}^{\alpha}([\underline{r}_{jk}\underline{r}_{jk}]_{\alpha}(\delta_{\ell j}-\delta_{\ell i})+ [\underline{r}_{ij}\underline{r}_{jk}]_{\alpha}(\delta_{\ell j}-\delta_{\ell k}))+\nonumber \\ & 2r_{jk}^{\alpha}([\underline{r}_{ij}\underline{r}_{ij}]_{\alpha} (\delta_{\ell k}-\delta_{\ell j}) +[\underline{r}_{ij}\underline{r}_{jk}]_{\alpha} (\delta_{\ell i}-\delta_{\ell j}))~~,\end{aligned}\end{split}\]

\[\begin{split}\begin{aligned} \frac{\partial}{\partial r_{\ell}^{\alpha}} |\underline{r}_{jk}\times\underline{r}_{kn}|^{2} = & 2r_{kn}^{\alpha}([\underline{r}_{jk}\underline{r}_{jk}]_{\alpha}(\delta_{\ell n}-\delta_{\ell k})+ [\underline{r}_{jk}\underline{r}_{kn}]_{\alpha}(\delta_{\ell j}-\delta_{\ell k}))+\nonumber \\ & 2r_{jk}^{\alpha}([\underline{r}_{kn}\underline{r}_{kn}]_{\alpha} (\delta_{\ell k}-\delta_{\ell j})+ [\underline{r}_{jk}\underline{r}_{kn}]_{\alpha}(\delta_{\ell k}-\delta_{\ell n}))~~.\end{aligned}\end{split}\]

Where we have used the following definition:

\[_{\alpha} = \sum_{\beta}(1-\delta_{\alpha\beta})a^{\beta}b^{\beta}~~.\]

Formally, the contribution to be added to the atomic virial is given by

\[{\cal W}=-\sum_{i=1}^{4} \underline{r}_{i} \cdot \underline{f}_{i}~~.\]

However, it is possible to show (by tedious algebra using the above formulae, or more elegantly by thermodynamic arguments ⁹⁷,) that the dihedral makes no contribution to the atomic virial.

The contribution to be added to the atomic stress tensor is given by

(37)\[\begin{split}\begin{aligned} \sigma^{\alpha \beta} = &r_{ij}^{\alpha}p_{i}^{\beta}+ r_{jk}^{\alpha}p_{jk}^{\beta}+r_{kn}^{\alpha}p_{n}^{\beta} \\ & -\frac{\cos(\phi_{ijkn})}{2}\left \{r_{ij}^{\alpha}g_{i}^{\beta}+ r_{jk}^{\alpha}g_{k}^{\beta}+r_{jk}^{\alpha}h_{j}^{\beta}+r_{kn}^{\alpha}h_{n}^{\beta}\right\} ~~,\nonumber\end{aligned}\end{split}\]

with

\[\begin{split}\begin{aligned} p_{i}^{\alpha} =& (r_{jk}^{\alpha}[\underline{r}_{jk}\underline{r}_{kn}]_{\alpha}- r_{kn}^{\alpha}[\underline{r}_{jk}\underline{r}_{jk}]_{\alpha})/ (|\underline{r}_{ij}\times\underline{r}_{jk}||\underline{r}_{jk}\times\underline{r}_{kn}|) \nonumber \\ p_{n}^{\alpha} =& (r_{jk}^{\alpha}[\underline{r}_{ij}\underline{r}_{jk}]_{\alpha}- r_{ij}^{\alpha}[\underline{r}_{jk}\underline{r}_{jk}]_{\alpha})/ (|\underline{r}_{ij}\times\underline{r}_{jk}||\underline{r}_{jk}\times\underline{r}_{kn}|) \nonumber \\ p_{jk}^{\alpha} =& (r_{ij}^{\alpha}[\underline{r}_{jk}\underline{r}_{kn}]_{\alpha}+ r_{kn}^{\alpha}[\underline{r}_{ij}\underline{r}_{jk}]_{\alpha}- 2r_{jk}^{\alpha}[\underline{r}_{ij}\underline{r}_{kn}]_{\alpha})/ (|\underline{r}_{ij}\times\underline{r}_{jk}||\underline{r}_{jk}\times\underline{r}_{kn}|) \nonumber \\ g_{i}^{\alpha} =& 2(r_{ij}^{\alpha}[\underline{r}_{jk}\underline{r}_{jk}]_{\alpha}- r_{jk}^{\alpha}[\underline{r}_{ij}\underline{r}_{jk}]_{\alpha})/ |\underline{r}_{ij}\times\underline{r}_{jk}|^{2} \\ g_{k}^{\alpha} =& 2(r_{jk}^{\alpha}[\underline{r}_{ij}\underline{r}_{ij}]_{\alpha}- r_{ij}^{\alpha}[\underline{r}_{ij}\underline{r}_{jk}]_{\alpha})/ |\underline{r}_{ij}\times\underline{r}_{jk}|^{2} \nonumber \\ h_{j}^{\alpha} =& 2(r_{jk}^{\alpha}[\underline{r}_{kn}\underline{r}_{kn}]_{\alpha}- r_{kn}^{\alpha}[\underline{r}_{jk}\underline{r}_{kn}]_{\alpha})/ |\underline{r}_{jk}\times\underline{r}_{kn}|^{2} \nonumber \\ h_{n}^{\alpha} =& 2(r_{kn}^{\alpha}[\underline{r}_{kn}\underline{r}_{kn}]_{\alpha}- r_{jk}^{\alpha}[\underline{r}_{jk}\underline{r}_{kn}]_{\alpha})/ |\underline{r}_{jk}\times\underline{r}_{kn}|^{2}~~. \nonumber\end{aligned}\end{split}\]

The sum of the diagonal elements of the stress tensor is zero (since the virial is zero) and the matrix is symmetric.

Lastly, it should be noted that the above description does not take into account the possible inclusion of distance-dependent 1-4 interactions, as permitted by some force field`s. Such interactions are permissible in DL_POLY_4 and are described in the section on pair potentials below. DL_POLY_4 also permits scaling of the 1-4 .:index:`van der Waals and Coulomb interactions by a numerical factor (see Table (6)).

Note

Scaling is abandoned when the 1-4 members are also 1-3 members in a valence angle intercation (1-4 checks are performed in dihedrals_14_check routine). 1-4 interactions do, of course, contribute to the atomic virial.

In DL_POLY_4 dihedral forces are handled by the routine dihedrals_forces (and intra_coul and dihedrals_14_vdw called within).

Improper Dihedral Angle Potentials

Improper dihedrals are used to restrict the geometry of molecules and as such need not have a simple relation to conventional chemical bonding. DL_POLY_4 makes no distinction between dihedral and improper dihedral angle functions (both are calculated by the same subroutines) and all the comments made in the preceding section apply.

An important example of the use of the improper dihedral is to conserve the structure of chiral centres in molecules modelled by united-atom centres. For example \(\alpha\)-amino acids such as alanine (CH:math:_{3}CH(NH\(_{2}\))COOH), in which it is common to represent the CH\(_{3}\) and CH groups as single centres. Conservation of the chirality of the \(\alpha\) carbon is achieved by defining a harmonic improper dihedral angle potential with an equilibrium angle of 35.264\(^{o}\). The angle is defined by vectors \(\underline{r}_{12}\), \(\underline{r}_{23}\) and \(\underline{r}_{34}\), where the atoms 1,2,3 and 4 are shown in the following figure. The figure defines the D and L enantiomers consistent with the international (IUPAC) convention. When defining the dihedral, the atom indices are entered in DL_POLY_4 in the order 1-2-3-4.

Fig. 4 The L and D enantiomers and defining vectors

In DL_POLY_4 improper dihedral forces are handled by the routine dihedrals_forces.

Torsional Restraints

In DL_POLY_4 the torsional restraints, in which the dihedral angle as defined by a quadruplet of atoms, is maintained around some preset value \(\phi_{0}\) is handled as a special case of dihedral potential. As a consequence angle restraints may be applied only between atoms in the same molecule. Unlike with application of the “pure” dihedral potentials, the electrostatic and van der Waals interactions between the pair of atoms are still evaluated when distance restraints are applied. All the potential forms of the previous section are available as torsional restraints, although they have different key words:

1. Cosine potential: (-cos)

2. Harmonic: (-hrm)

3. Harmonic cosine: (-hcs)

4. Triple cosine: (-cs3)

5. Ryckaert-Bellemans ⁸⁴ with fixed constants a-f: (-rck)

6. Fluorinated Ryckaert-Bellemans ⁸⁸ with fixed constants a-h: (-rbf)

7. OPLS torsion potential: (-opl)

8. Tabulated potential: (-tab). The potential is defined numerically in TABDIH (see Section Setting up Tabulated Intramolecular Force-Field Files and Section  The TABBND, TABANG, TABDIH & TABINV Files).

In DL_POLY_4 torsional restraints are handled by the routine dihedrals_forces.

Inversion Angle Potentials

Fig. 5 The inversion angle and associated vectors

The inversion angle potentials describe the interaction arising from a particular geometry of three atoms around a central atom. The best known example of this is the arrangement of hydrogen atoms around nitrogen in ammonia to form a trigonal pyramid. The hydrogens can ‘flip’ like an inverting umbrella to an alternative structure, which in this case is identical, but in principle causes a change in chirality. The force restraining the ammonia to one structure can be described as an inversion potential (though it is usually augmented by valence angle potentials also). The inversion angle is defined in the figure above. It resembles a dihedral potential in that it requires the specification of four atomic positions.

Note

The inversion angle potential is a sum of the three possible inversion angle terms

The potential functions available in DL_POLY_4 are as follows:

1. Harmonic: (harm)

   (38)\[U(\phi_{ijkn}) = {k \over 2}~(\phi_{ijkn} - \phi_{0})^{2}\]

2. Harmonic cosine: (hcos)

   (39)\[U(\phi_{ijkn}) = {k \over 2}~(\cos(\phi_{ijkn}) - \cos(\phi_{0}))^{2}\]

3. Planar potential: (plan)

   (40)\[U(\phi_{ijkn}) = A~\left[ 1 - \cos(\phi_{ijkn}) \right]\]

4. Extended planar potential: (xpln)

   (41)\[U(\phi_{ijkn}) = {k \over 2}~\left[ 1 - \cos(m~\phi_{ijkn} - \phi_{0}) \right]\]

5. Tabulated potential: (tab). The potential is defined numerically in TABINV (see Section Setting up Tabulated Intramolecular Force-Field Files and Section  The TABBND, TABANG, TABDIH & TABINV Files).

In these formulae \(\phi_{ijkn}\) is the inversion angle defined by

\[\phi_{ijkn}=\cos^{-1} \left\{ \frac{\underline{r}_{ij}\cdot\underline{w}_{kn}}{r_{ij}w_{kn}} \right\}~~,\]

with

\[\underline{w}_{kn}=(\underline{r}_{ij}\cdot\underline{\hat{u}}_{kn})\underline{\hat{u}}_{kn}+ (\underline{r}_{ij}\cdot\underline{\hat{v}}_{kn})\underline{\hat{v}}_{kn}\]

and the unit vectors

\[\begin{split}\begin{aligned} \underline{\hat{u}}_{kn}=&(\underline{\hat{r}}_{ik}+\underline{\hat{r}}_{in})/ |\underline{\hat{r}}_{ik}+\underline{\hat{r}}_{in}| \nonumber \\ \underline{\hat{v}}_{kn}=&(\underline{\hat{r}}_{ik}-\underline{\hat{r}}_{in})/ |\underline{\hat{r}}_{ik}-\underline{\hat{r}}_{in}|~~.\end{aligned}\end{split}\]

As usual, \(\underline{r}_{ij}=\underline{r}_{j}-\underline{r}_{i}\) etc. and the hat \(\underline{\hat{r}}\) indicates a unit vector in the direction of \(\underline{r}\). The total inversion potential requires the calculation of three such angles, the formula being derived from the above using the cyclic permutation of the indices \(j \rightarrow k \rightarrow n \rightarrow j\) etc.

Equivalently, the angle \(\phi_{ijkn}\) may be written as

\[\phi_{ijkn} = \cos^{-1} \left\{ \frac{ [ (\underline{r}_{ij}\cdot\underline{\hat{u}}_{kn})^{2} + (\underline{r}_{ij}\cdot\underline{\hat{v}}_{kn})^{2} ]^{1/2}}{r_{ij}} \right\}~~.\]

Formally, the force on an atom arising from the inversion potential is given by

(42)\[f_{\ell}^{\alpha} = -\frac{\partial}{\partial r_{\ell}^{\alpha}}U(\phi_{ijkn})~~,\]

with \(\ell\) being one of \(i,j,k,n\) and \(\alpha\) one of \(x,y,z\). This may be expanded into

\[\begin{split}\begin{aligned} -\frac{\partial}{\partial r_{\ell}^{\alpha}}U(\phi_{ijkn})=& \left\{ \frac{1}{\sin(\phi_{ijkn})} \right\} \frac{\partial}{\partial \phi_{ijkn}}U(\phi_{ijkn})\times \nonumber \\ & \frac{\partial}{\partial r_{\ell}^{\alpha}} \left\{ \frac{[(\underline{r}_{ij}\cdot\underline{\hat{u}}_{kn})^{2} +(\underline{r}_{ij}\cdot\underline{\hat{v}}_{kn})^{2}]^{1/2}} {r_{ij}} \right\}~~.\end{aligned}\end{split}\]

Following through, the (extremely tedious!) differentiation gives the result:

\[\begin{split}\begin{aligned} f_{\ell}^{\alpha} =& \left\{\frac{1}{\sin(\phi_{ijkn})}\right\} \frac{\partial}{\partial \phi_{ijkn}}U(\phi_{ijkn})\times \\ & \left\{-(\delta_{\ell j}-\delta_{\ell i})\frac{\cos(\phi_{ijkn})} {r_{ij}^{2}}r_{ij}^{\alpha} +\frac{1}{r_{ij}w_{kn}}\left [ (\delta_{\ell j}-\delta_{\ell i}) \{(\underline{r}_{ij}\cdot\underline{\hat{u}}_{kn})\hat{u}_{kn}^{\alpha}+ (\underline{r}_{ij}\cdot\underline{\hat{v}}_{kn})\hat{v}_{kn}^{\alpha}\} \phantom{\left\{\frac{a_{a}^{a}}{a_{a}^{a}}\right\}} \right. \right. \nonumber \\ & + (\delta_{\ell k}-\delta_{\ell i}) \frac{\underline{r}_{ij}\cdot\underline{\hat{u}}_{kn}}{u_{kn}r_{ik}}\left \{ r_{ij}^{\alpha}-(\underline{r}_{ij}\cdot\underline{\hat{u}}_{kn})\hat{u}_{kn}^{\alpha} -(\underline{r}_{ij}\cdot\underline{r}_{ik}-(\underline{r}_{ij}\cdot\underline{\hat{u}}_{kn}) (\underline{r}_{ik}\cdot\underline{\hat{u}}_{kn}))\frac{r_{ik}^{\alpha}}{r_{ik}^{2}} \right \} \nonumber \\ & + (\delta_{\ell k}-\delta_{\ell i}) \frac{\underline{r}_{ij}\cdot\underline{\hat{v}}_{kn}}{v_{kn}r_{ik}}\left \{ r_{ij}^{\alpha}-(\underline{r}_{ij}\cdot\underline{\hat{v}}_{kn})\hat{v}_{kn}^{\alpha} -(\underline{r}_{ij}\cdot\underline{r}_{ik}-(\underline{r}_{ij}\cdot\underline{\hat{v}}_{kn}) (\underline{r}_{ik}\cdot\underline{\hat{v}}_{kn}))\frac{r_{ik}^{\alpha}}{r_{ik}^{2}} \right \} \nonumber \\ &+ (\delta_{\ell n}-\delta_{\ell i}) \frac{\underline{r}_{ij}\cdot\underline{\hat{u}}_{kn}}{u_{kn}r_{in}}\left \{ r_{ij}^{\alpha}-(\underline{r}_{ij}\cdot\underline{\hat{u}}_{kn})\hat{u}_{kn}^{\alpha} -(\underline{r}_{ij}\cdot\underline{r}_{in}-(\underline{r}_{ij}\cdot\underline{\hat{u}}_{kn}) (\underline{r}_{in}\cdot\underline{\hat{u}}_{kn}))\frac{r_{in}^{\alpha}}{r_{in}^{2}} \right \} \nonumber \\ & \left . \left .- (\delta_{\ell n}-\delta_{\ell i}) \frac{\underline{r}_{ij}\cdot\underline{\hat{v}}_{kn}}{v_{kn}r_{in}}\left \{ r_{ij}^{\alpha}-(\underline{r}_{ij}\cdot\underline{\hat{v}}_{kn})\hat{v}_{kn}^{\alpha} -(\underline{r}_{ij}\cdot\underline{r}_{in}-(\underline{r}_{ij}\cdot\underline{\hat{v}}_{kn}) (\underline{r}_{in}\cdot\underline{\hat{v}}_{kn}))\frac{r_{in}^{\alpha}}{r_{in}^{2}} \right \} \right ] \right \}~~. \nonumber\end{aligned}\end{split}\]

This general formula applies to all atoms \(\ell=i,j,k,n\). It must be remembered however, that these formulae apply to just one of the three contributing terms (i.e. one angle \(\phi\)) of the full inversion potential: specifically the inversion angle pertaining to the out-of-plane vector \(\underline{r}_{ij}\) . The contributions arising from the other vectors \(\underline{r}_{ik}\) and \(\underline{r}_{in}\) are obtained by the cyclic permutation of the indices in the manner described above. All these force contributions must be added to the final atomic forces.

Formally, the contribution to be added to the atomic virial is given by

\[{\cal W} = -\sum_{i=1}^{4}\underline{r}_{i} \cdot \underline{f}_{i}~~.\]

However, it is possible to show by thermodynamic arguments ( cf ⁹⁷,) or simply from the fact that the sum of forces on atoms j,k and n is equal and opposite to the force on atom i, that the inversion potential makes no contribution to the atomic virial.

If the force components \(f_{\ell}^{\alpha}\) for atoms \(\ell=i,j,k,n\) are calculated using the above formulae, it is easily seen that the contribution to be added to the atomic stress tensor` is given by

(43)\[\sigma^{\alpha \beta} = r_{ij}^{\alpha} f_{j}^{\beta} + r_{ik}^{\alpha} f_{k}^{\beta}+r_{in}^{\alpha}f_{n}^{\beta}~~.\]

The sum of the diagonal elements of the stress tensor is zero (since the virial is zero) and the matrix is symmetric.

In DL_POLY_4 inversion forces are handled by the routine inversions_forces.

The Calcite Four-Body Potential

Fig. 6 The vectors of the calcite potential

This potential ^78,83 is designed to help maintain the planar structure of the carbonate anion \([CO_{3}]^{2-}\) in a similar manner to the planar inversion potential described above. However, it is not an angular potential. It is dependent on the perpendicular displacement (\(u\)) of an atom \(a\) from a plane defined by three other atoms \(b\), \(c\), and \(d\) (see Figure 6) and has the form:

(44)\[U_{abcd}(u)=Au^{2}+Bu^{4}~~, \label{calcite1}\]

where the displacement \(u\) is given by

(45)\[u=\frac{\underline{r}_{ab}\cdot\underline{r}_{bc}\times\underline{r}_{bd}}{|\underline{r}_{bc}\times\underline{r}_{bd}|}~~. \label{calcite2}\]

Vectors \(\underline{r}_{ab}\),\(\underline{r}_{ac}\) and \(\underline{r}_{ad}\) define bonds between the central atom \(a\) and the peripheral atoms \(b\), \(c\) and \(d\). Vectors \(\underline{r}_{bc}\) and \(\underline{r}_{bd}\) define the plane and are related to the bond vectors by

\[\begin{split}\begin{aligned} \underline{r}_{bc}=&\underline{r}_{ac}-\underline{r}_{ab} \nonumber \\ \underline{r}_{bd}=&\underline{r}_{ad}-\underline{r}_{ab}~~.\end{aligned}\end{split}\]

In what follows it is convenient to define the vector product appearing in both the numerator and denominator of equation (45) as the vector \(\underline{w}_{cd}\) vis.

\[\underline{w}_{cd}=\underline{r}_{bc}\times\underline{r}_{bd}~~.\]

We also define the quantity \(\gamma(u)\) as

\[\gamma(u)=-(2Au+4Bu^{3})~~.\]

The forces on the individual atoms due to the calcite potential are then given by

\[\begin{split}\begin{aligned} \underline{f}_{a}=&-\gamma(u)~\hat{\underline{w}}_{cd} \nonumber \\ \underline{f}_{c}=&\phantom{+}\underline{r}_{bd}\times(\underline{r}_{ab}-u\hat{\underline{w}}_{cd})~\gamma(u)/w_{cd} \nonumber \\ \underline{f}_{d}=&-\underline{r}_{bc}\times(\underline{r}_{ab}- u\hat{\underline{w}}_{cd})~\gamma(u)/w_{cd} \\ \underline{f}_{b}=&-(\underline{f}_{a}+\underline{f}_{c}+\underline{f}_{d})~~,\nonumber\end{aligned}\end{split}\]

where \(w_{cd}=|\underline{w}_{cd}|\) and \(\hat{\underline{w}}_{cd}=\underline{w}_{cd}/w_{cd}\). The virial contribution \(\psi_{abcd}(u)\) is given by

\[\psi_{abcd}(u)=2Au^{2}+4Bu^{4}\]

and the stress tensor contribution \(\sigma_{abcd}^{\alpha\beta}(u)\) by

\[\sigma_{abcd}^{\alpha\beta}(u)=\frac{u~\gamma(u)}{w_{cd}^{2}}~w_{cd}^{\alpha}~w_{cd}^{\beta}~.\]

In DL_POLY_4 the calcite forces are handled by the routine inversions_forces, which is a convenient intramolecular four-body force routine. However, it is manifestly not an inversion potential as such.

Inversional Restraints

In DL_POLY_4 the inversional restraints, in which the inversion angle, as defined by a quadruplet of atoms, is maintained around some preset value \(\phi_{0}\), is handled as a special case of inversion potential. As a consequence angle restraints may be applied only between atoms in the same molecule. Unlike with application of the “pure” dihedral potentials, the electrostatic and van der Waals interactions between the pair of atoms are still evaluated when distance restraints are applied. All the potential forms of the previous section are available as torsional restraints, although they have different key words:

1. Harmonic: (-hrm)

2. Harmonic cosine: (-hcs)

3. Planar potential: (-pln)

4. Extended planar potential: (-xpl)

5. Tabulated potential: (-tab). The potential is defined numerically in TABINV (see Section Setting up Tabulated Intramolecular Force-Field Files and Section  The TABBND, TABANG, TABDIH & TABINV Files).

In DL_POLY_4 inversional restraints are handled by the routine inversions_forces.

Tethering Forces

DL_POLY_4 also allows atomic sites to be tethered to a fixed point in space, \(\vec{r_{0}}\), taken as their position at the beginning of the simulation (t = 0). This is also known as position restraining. The specification, which comes as part of the molecular description, requires a tether potential type and the associated interaction parameters.

Note

Firstly, that application of tethering potentials means that the momentum will no longer be a conserved quantity of the simulation. Secondly, in constant pressure simulations, where the MD cell changes size or shape, the tethers’ reference positions are scaled with the cell vectors.

The tethering potential functions available in DL_POLY_4 are as follows:

1. Harmonic: (harm)

   (46)\[U(r_{ij}) = \frac{1}{2}k(r_{i0})^{2}\]

2. Restrained harmonic: ( rhrm)

   (47)\[\begin{split}U(r_{ij}) = \left\{ \begin{array} {l@{\quad:\quad}l} \frac{1}{2}k(r_{i0})^{2} & |r_{i0}| \le r_{c} \\ \frac{1}{2}kr_{c}^{2} + kr_{c}(r_{i0}-r_{c}) & |r_{i0}| > r_{c} \end{array} \right.\end{split}\]

3. Quartic potential: (quar)

   (48)\[U(r_{ij}) = \frac{k}{2}(r_{i0})^{2}+\frac{k'}{3}(r_{i0})^{3}+\frac{k''}{4}(r_{i0})^{4}\]

as in each case \(r_{io}\) is the distance between the atom positions at moment \(t~=~t1\) and \(t~=~0\).

The force on the atom \(i\) arising from a tether potential potential is obtained using the general formula:

\[\underline{f}_{i} = -\frac{1}{{r}_{i0}}\left[ \frac{\partial }{\partial r_{i0}}U(r_{i0})\right]\underline{r}_{i0}~~.\]

The contribution to be added to the atomic virial is given by

\[{\cal W} = \underline{r}_{i0} \cdot \underline{f}_{i}~~.\]

The contribution to be added to the atomic stress tensor is given by

(49)\[\sigma^{\alpha \beta} = -r_{i0}^{\alpha} f_{i}^{\beta}~~,\]

where \(\alpha\) and \(\beta\) indicate the \(x,y,z\) components. The atomic stress tensor derived in this way is symmetric.

In DL_POLY_4 tether forces are handled by the routine tethers_forces.