Long Ranged Electrostatic (coulombic) Potentials

DL_POLY_4 incorporates several techniques for dealing with long-ranged electrostatic potentials [1]. These are as follows:

1. Direct Coulomb sum

2. Force-shifted Coulomb sum

3. Coulomb sum with distance dependent dielectric

4. Reaction field

5. Smoothed Particle Mesh Ewald (SPME)

All of these can be used in conjunction with the shell model technique used to account for ions polarisation.

The SPME technique is restricted to periodic systems only. (Users must exercise care when using pseudo-periodic boundary conditions.) The other techniques can be used with either periodic or non-periodic systems safely, although in the case of the direct Coulomb sum there are likely to be problems with convergence.

DL_POLY_4 will correctly handle the electrostatics of both molecular and atomic species. However, it is assumed that the system is electrically neutral. A warning message is printed if the system is found to be charged, but otherwise the simulation proceeds as normal.

Note

DL_POLY_4 does not use the basic Ewald method, which is an option in , on account of it being too slow for large scale systems. The SPME method is the standard Ewald method in .

Default (Point Charges) Electrostatics

Direct Coulomb Sum

Use of the direct Coulomb sum is sometimes necessary for accurate simulation of isolated (non-periodic) systems. It is not recommended for periodic systems.

The interaction potential for two charged ions is

\[U(r_{ij}) = \frac{1}{4\pi\epsilon_{0}\epsilon}\frac{q_{i}q_{j}}{r_{ij}}~~,\]

with \(q_{\ell}\) the charge on an atom labelled \(\ell\), and \(r_{ij}\) the magnitude of the separation vector \(\underline{r}_{ij}=\underline{r}_{j}-\underline{r}_{i}\) .

The force on an atom \(j\) derived from this force is

\[\underline{f}_{j} = \frac{1}{4\pi\epsilon_{0}\epsilon} \frac{q_{i}q_{j}}{r_{ij}^{3}} \underline{r}_{ij}~~,\]

with the force on atom \(i\) the negative of this.

The contribution to the atomic virial is

\[{\cal W} = -\frac{1}{4\pi\epsilon_{0}\epsilon} \frac{q_{i}q_{j}}{r_{ij}}~~,\]

which is simply the negative of the potential term.

The contribution to be added to the atomic stress tensor is

\[\sigma^{\alpha \beta} = r_{ij}^{\alpha} f_{j}^{\beta}~~,\]

where \(\alpha,\beta\) are \(x,y,z\) components. The atomic stress tensor is symmetric.

In DL_POLY_4 these forces are handled by the subroutine coul_cp_forces.

Force-Shifted Coulomb Sum

This form of the Coulomb sum has the advantage that it drastically reduces the range of electrostatic interactions, without giving rise to a violent step in the potential energy at the cutoff. Its main use is for preliminary preparation of systems and it is not recommended for realistic models.

The form of the simple truncated and shifted potential function is

\[U(r_{ij}) = \frac{q_{i}q_{j}}{4\pi\epsilon_{0}\epsilon} \left\{\frac{1}{r_{ij}} - \frac{1}{r_{\rm cut}}\right\}~~,\]

with \(q_{\ell}\) the charge on an atom labelled \(\ell\), \(r_{\rm cut}\) the cutoff radius and \(r_{ij}\) the magnitude of the separation vector \(\underline{r}_{ij}=\underline{r}_{j}-\underline{r}_{i}\) .

A further refinement of this approach is to truncate the \(1/r\) potential at \(r_{\rm cut}\) and add a linear term to the potential in order to make both the energy and the force zero at the cutoff. This removes the heating effects that arise from the discontinuity in the forces at the cutoff in the simple truncated and shifted potential (the formula above). (The physics of this potential, however, is little better. It is only recommended for very crude structure optimizations.)

The force-shifted potential is thus

\[U(r_{ij}) = {q_{i} q_{j} \over 4\pi\epsilon_{0}\epsilon} \left[ \left\{\frac{1}{r_{ij}}+\frac{1}{r_{\rm cut}^{2}}~r_{ij}\right\} - \left\{\frac{1}{r_{\rm cut}}+\frac{1}{r_{\rm cut}^{2}}~r_{\rm cut}\right\} \right] = {q_{i} q_{j} \over 4\pi\epsilon_{0}\epsilon} \left[ \frac{1}{r_{ij}} + \frac{r_{ij}}{r_{\rm cut}^{2}} - \frac{2}{r_{\rm cut}} \right]~~,\]

with the force on an atom \(j\) given by

\[\underline{f}_{j} = \frac{q_{i}q_{j}}{4\pi\epsilon_{0}\epsilon} \left[ \frac{1} {r_{ij}^{3}} - \frac{1}{r_{ij}r_{\rm cut}^{2}} \right] \underline{r}_{ij}~~,\]

with the force on atom \(i\) the negative of this.

The force-shifted Coulomb potential can be elegantly extended to emulate long-range ordering by including distance depending damping function \({\rm erfc}(\alpha~r_{ij})\) (identical to that seen in the real-space portion of the Ewald sum) and thus mirror the effective charge screening ³¹ as shown below

\[\begin{split}\begin{aligned} U(r_{ij}) = {q_{i} q_{j} \over 4\pi\epsilon_{0}\epsilon} \left[ \left\{ \frac{{\rm erfc}(\alpha~r_{ij})}{r_{ij}} + \left(\frac{{\rm erfc}(\alpha~r_{\rm cut})}{r_{\rm cut}^{2}}+ \frac{2\alpha}{\sqrt{\pi}}~\frac{\exp(-\alpha^{2}~r_{\rm cut}^{2})}{r_{\rm cut}}\right)~r_{ij}\right\} - \phantom{xx} \right. \nonumber \\ \left. \left\{ \frac{{\rm erfc}(\alpha~r_{\rm cut})}{r_{\rm cut}} + \left(\frac{{\rm erfc}(\alpha~r_{\rm cut})}{r_{\rm cut}^{2}}+ \frac{2\alpha}{\sqrt{\pi}}~\frac{\exp(-\alpha^{2}~r_{\rm cut}^{2})}{r_{\rm cut}}\right)~r_{\rm cut}\right\} \right]~~,\end{aligned}\end{split}\]

with the force on an atom \(j\) given by

\[\begin{split}\begin{aligned} \underline{f}_{j} = {q_{i} q_{j} \over 4\pi\epsilon_{0}\epsilon} \left[ \left( \frac{{\rm erfc}(\alpha~r_{ij})}{r_{ij}^{2}} + \frac{2\alpha}{\sqrt{\pi}}~\frac{\exp(-\alpha^{2}~r_{ij}^{2})}{r_{ij}} \right) - \phantom{xxxxx} \right. \nonumber \\ \left. \left( \frac{{\rm erfc}(\alpha~r_{\rm cut})}{r_{\rm cut}^{2}} + \frac{2\alpha}{\sqrt{\pi}}~\frac{\exp(-\alpha^{2}~r_{\rm cut}^{2})}{r_{\rm cut}} \right) \right]~\frac{\underline{r}_{ij}}{r_{ij}}~~, \end{aligned}\end{split}\]

with the force on atom \(i\) the negative of this.

It is worth noting that, as discussed in ³¹ and references therein, this is only an approximation of the Ewald sum and its accuracy and effectiveness become better when the cutoff is large (\(>\) 10 preferably 12 Å).

The contribution to the atomic virial is

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

which is not the negative of the potential term in this case.

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

\[\sigma^{\alpha \beta} = r_{ij}^{\alpha} f_{j}^{\beta}~~,\]

where \(\alpha,\beta\) are \(x,y,z\) components. The atomic stress tensor is symmetric.

In DL_POLY_4 these forces are handled by the routine coul_fscp_forces.

Coulomb Sum with Distance Dependent Dielectric

This potential attempts to address the difficulties of applying the direct Coulomb sum, without the brutal truncation of the previous case. It hinges on the assumption that the electrostatic forces are effectively ‘screened’ in real systems - an effect which is approximated by introducing a dielectric term that increases with distance.

The interatomic potential for two charged ions is

\[U(r_{ij}) = \frac{1}{4\pi\epsilon_{0}\epsilon(r_{ij})} \frac{q_{i}q_{j}}{r_{ij}}~~,\]

with \(q_{\ell}\) the charge on an atom labelled \(\ell\), and \(r_{ij}\) the magnitude of the separation vector \(\underline{r}_{ij}=\underline{r}_{j}-\underline{r}_{i}\) . \(\epsilon(r)\) is the distance dependent dielectric function. In DL_POLY_4 it is assumed that this function has the form

\[\epsilon(r)~=~\epsilon~r~~,\]

where \(\epsilon\) is a constant. Inclusion of this term effectively accelerates the rate of convergence of the Coulomb sum.

The force on an atom \(j\) derived from this potential is

\[\underline{f}_{j} = \frac{1}{2\pi\epsilon_{0}\epsilon} \frac{q_{i}q_{j}}{r_{ij}^{4}} \underline{r}_{ij}~~,\]

with the force on atom \(i\) the negative of this.

The contribution to the atomic virial is

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

which is \(-2\) times the potential term.

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

\[\sigma^{\alpha \beta} = r_{ij}^{\alpha} f_{j}^{\beta}~~,\]

where \(\alpha,\beta\) are \(x,y,z\) components. The atomic stress tensor is symmetric.

In DL_POLY_4 these forces are handled by the routine coul_dddp_forces.

Reaction Field

In the reaction field method it is assumed that any given molecule is surrounded by a spherical cavity of finite radius within which the electrostatic interactions are calculated explicitly. Outside the cavity the system is treated as a dielectric continuum. The occurrence of any net dipole within the cavity induces a polarisation in the dielectric, which in turn interacts with the given molecule. The model allows the replacement of the infinite Coulomb sum by a finite sum plus the reaction field.

The reaction field model coded into DL_POLY_4 is the implementation of Neumann based on charge-charge interactions ⁶⁸. In this model, the total coulombic potential is given by

\[U_{c} = \frac{1}{4\pi\epsilon_{0}\epsilon} \sum_{j<n} q_{j}q_{n} \left[ \frac{1}{r_{nj}} + \frac{B_{0}r_{nj}^{2}}{2 R_{c}^{3}} \right]~~,\]

where the second term on the right is the reaction field correction to the explicit sum, with \(R_{c}\) the radius of the cavity. The constant \(B_{0}\) is defined as

\[B_{0} = \frac{2(\epsilon_{1}-1)}{(2\epsilon_{1}+1)}~~,\]

with \(\epsilon_{1}\) the dielectric constant outside the cavity. The effective pair potential is therefore

\[U(r_{ij}) = \frac{1}{4\pi\epsilon_{0}\epsilon} q_{i}q_{j} \left[ \frac{1}{r_{ij}} + \frac{B_{0}r_{ij}^{2}}{2 R_{c}^{3}} \right]~~.\]

This expression unfortunately leads to large fluctuations in the system coulombic energy, due to the large ‘step’ in the function at the cavity boundary. In DL_POLY_4 this is countered by subtracting the value of the potential at the cavity boundary from each pair contribution. The term subtracted is

\[\frac{1}{4\pi\epsilon_{0}\epsilon} \frac{q_{i}q_{j}}{R_{c}} \left[ 1+\frac{B_{0}}{2} \right]~~.\]

The effective pair force on an atom \(j\) arising from another atom \(n\) within the cavity is given by

\[\underline{f}_{j}=\frac{q_{i}q_{j}}{4\pi\epsilon_{0}\epsilon}\left[ \frac{1}{r_{ij}^{3}}-\frac{B_{0}}{R_{c}^{3}}\right] \underline{r}_{ij}~~.\]

In DL_POLY_4 the reaction field is optionally extended to emulate long-range ordering in a force-shifted manner by countering the reaction term and using a distance depending damping function \({\rm erfc}(\alpha~r_{ij})\) (identical to that seen in the real-space portion of the Ewald sum) and thus mirror the effective charge screening ³¹:

\[\begin{split}\begin{aligned} U(r_{ij}) = {q_{i} q_{j} \over 4\pi\epsilon_{0}\epsilon}& \left[ \left\{ \frac{{\rm erfc}(\alpha~r_{ij})}{r_{ij}} + \left(\frac{{\rm erfc}(\alpha~r_{\rm cut})}{r_{\rm cut}^{2}}+ \frac{2\alpha}{\sqrt{\pi}}~\frac{\exp(-\alpha^{2}~r_{\rm cut}^{2})}{r_{\rm cut}}\right)~r_{ij}\right\} \right. \\ &- \left\{ \frac{{\rm erfc}(\alpha~r_{\rm cut})}{r_{\rm cut}} + \left(\frac{{\rm erfc}(\alpha~r_{\rm cut})}{r_{\rm cut}^{2}}+ \frac{2\alpha}{\sqrt{\pi}}~\frac{\exp(-\alpha^{2}~r_{\rm cut}^{2})}{r_{\rm cut}}\right)~r_{\rm cut}\right\} \nonumber \\ &+ \left. \frac{B_{0}(r_{ij}^{2}-r_{\rm cut}^{2})}{2 r_{\rm cut}^{3}} \right]~~, \nonumber \end{aligned}\end{split}\]

with the force on an atom \(j\) given by

\[\begin{split}\begin{aligned} \underline{f}_{j} = {q_{i} q_{j} \over 4\pi\epsilon_{0}\epsilon} \left[ \left( \frac{{\rm erfc}(\alpha~r_{ij})}{r_{ij}^{2}} + \frac{2\alpha}{\sqrt{\pi}}~\frac{\exp(-\alpha^{2}~r_{ij}^{2})}{r_{ij}} \right) - \phantom{xxxxxxxxxxx} \right. \\ \left. \left( \frac{{\rm erfc}(\alpha~r_{\rm cut})}{r_{\rm cut}^{2}} + \frac{2\alpha}{\sqrt{\pi}}~\frac{\exp(-\alpha^{2}~r_{\rm cut}^{2})}{r_{\rm cut}} \right) - \frac{B_{0}r_{ij}}{r_{\rm cut}^{3}} \right]~\frac{\underline{r}_{ij}}{r_{ij}}~~, \nonumber\end{aligned}\end{split}\]

with the force on atom \(i\) the negative of this.

It is worth noting that, as discussed in ³¹ and references therein, this is only an approximation of the Ewald sum and its accuracy and effectiveness become better when the cutoff is large (\(>\) 10 preferably 12 Å).

The contribution of each effective pair interaction to the atomic virial is

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

and the contribution to the atomic stress tensor is

\[\sigma^{\alpha \beta} = r_{ij}^{\alpha} f_{j}^{\beta}~~,\]

where \(\alpha,\beta\) are \(x,y,z\) components. The atomic stress tensor is symmetric.

In DL_POLY_4 the reaction field is handled by the subroutine coul_rfp_forces.

Smoothed Particle Mesh Ewald

The Ewald sum ⁴ is the best technique for calculating electrostatic interactions in a periodic (or pseudo-periodic) system.

The basic model for a neutral periodic system is a system of charged point ions mutually interacting via the Coulomb potential. The Ewald method makes two amendments to this simple model. Firstly, each ion is effectively neutralised (at long-ranged) by the superposition of a spherical Gaussian cloud of opposite charge centred on the ion. The combined assembly of point ions and Gaussian charges becomes the Real Space part of the Ewald sum, which is now short ranged and treatable by the methods described above (Chapter Force Fields) [2]. The second modification is to superimpose a second set of Gaussian charges, this time with the same charges as the original point ions and again centred on the point ions (so nullifying the effect of the first set of Gaussians). The potential due to these Gaussian sum requires an additional correction, known as the self energy correction, which arises from a Gaussian acting on its own site, and is constant. Ewald’s method, therefore, replaces a potentially infinite sum in real space by two finite sums: one in real space and one in reciprocal space; and the self energy correction.

For molecular systems, as opposed to systems comprised simply of point ions, additional modifications are necessary to correct for the excluded (intra-molecular) coulombic interactions. In the real space sum these are simply omitted. In reciprocal space however, the effects of individual Gaussian charges cannot easily be extracted, and the correction is made in real space. It amounts to removing terms corresponding to the potential energy of an ion \(\ell\) due to the Gaussian charge on a neighbouring charge \(m\) (or vice versa). This correction appears in the term noting a summation over \(molecules\) in the full Ewald formula below.

The same considerations and modifications ewald_frzn_forces are taken into account for frozen atoms, which mutual coulombic interaction must also be excluded. This correction appears in the term noting a summation over \(F^{*}\) (all frozen-frozen pairs in the MD cell) in the full Ewald formula below.

Note the distinction between the error function erf and the more usual complementary error function erfc found in the real space sums below.

The total electrostatic energy is given by the following formula:

\[\begin{split}\begin{aligned} U_{c}=&\frac{1}{2V_{o} \epsilon_{0}\epsilon} \sum_{\underline{k} \neq \underline{0}}^{\underline{\infty}} \frac{\exp(-k^{2}/4\alpha^{2})}{k^{2}} \left|\sum_{j}^{N} q_{j}\exp(-i\underline{k}\cdot\underline{r}_{j})\right|^{2} \nonumber \\ &-\frac{1}{4\pi\epsilon_{0}\epsilon} \frac{\alpha}{\sqrt{\pi}} \sum_{j}^{N} q_{j}^2 \nonumber \\ &+ \frac{1}{4\pi\epsilon_{0}\epsilon} \sum_{n<j}^{N^{*}} \frac{q_{j}q_{n}} {r_{nj}} {\rm erfc}(\alpha r_{nj}) \nonumber \\ &- \frac{1}{4\pi\epsilon_{0}\epsilon} \sum_{molecules} \sum_{\ell\le m}^{M^{*}} q_{\ell}q_{m} \left\{ \delta_{\ell m} \frac{\alpha}{\sqrt{\pi}} + \frac{{\rm erf}(\alpha r_{\ell m})}{r_{\ell m}^{1-\delta_{\ell m}}} \right\} \\ &- \frac{1}{4\pi\epsilon_{0}\epsilon} \sum_{\ell\le m}^{F^{*}} q_{\ell}q_{m} \left\{ \delta_{\ell m} \frac{\alpha}{\sqrt{\pi}} + \frac{{\rm erf}(\alpha r_{\ell m})}{r_{\ell m}^{1-\delta_{\ell m}}} \right\} \nonumber \\ &- \frac{1}{4\pi\epsilon_{0}\epsilon} \frac{\pi}{2V_{o} \alpha^{2}} \left\{ \sum_{j}^{N} q_{j} \right\}^{2}~~, \nonumber \end{aligned}\end{split}\]

where \(N\) is the number of ions in the system and \(N^{*}\) the same number discounting any excluded (intramolecular and frozen) interactions. \(M^{*}\) represents the number of excluded atoms in a given molecule. \(F^{*}\) represents the number of frozen atoms in the MD cell. \(V_{o}\) is the simulation cell volume and \(\underline{k}\) is a reciprocal lattice vector defined by

\[\underline{k} = \ell \underline{u} + m \underline{v} + n \underline{w} \label{k-vector}~~,\]

where \(\ell,m,n\) are integers and \(\underline{u},\underline{v},\underline{w}\) are the reciprocal space basis vectors. Both \(V_{o}\) and \(\underline{u},\underline{v},\underline{w}\) are derived from the vectors (\(\underline{a},\underline{b},\underline{c}\)) defining the simulation cell. Thus

\[V_{o} = |\underline{a} \cdot \underline{b} \times \underline{c}|\]

and

\[\begin{split}\begin{aligned} \underline{u} =& 2 \pi \frac{\underline{b} \times \underline{c}}{\underline{a} \cdot \underline{b} \times \underline{c}} \nonumber \\ \underline{v} =& 2 \pi \frac{\underline{c} \times \underline{a}}{\underline{a} \cdot \underline{b} \times\underline{c}} \\ \underline{w} =& 2 \pi \frac{\underline{a} \times \underline{b}}{\underline{a} \cdot \underline{b} \times \underline{c}}~~. \nonumber\end{aligned}\end{split}\]

With these definitions, the Ewald formula above is applicable to general periodic systems. The last term in the Ewald formula above is the Fuchs correction ³⁶ for electrically non-neutral MD cells which prevents the build-up of a charged background and the introduction of extra pressure due to it.

In practice the convergence of the Ewald sum is controlled by three variables: the real space cutoff \(r_{\rm cut}\); the convergence parameter \(\alpha\) and the largest reciprocal space vector \(\underline{k}_{max}\) used in the reciprocal space sum. These are discussed more fully in Section Choosing Ewald Sum Variables. DL_POLY_4 can provide estimates if requested (see CONTROL file description The CONTROL File).

As its name implies the Smoothed Particle Mesh Ewald (SPME) method is a modification of the standard Ewald method. DL_POLY_4 implements the SPME method of Essmann et al. ²⁹. Formally, this method is capable of treating van der Waals forces also, but in DL_POLY_4 it is confined to electrostatic forces only. The main difference from the standard Ewald method is in its treatment of the reciprocal space terms. By means of an interpolation procedure involving (complex) B-splines, the sum in reciprocal space is represented on a three dimensional rectangular grid. In this form the Fast Fourier Transform (FFT) may be used to perform the primary mathematical operation, which is a 3D convolution. The efficiency of these procedures greatly reduces the cost of the reciprocal space sum when the range of \(\underline{k}\) vectors is large. The method (briefly) is as follows (for full details see ²⁹):

1. Interpolation of the \(\exp(-i~\underline{k}\cdot\underline{r}_{j})\) terms (given here for one dimension):

   \[\exp(2\pi i~u_{j} k/L) \approx b(k) \sum_{\ell=-\infty}^{\infty} M_{n}(u_{j}-\ell)~\exp(2\pi i~k\ell/K)~~,\]

   in which \(k\) is the integer index of the \(\underline{k}\) vector in a principal direction, \(K\) is the total number of grid points in the same direction and \(u_{j}\) is the fractional coordinate of ion \(j\) scaled by a factor \(K\) (i.e. \(u_{j}=K s_{j}^{x}\)) . Note that the definition of the B-splines implies a dependence on the integer \(K\), which limits the formally infinite sum over \(\ell\). The coefficients \(M_{n}(u)\) are B-splines of order \(n\) and the factor \(b(k)\) is a constant computable from the formula:

   \[b(k) = \exp(2\pi i~(n-1)k/K) \left[ \sum_{\ell=0}^{n-2} M_{n}(\ell+1)~\exp(2\pi i~k\ell/K) \right]^{-1}~.\]

2. Approximation of the structure factor \(S(\underline{k})\):

   \[S(\underline{k}) \approx b_{1}(k_{1})~b_{2}(k_{2})~b_{3}(k_{3})~Q^{\dagger}(k_{1},k_{2},k_{3})~~,\]

   where \(Q^{\dagger}(k_{1},k_{2},k_{3})\) is the discrete Fourier transform of the charge array \(Q(\ell_{1},\ell_{2},\ell_{3})\) defined as

   \[\begin{split}\begin{aligned} Q(\ell_{1},\ell_{2},\ell_{3})=& \sum_{j=1}^{N}q_{j} \sum_{n_{1},n_{2},n_{3}} M_{n}(u_{1j}-\ell_{1}-n_{1}L_{1})~\times~ M_{n}(u_{2j}-\ell_{2}-n_{2}L_{2}) \phantom{xxxx} \nonumber \\ & \phantom{xxxxxxxxxx}~\times~ M_{n}(u_{3j}-\ell_{3}-n_{3}L_{3})~~,\end{aligned}\end{split}\]

   in which the sums over \(n_{1,2,3}\) etc are required to capture contributions from all relevant periodic cell images (which in practice means the nearest images).

3. Approximating the reciprocal space energy \(U_{recip}\):

   \[U_{recip} = \frac{1}{2V_{o} \epsilon_{0}\epsilon} \sum_{k_{1},k_{2},k_{3}} G^{\dagger}(k_{1},k_{2},k_{3})~Q(k_{1},k_{2},k_{3})~~,\]

   where \(G^{\dagger}\) is the discrete Fourier transform of the function

   \[G(k_{1},k_{2},k_{3}) = \frac{\exp(-k^{2}/4\alpha^{2})}{k^{2}}~ B(k_{1},k_{2},k_{3})~(Q^{\dagger}(k_{1},k_{2},k_{3}))^{*}~~,\]

   in which \((Q^{\dagger}(k_{1},k_{2},k_{3}))^{*}\) is the complex conjugate of \(Q^{\dagger}(k_{1},k_{2},k_{3})\) and

   \[B(k_{1},k_{2},k_{3}) = |b_{1}(k_{1})|^{2}~|b_{2}(k_{2})|^{2}~|b_{3}(k_{3})|^{2}~~.\]

   The function \(G(k_{1},k_{2},k_{3})\) is thus a relatively simple product of the Gaussian screening term appearing in the conventional Ewald sum, the function \(B(k_{1},k_{2},k_{3})\) and the discrete Fourier transform of \(Q(k_{1},k_{2},k_{3})\).

4. Calculating the atomic forces, which are given formally by:

   \[f_{j}^{\alpha} = -\frac{\partial U_{recip}}{\partial r_{j}^{\alpha}} = -\frac{1}{V_{o} \epsilon_{0}\epsilon} \sum_{k_{1},k_{2},k_{3}} G^{\dagger}(k_{1},k_{2},k_{3}) ~\frac{\partial Q(k_{1},k_{2},k_{3})}{\partial r_{j}^{\alpha}}~~.\]

Fortunately, due to the recursive properties of the B-splines, these formulae are easily evaluated.

The virial and the stress tensor are calculated in the same manner as for the conventional Ewald sum.

The DL_POLY_4 subroutines required to calculate the SPME contributions are:

1. spme_container containing

   1. bspgen, which calculates the B-splines

   2. bspcoe, which calculates B-spline coefficients

   3. spl_cexp, which calculates the FFT and B-spline complex exponentials

2. parallel_fft and gpfa_module (native DL_POLY_4 subroutines that respect the domain decomposition concept) which calculate the 3D complex fast Fourier transforms

3. ewald_spme_forces, which calculates the reciprocal space contributions (uncorrected)

4. ewald_real_forces, which calculates the real space contributions (corrected)

5. ewald_excl_forces, which calculates the reciprocal space corrections due to the coulombic exclusions in intramolecular interactions

6. ewald_frzn_forces, which calculates the reciprocal space corrections due to the exclusion interactions between frozen atoms

7. two_body_forces, in which all of the above subroutines are called sequentially and also the Fuchs correction ³⁶ for electrically non-neutral MD cells is applied if needed.

Multipolar Electrostatics

DL_POLY_4 offers advanced potential energy calculations through multipolar electrostatics. This is an extension to the point-charge model where the charge density of chemical species are described by higher order point multipoles. The generic algorithms in DL_POLY_4 are designed to allow for arbitrary order ¹¹ multipoles but for practical reasons the functionality is limited to hexadecapoles only.

Multipoles

Define the multipolar operator, \(\hat{L}_i\) as

\[\hat{L}_i= (q_i + {\mathbf{p}}_i\cdot \nabla_i + {\mathbf{Q}}_i: \nabla_i\nabla_i + {\mathbf{O}}_i{\vdots} \nabla_i\nabla_i\nabla_i + {\mathbf{H}}_i:: \nabla_i\nabla_i\nabla_i\nabla_i + \dots)~~,\]

where \(q_i\), \({\mathbf{p}}_i\), \({\mathbf{Q}}_i\), \({\mathbf{O}}_i\), and \({\mathbf{H}}_i\) are the point charge, dipole, quadrupole, octupole, and hexadecapole tensors, respectively of atom i, \(\nabla_i\) refers to the three-dimensional gradient with respect to the position of atom i and the “dot” products stand for tensor contraction. By defining a unidimensional vector of independent (non-degenerate) multipole moments, \(\mathcal{M}_i\), for atom i, the corresponding multipolar operator to an arbitrary order \(p\) can be written in a more compact form as

(86)\[\hat{L}_i= \sum_{||\mathbf{s}||= 0}^{p}\mathcal{M}_{i}^{\mathbf{s}}\partial_{i}^{\mathbf{s}} = \sum_{s_3 = 0}^{p}\sum_{s_2 = 0}^{p-s_3}\sum_{s_1=0}^{p-s_3-s_2} \mathcal{M}_{j}^{s_1 s_2 s_3} {\partial}_{z_i}^{s_3}{\partial}_{y_i}^{s_2}{\partial}_{x_i}^{s_1}~~.\label{defLi}\]

Here, \(\mathbf{s}= (s_1,s_2,s_3)\) is the triplet that runs over all independent multipoles, \(||\mathbf{s}||= s_1 + s_2 + s_3\), \(\mathcal{M}_{i}^{\mathbf{s}}=\mathcal{M}_{i}^{s_1 s_2 s_3}\) and \(\partial_{i}^{\mathbf{s}} = {\partial}_{z_i}^{s_3}{\partial}_{y_i}^{s_2}{\partial}_{x_i}^{s_1}\) is the multidimensional derivative with respect to the position \(\langle x_i, y_i, z_i \rangle\) of atom i with orders \(s_1\), \(s_2\) and \(s_3\) in the \(x\), \(y\) and \(z\) directions respectively. Individual components of \(\mathcal{M}\) contain the sum of all degenerate original multipole components. As an example, the octupole \(\mathcal{M}^{111}\), is a sum of all six degenerate original octupole components formed from the permutation of the triplet \(\{x,y,z\}\) . If the original octupole vector with degnerate components is labelled as \(O'\), then \(\mathcal{M}^{111}= O_{xyz}' + O_{xzy}' + O_{yxz}' + O_{yzx}' + O_{zxy}' + O_{zyx}' = 6~O_{xyz}'\) . For pair potentials it is often convenient to redefine the multipolar operator for atom j in terms of the derivatives with respect to the position of atom i to arrive at

(87)\[\hat{L}_{j_i}= \sum_{||\mathbf{s}||= 0}^{p}\mathcal{M}_{j}^{\mathbf{s}}\partial_{j}^{\mathbf{s}} = \sum_{||\mathbf{s}||= 0}^{p}(-1)^{||\mathbf{s}||}\mathcal{M}_{j}^{\mathbf{s}}\partial_{i}^{\mathbf{s}} = \sum_{s_3 = 0}^{p}\sum_{s_2 = 0}^{p-s_3}\sum_{s_1=0}^{p-s_3-s_2} (-1)^{s_1+s_2+s_3}\mathcal{M}_{j}^{s_1 s_2 s_3} {\partial}_{z_i}^{s_3}{\partial}_{y_i}^{s_2}{\partial}_{x_i}^{s_1}~~.\label{defLj}\]

Application to Pair Potentials

In DL_POLY_4 for \(N\) point-multipoles interacting via a pair potential function \(\psi\), the multipolar electrostatic potential at position \(\mathbf{r_i}\) is computed as

(88)\[\phi(\mathbf{r_i}) = \sum_{j \ne i}^{N}\hat{L}_{j_i}\psi(\mathbf{r_{ji}}) = \sum_{j \ne i}^{N} \sum_{\mathbf{s}= \mathbf{0}}^{p} (-1)^{||\mathbf{s}||}\mathcal{M}_{j}^{\mathbf{s}}\partial_{i}^{\mathbf{s}}\psi(r_{ij})~~,\label{eqn:mpolpot}\]

the electrostatic field at \(\mathbf{r_i}\) is

(89)\[\begin{split}\mathbf{E}(\mathbf{r_{ij}}) = -\nabla_i \phi(r_{ij})= -\sum_{j \ne i}^{N} \sum_{\mathbf{s}= \mathbf{0}}^{p} (-1)^{||\mathbf{s}||}\mathcal{M}_{j}^{\mathbf{s}}\left[\begin{array}{l} \partial_{i}^{\mathbf{s}+ \mathbf{e}_1} \\ \partial_{i}^{\mathbf{s}+ \mathbf{e}_2} \\ \partial_{i}^{\mathbf{s}+ \mathbf{e}_3} \end{array}\right]\psi(r_{ij})~~,\label{eqn:mpolefield}\end{split}\]

where \(\mathbf{e}_1=\langle1,0, 0\rangle\), \(\mathbf{e}_2=\langle0,1,0\rangle\), and \(\mathbf{e}_3=\langle0,0,1\rangle\) and the torque ⁸⁶ on particle \(i\) in the \(\alpha\)-direction, \(\tau_{i,\alpha}\), is obtained as

\[\tau_{i,\alpha} = \sum_{\mathbf{s}= \mathbf{0}}^{p} \mathcal{M}_{i,\alpha}^{\mathbf{s}} \partial_{i}^{\mathbf{s}} \phi(\mathbf{r_{ij}}) = \sum_{\mathbf{s}= \mathbf{0}}^{p} \mathcal{M}_{i,\alpha}^{\mathbf{s}} \sum_{j \ne i}^{N} \sum_{\mathbf{k}= \mathbf{0}}^{p} (-1)^{||\mathbf{k}||}\mathcal{M}_{j}^{\mathbf{k}}\partial_{i}^{\mathbf{s}+\mathbf{k}}\psi(r_{ij})~~,\]

where \(\mathcal{M}_{i,\alpha}\) is the infinitesimal counter-clockwise rotation of multipole vector \(\mathcal{M}_i\) about the \(\alpha\)-axis. The total electrostatic potential energy is given by

(90)\[U = \sum_{i < j}^{N} \hat{L}_i\hat{L}_{j_i}\psi(r_{ij}) = \sum_{i < j}^{N} \sum_{\mathbf{s}= \mathbf{0}}^{p} (-1)^{||\mathbf{s}||}\mathcal{M}_{j}^{\mathbf{s}}\sum_{\mathbf{k}= \mathbf{0}}^{p}\mathcal{M}_{i}^{\mathbf{k}}\partial_{i}^{\mathbf{s}+\mathbf{k}}\psi(r_{ij})~~,\label{eqn:mpolene}\]

where \(\mathbf{s}+ \mathbf{k}= (s_1+k_1,s_2+k_2,s_3+k_3)\) and the force on atom \(i\) is

(91)\[\begin{split}\mathbf{f}_i = -\nabla_i \sum_{j \ne i}^{N} \hat{L}_i\hat{L}_{j_i}\psi(r_{ij}) = -\sum_{j \ne i}^{N} \sum_{\mathbf{s}= \mathbf{0}}^{p}(-1)^{||\mathbf{s}||}\mathcal{M}_{j}^{\mathbf{s}} \sum_{\mathbf{k}= \mathbf{0}}^{p}\mathcal{M}_{i}^{\mathbf{k}} \left[\begin{array}{l} \partial_{i}^{\mathbf{s}+\mathbf{k}+ \mathbf{e}_1} \\ \partial_{i}^{\mathbf{s}+\mathbf{k}+ \mathbf{e}_2} \\ \partial_{i}^{\mathbf{s}+\mathbf{k}+ \mathbf{e}_3} \end{array}\right] \psi(r_{ij})~~.\label{eqn:mpolforce}\end{split}\]

To implement equations (88)-(91) for the variety of potentials in DL_POLY_4 a number of recurrence relations are used to compute the multi-dimensional derivatives of the kernels corresponding to the potentials. These kernels are

\[\begin{split}\begin{aligned} \theta(|\mathbf{x}|)=& \frac{1}{{|\mathbf{x}|}^{\nu}},\mbox{\hskip 10pt}\Omega(|\mathbf{x}|) \\ =& \frac{1}{2}\textrm{exp}(-\alpha^2|\mathbf{x}|^2), \mbox{\hskip 10pt } \psi(|\mathbf{x}|) \\ =&\frac{\sqrt{\pi}}{2} \frac{\textrm{erfc}(\alpha|\mathbf{x}|)}{|\mathbf{x}|}, \mbox{\hskip 10pt and \hskip 10pt }\Gamma(|\bar{\mathbf{x}}|) \\ =& \frac{\sqrt{\pi}}{2} \frac{\textrm{erf}(\alpha |\mathbf{x}|)}{|\mathbf{x}|}~~; \end{aligned}\end{split}\]

with

\[\begin{split}\begin{aligned} a_{\mathbf{s}}(\nu)&=\frac{\partial^{||\mathbf{s}||}\theta(|\mathbf{x}|)}{\partial_{x_1}^{s_1}\partial_{x_2}^{s_2}\partial_{x_3}^{s_3}},\mbox{\hskip 10pt} b_{\mathbf{s}} \\ &=\frac{\partial^{||\mathbf{s}||}\Omega(|\mathbf{x}|)}{\partial_{x_1}^{s_1}\partial_{x_2}^{s_2}\partial_{x_3}^{s_3}},\mbox{\hskip 10pt} c_{\mathbf{s}} \\ &=\frac{\partial^{||\mathbf{s}||}\psi(|\mathbf{x}|)}{\partial_{x_1}^{s_1}\partial_{x_2}^{s_2}\partial_{x_3}^{s_3}},\mbox{\hskip 10pt and \hskip 10pt} d_{\mathbf{s}} \\ &=\frac{\partial^{||\mathbf{s}||}\Gamma(|\bar{\mathbf{x}}|)}{\partial_{x_1}^{s_1}\partial_{x_2}^{s_2}\partial_{x_3}^{s_3}}~~. \end{aligned}\end{split}\]

The recurrence relations used in DL_POLY_4 are

(92)\[a_{\mathbf{s}}(\nu) = \frac{1}{|\mathbf{x}|^2}\left\{\left(\frac{2-\nu}{||\mathbf{s}||} - 2\right) \sum_{i=1}^{3}s_i x_i a_{\mathbf{s}-\mathbf{e}_i} + \left(\frac{2-\nu}{||\mathbf{s}||} - 1\right) \sum_{i=1}^{3}s_i (s_i-1) a_{\mathbf{s}-2\mathbf{e}_i} \right\}~~,\label{coulrecur}\]

(93)\[b_{\mathbf{s}} = \frac{-2\alpha^2}{||\mathbf{s}||} \sum_{i=1}^{3} \left[ s_i x_i b_{\mathbf{s}-\mathbf{e}_i} + s_i (s_i-1) b_{\mathbf{s}-2\mathbf{e}_i} \right]~~,\label{eqn:exprecur}\]

(94)\[c_{\mathbf{s}} = \frac{1}{|\mathbf{x}|^2} \left\{\left(\frac{1}{||\mathbf{s}||} - 2\right) \sum_{i=1}^{3}s_i x_i c_{\mathbf{s}-\mathbf{e}_i} + \left(\frac{1}{||\mathbf{s}||} - 1\right) \sum_{i=1}^{3}s_i (s_i-1) c_{\mathbf{s}-2\mathbf{e}_i} + \frac{1}{\alpha} b_{\mathbf{s}} \right\}~~,\label{eqn:erfcrecur}\]

and

(95)\[d_{\mathbf{s}} = \frac{1}{|\mathbf{x}|^2}\left\{\left(\frac{1}{||\mathbf{s}||} - 2\right) \sum_{i=1}^{3}s_i x_i d_{\mathbf{s}-\mathbf{e}_i} + \left(\frac{1}{||\mathbf{s}||} - 1\right) \sum_{i=1}^{3}s_i (s_i-1) d_{\mathbf{s}-2\mathbf{e}_i} -\frac{1}{\alpha} b_{\mathbf{s}} \right\}~~.\label{eqn:erfrecur}\]

Direct Coulomb Sum

For two interacting ions \(i\) and \(j\), the potential energy is given as

\[U(r_{ij}) = \frac{1}{4\pi\epsilon_0\epsilon}\hat{L}_i\hat{L}_{j_i}\left[\frac{1}{r_{ij}}\right]~~,\]

and the relevant kernel is \(\psi(r_{ij}) = \frac{1}{r_{ij}}\) . The derivatives for this kernel are obtained by using equation (92) with \(\nu = 1\) . Thus,

\[\partial_i^{\mathbf{s}}\psi(r_{ij}) = a_{\mathbf{s}}(1)~~.\]

In DL_POLY_4 the multipolar direct Coulomb sum is handled by the routine coul_cp_mforces.

Force-Shifted Coulomb Sum

DL_POLY_4 employs two forms of the force-shifted Coulomb sum. In the first form, the potential energy due to two interacting ions \(i\) and \(j\) is

\[U(r_{ij}) = \frac{1}{4\pi\epsilon_0\epsilon} \hat{L}_i\hat{L}_{j_i}\left[\frac{1}{r_{ij}}+\frac{r_{ij}}{r_{\textrm{cut}}^2}-\frac{2}{r_{\textrm{cut}}}\right]~~,\]

where \(r_{\textrm{cut}}\) is the cutoff radius. The kernel is \(\psi(r_{ij}) = \frac{1}{r_{ij}}+\frac{r_{ij}}{r_{\textrm{cut}}^2}-\frac{2}{r_{\textrm{cut}}}\) . The last term, \(\frac{2}{r_{\textrm{cut}}}\), is a constant which has a zero derivative, hence the derivatives of the kernel are obtained as a sum of the derivatives of the first term and second terms. Thus,

\[\partial_i^{\mathbf{s}}\psi(r_{ij}) = a_{\mathbf{s}}(1)+\frac{a_{\mathbf{s}}(-1)}{r_{\textrm{cut}}^2}~~.\]

The potential energy due to two point-multipoles \(i\) and \(j\) interacting via the second form of the force-shifted Coulomb sum is

\[\begin{split}\begin{aligned} U(r_{ij}) =& \frac{1}{4\pi\epsilon_0\epsilon}\hat{L}_i\hat{L}_{j_i}\left[ \left\{ \frac{\textrm{erfc}(\alpha \cdot r_{ij})}{r_{ij}} + \left( \frac{\textrm{erfc}(\alpha \cdot r_{\textrm{cut}})}{r_{\textrm{cut}}^2} + \frac{2\alpha}{\sqrt{\pi}} \frac{\textrm{exp}(-\alpha^2r_{\textrm{cut}}^2)}{r_{\textrm{cut}}}\right)r_{ij}\right\} \right. \nonumber \\ & \left. \phantom{xxxxxxxxx} - \left\{ \frac{\textrm{erfc}(\alpha \cdot r_{\textrm{cut}})}{r_{\textrm{cut}}} + \left( \frac{\textrm{erfc}(\alpha \cdot r_{\textrm{cut}})}{r_{\textrm{cut}}^2} + \frac{2\alpha}{\sqrt{\pi}}\frac{\textrm{exp}(-\alpha^2r_{\textrm{cut}}^2)}{r_{\textrm{cut}}}\right)r_{\textrm{cut}}\right\} \right]~~. \end{aligned}\end{split}\]

The kernel, \(\psi(r_{ij})\) is the terms in the square bracket but the only terms which contribute to the derivatives are the first and second terms which are functions of \(r_{ij}\) . The derivative of the first term is obtained from equations (94) and the derivative for \(r_{ij}\) in the second term is given by \(d_{\mathbf{s}}(-1)\) . Thus,

\[D_i^{\mathbf{s}}\psi(r_{ij}) = \frac{2}{\sqrt{\pi}}c_{\mathbf{s}} + \left(\frac{\textrm{erfc}(\alpha \cdot r_{\textrm{cut}})}{r_{\textrm{cut}}^2} + \frac{2\alpha}{\sqrt{\pi}} \frac{\textrm{exp}(-\alpha^2r_{\textrm{cut}}^2)}{r_{\textrm{cut}}} \right) \cdot a_{\mathbf{s}}(-1)~~.\]

In DL_POLY_4 the multipolar force-shifted Coulomb sum is handled by the routine coul_fscp_mforces.

Coulomb Sum with Distance Dependent Dielectric

The potential energy between two interacting ions \(i\) and \(j\) is

\[U(r_{ij}) = \frac{1}{4\pi\epsilon_0\epsilon}\hat{L}_i\hat{L}_{j_i}\left[\frac{1}{r_{ij}^2}\right]~~,\]

and the kernel is \(\psi(r_{ij}) = \frac{1}{r_{ij}^2}\) . The derivatives for this kernel are obtained by using equation (92) with \(\nu = 2\) . Hence,

\[\partial_i^{\mathbf{s}}\psi(r_{ij}) = a_{\mathbf{s}}(2)~~.\]

In DL_POLY_4 the multipolar Coulomb sum with distance dependent dielectric is handled by the routine coul_dddp_mforces.

Reaction Field

DL_POLY_4 provides two forms of a multipolar reaction field potential. In the first form, the effective pair potential energy due to two interacting point multipoles \(i\) and \(j\) is given as

\[U(r_{ij}) = \frac{1}{4\pi\epsilon_0\epsilon}\hat{L}_i\hat{L}_{j_i}\left[ \frac{1}{r_{ij}} + \frac{B_0r_{ij}^2}{2R_c^3} - 1 - \frac{B_0}{2}\right]~~,\]

where

\[B_0 = \frac{2(\epsilon_1 - 1)}{(2\epsilon_1 + 1)}~~,\]

\(R_c\) is the radius of the spherical cavity and \(\epsilon_1\) is the dielectric constant outside the cavity. Again the kernel \(\psi(r_{ij})\) is the terms in the square bracket and only the first and second terms contribute to its derivatives. The derivatives of the first and second terms are given by equation (92) with \(\nu = 1\) and \(\nu = -2\) respectively. Thus,

\[\partial_i^{\mathbf{s}}\psi(r_{ij}) = a_{\mathbf{s}}(1)+\frac{B_0}{2R_c^3} \cdot a_{\mathbf{s}}(-2)~~.\]

The second form of the reaction field method is similar to that of the force-shifted Coulomb sum. The potential energy due to interacting ions \(i\) and \(j\) is

\[\begin{split}\begin{aligned} U(r_{ij}) =& \frac{1}{4\pi\epsilon_0\epsilon}\hat{L}_i\hat{L}_{j_i}\left[ \left\{ \frac{\textrm{erfc}(\alpha \cdot r_{ij})}{r_{ij}} + \left( \frac{\textrm{erfc}(\alpha \cdot r_{\textrm{cut}})}{r_{\textrm{cut}}^2} + \frac{2\alpha}{\sqrt{\pi}} \frac{\textrm{exp}(-\alpha^2r_{\textrm{cut}}^2)}{r_{\textrm{cut}}}\right)r_{ij}\right\} - \phantom{xxxxxxxxx} \right. \\ \nonumber & \left. \phantom{xxxxxxxxx} \left\{ \frac{\textrm{erfc}(\alpha \cdot r_{\textrm{cut}})}{r_{\textrm{cut}}} + \left( \frac{\textrm{erfc}(\alpha \cdot r_{\textrm{cut}})}{r_{\textrm{cut}}^2} + \frac{2\alpha}{\sqrt{\pi}}\frac{\textrm{exp}(-\alpha^2r_{\textrm{cut}}^2)}{r_{\textrm{cut}}}\right)r_{\textrm{cut}}\right\} - \frac{B_0 r_{\textrm{cut}}^2}{2r_{\textrm{cut}}^3}+\frac{B_0r_{ij}^2}{2r_{\textrm{cut}}^3} \right]~~. \end{aligned}\end{split}\]

The kernel, \(\psi(r_{ij})\) is the terms in the square bracket and the only terms which contribute to the derivatives are the first, second and last terms which are functions of \(r_{ij}\) . The derivative of the first term is obtained from equation (94) and the derivative for \(r_{ij}\) in the second term is given by \(a_{\mathbf{s}}(-1)\) and the derivative for \(r_{ij}^2\) in the last term is given by \(d_{\mathbf{s}}(-2)\) . Thus,

\[D_i^{\mathbf{s}}\psi(r_{ij}) = \frac{2}{\sqrt{\pi}}c_{\mathbf{s}}+\left( \frac{\textrm{erfc}(\alpha \cdot r_{\textrm{cut}})}{r_{\textrm{cut}}^2} + \frac{2\alpha}{\sqrt{\pi}}\frac{\textrm{exp}(-\alpha^2r_{\textrm{cut}}^2)}{r_{\textrm{cut}}}\right) \cdot a_{\mathbf{s}}(-1) + \frac{B_0}{2r_{\textrm{cut}}^3} \cdot a_{\mathbf{s}}(-2)~~.\]

In DL_POLY_4 the multipolar reaction field is handled by the routine coul_rfp_mforces.

Smoothed Particle Mesh Ewald

DL_POLY_4 provides two different smooth particle Mesh Ewald implementations for multipolar electrostatics. The first implementation is for systems with charges, dipoles and quadrupoles and does not use recurrence relations. The second implementation, which uses recurrence relations, is more general and allows for specification of an arbitrary order up to hexadecapoles.

When the multipolar form of SPME is employed, the total electrostatic energy for a system on \(N\) point ions is given as

(96)\[U_c = U_{\textrm{dir}} + U_{\textrm{rec}} - U_{\textrm{excl}} - U_{\textrm{frzn}} - U_{\textrm{self}}~~,\label{eqn:totE}\]

where

(97)\[U_{\textrm{dir}} = \sum_{i < j}^{N^*} \sum_{\mathbf{n}}^{'}\hat{L}_i\hat{L}_{j_i} \frac{\textrm{erfc}(\alpha \cdot |\mathbf{r_{ij}}+ \mathbf{n}|)}{4\pi\epsilon_0\epsilon|\mathbf{r_{ij}}+ \mathbf{n}|}~~,\label{direwald}\]

(98)\[U_{\textrm{excl}} = \frac{1}{4\pi\epsilon_0\epsilon}\sum_{(i,j)\in M^*} \hat{L}_i\hat{L}_{j_i}\frac{\textrm{erf}(\alpha \cdot r_{ij})}{r_{ij}}~~,\label{exclewald}\]

(99)\[U_{\textrm{frzn}} = \frac{1}{4\pi\epsilon_0\epsilon}\sum_{(i,j)\in F^*} \hat{L}_i\hat{L}_{j_i}\frac{\textrm{erf}(\alpha \cdot r_{ij})}{r_{ij}}~~,\label{frznewald}\]

(100)\[U_{\textrm{self}} = \frac{1}{8\pi\epsilon_0\epsilon}\lim_{| \mathbf{r_i}|\to 0}\sum_{i = 1}^{N} \hat{L}_i\hat{L}_i\frac{\textrm{erf}(\alpha \cdot |\mathbf{r_i}|)}{|\mathbf{r_i}|}~~,\label{selfewald}\]

and

(101)\[U_{\textrm{rec}} = \frac{1}{2V_o\epsilon_0\epsilon}\displaystyle\sum_{\mathbf{k} \ne 0} \frac{\textrm{exp}(-k^2/4\alpha^2)}{k^2}\left |S(\mathbf{k})\right|^2~~,\label{recewald}\]

with

(102)\[S(\mathbf{k}) = \sum_{i=1}^N \hat{L}_i\textrm{exp}(\imath \mathbf{k}\cdot \mathbf{r_i})~~.\label{eqn:sfac}\]

In the expressions above, \(M^*\) is the set of all excluded interactions due to intramolecular bonds in the simulation cell, \(F^*\) the set of frozen-frozen interactions in the simulation cell, \(N^* = N - M^* - F^*\), \(V_o\) is the volume of the simulation cell and \(S(\mathbf{k})\) is the structure factor.

Real Space Sum

The relevant kernel for the real space from equation (97) is \(\displaystyle \psi(r_{ij}) = \frac{\textrm{erfc}(\alpha |\mathbf{r_{ij}}+ \mathbf{n}|)}{|\mathbf{r_{ij}}+ \mathbf{n}|}\) . DL_POLY_4 uses the recurrence giving in equation (94) to generate the multidimensional derivatives of the kernel. Thus, the derivatives of the kernel are computed as

\[\mathbf{D}_i^{\mathbf{s}}\psi(r_{ij}) = \frac{2}{\sqrt{\pi}}c_{\mathbf{s}}~~.\]

In DL_POLY_4 the routine ewald_real_mforces_d computes the real space interactions explicitly for simulations with multipoles of order 2 without using the recurrence relation. The routine ewald_real_mforces handles the general version of up to order 4 using recurrence relations.

Excluded Sum

The relevant kernel for the real space from equation (98) is \(\displaystyle \psi(r_{ij}) = \frac{\textrm{erf}(\alpha \cdot r_{ij})}{r_{ij}}\) . DL_POLY_4 uses the recurrence giving in equation (95) to generate the multidimensional derivatives of the kernel. Thus, the derivatives of the kernel are computed as

\[\mathbf{D}_i^{\mathbf{s}}\psi(r_{ij}) = \frac{2}{\sqrt{\pi}}d_{\mathbf{s}}~~.\]

In DL_POLY_4 the routine ewald_excl_mforces_d computes the reciprocal space corrections due to the exclusions between intramolecularly related atoms explicitly for simulations with multipoles of order 2 without using the recurrence relation. The routine ewald_excl_mforces handles the general version of up to order 4 using recurrence relations.

Frozen Sum

The relevant kernel for the real space from equation (99) is \(\displaystyle \psi(r_{ij}) = \frac{\textrm{erf}(\alpha \cdot r_{ij})}{r_{ij}}\) . DL_POLY_4 uses the recurrence giving in equation (95) to generate the multidimensional derivatives of the kernel. Thus, the derivatives of the kernel are computed as

\[\mathbf{D}_i^{\mathbf{s}}\psi(r_{ij}) = \frac{2}{\sqrt{\pi}}d_{\mathbf{s}}~~.\]

In DL_POLY_4 the routine ewald_frzn_mforces computes computes the reciprocal space corrections due to the exclusions between frozen atoms generically for simulations with multipoles up to order 4 using recurrence relations.

Self-Interaction

DL_POLY_4 computes \(U_{self}\) directly for interactions involving multipoles up to order 4 using the series representation of the kernel \(\displaystyle \psi(r_{ij}) = \frac{\textrm{erf}(\alpha \cdot r_i)}{r_i}\) . The self interaction is computed in ewald_real_mforces_d for simulations with multipoles of maximum order 2. For simulations of arbitrary order, the self-interaction is computed in the reciprocal space.

Reciprocal Space Sum

The key idea of SPME is in approximating the structure factor, in a uniform grid, with \(K_1 \times K_2 \times K_3\) dimensions, that fills the simulation cell. Define the fractional coordinates of an ion \(i\) as \(\langle s_{i_{1}}, s_{i_{2}}, s_{i_{3}} \rangle = \langle \mathbf{a}_1^*\cdot \mathbf{r_i}, \mathbf{a}_2^*\cdot \mathbf{r_i}, \mathbf{a}_3^*\cdot \mathbf{r_i}\rangle\), \(u_{\alpha_i} = K_{\alpha} \cdot s_i^{\alpha}\) and \(M_n\) is a B-spline of order \(n\) then the approximation of the structure factor is given as

(103)\[S(\mathbf{k}) \approx b_1(k_1)b_2(k_2)b_3(k_3) Q^{\mathcal{F}}(k_1,k_2,k_3)~~,\label{skapprox}\]

where \(\mathbf{k}= \langle k_1, k_2, k_3 \rangle\) is a reciprocal space vector,

\[b_i(k_i) = \textrm{exp}(2\pi \imath (n-1)k_i/K_i)\left[\sum_{l=0}^{n-2}M_n(l+1)\textrm{exp}(2\pi \imath k l/K_i)\right ]^{-1},\]

\(Q\) is the multipolar array defined on the uniform grid and \(Q^{\mathcal{F}}\) its discrete Fourier transform. At position \((l_1,l_2,l_3)\) on the grid, the multipolar array is defined by

(104)\[\begin{split}\begin{aligned} Q(l_1,l_2,l_3)=\sum_{i=1}^{N}\hat{L}_i\sum_{n_1,n_2,n_3} &M_n(u_{1_i}-l_1-n_1 K_1) \times M_n(u_{2_i}-l_2-n_2 K_2) \\ &\times M_n(u_{3_i}-l_3-n_3 K_3)~~, \end{aligned} \label{marray1}\end{split}\]

where, \(u_{\alpha_i}-l_{\alpha}-n_{\alpha}K_{\alpha}\) are evaluation points of the B-spline on the grid that spans the fundamental cell and the periodic images. Then from equation (86) and considering only the fundamental cell, the multipolar array can be written explicitly as

(105)\[\begin{aligned} Q(l_1,l_2,l_3) = \sum_{i=1}^{N}\sum_{s_3 = 0}^{p}\sum_{s_2 = 0}^{p-s_3} \sum_{s_1=0}^{p-s_3-s_2}\mathcal{M}_{i}^{s_1 s_2 s_3} {\partial}_{z_i}^{s_3} {\partial}_{y_i}^{s_2}{\partial}_{x_i}^{s_1} \left\{ M_n(u_{1_i}-l_1) M_n(u_{2_i}-l_2) M_n(u_{3_i}-l_3)\right\}~~. \end{aligned} \label{marray2}\]

To compute the arbitrary order multidimensional derivatives of the product of three b-splines in equation (105), DL_POLY_4 uses the closed form formula:

(106)\[\begin{split}\begin{aligned} & {\partial}_{z_i}^{s_3}{\partial}_{y_i}^{s_2}{\partial}_{x_i}^{s_1} \left\{ M_n(u_{1_i}-l_1) M_n(u_{2_i}-l_2) M_n(u_{3_i}-l_3)\right\} = \nonumber \\ & \displaystyle \sum_{k_3 = 0}^{s_3}\left(K_1 a_{13}^{*}\right)^{k_3} \binom{s_3}{k_3} \sum_{k_2=0}^{s_2}\left(K_1 a_{12}^{*}\right)^{k_2} \binom{s_2}{k_2}\sum_{k_1=0}^{s_1}\left(K_1 a_{11}^{*}\right)^{k_1} \binom{s_1}{k_1} {\partial}_{u_{1_i}}^{||\mathbf{k}||}M_n(u_{1_i}-l_1) \times \\ & \displaystyle \sum_{j_3=0}^{s_3-k_3}\left(K_2 a_{23}^{*}\right)^{j_3} \left(K_3 a_{33}^{*}\right)^{s_3-k_3-j_3}\binom{s_3-k_3}{j_3} \sum_{j_2=0}^{s_2-k_2}\left(K_2 a_{22}^{*}\right)^{j_2} \left(K_3 a_{32}^{*}\right)^{s_2-k_2-j_2} \binom{s_2-k_2}{j_2} \times \nonumber \\ & \displaystyle \sum_{j_1=0}^{s_1-k_1}\left(K_2 a_{21}^{*}\right)^{j_1} \left(K_3 a_{31}^{*}\right)^{s_1-k_1-j_1} \binom{s_1-k_1}{j_1} {\partial}_{u_{2_i}}^{||\mathbf{j}||}M_n(u_{2_i}-l_2){\partial}_{u_{3_i}}^{||\mathbf{s}-\mathbf{k}-\mathbf{j}||}M_n(u_{3_i}-l_3)~~,\nonumber \label{dprodmn} \end{aligned}\end{split}\]

where \(\mathbf{a}_{1}^{*} = \langle a_{11}^{*},a_{12}^{*},a_{13}^{*}\rangle\), \(\mathbf{a}_{2}^{*} = \langle a_{21}^{*},a_{22}^{*},a_{23}^{*}\rangle\), and \(\mathbf{a}_{3}^{*} = \langle a_{31}^{*},a_{32}^{*},a_{33}^{*}\rangle\) are the reciprocal space basis vectors and \(K_1\), \(K_2\), and \(K_3\), the maximum number of grid points in the fundamental cell in the \(x\), \(y\), and \(z\) directions respectively. For an orthogonal box, where

\[a_{12}^{*}=a_{13}^{*}=a_{21}^{*}=a_{23}^{*}=a_{31}^{*}=a_{32}^{*}=0~~,\]

DL_POLY_4 uses the simplification of equationv (106) to

(107)\[\begin{split}\begin{aligned} \label{dprodmnsimple} & {\partial}_{z_i}^{s_3}{\partial}_{y_i}^{s_2}{\partial}_{x_i}^{s_1} \left\{ M_n(u_{1_i}-l_1) M_n(u_{2_i}-l_2) M_n(u_{3_i}-l_3)\right\}= \\ & \left(K_1 a_{11}^{*}\right)^{s_1}\left(K_2 a_{22}^{*}\right)^{s_2} \left(K_3 a_{33}^{*}\right)^{s_3}{\partial}_{u_{1_i}}^{s_1}M_n(u_{1_i}-l_1) {\partial}_{u_{2_i}}^{s_2}M_n(u_{2_i}-l_2){\partial}_{u_{3_i}}^{s_3}M_n(u_{3_i}-l_3)~~. \nonumber \end{aligned}\end{split}\]

The formulas in equations (106) and (107) require derivatives of a b-spline. To compute an arbitrary \(p_\textrm{th}\) order derivative of a b-spline of order \(n\), \(M_n\), at an arbitrary grid point \(j\), DL_POLY_4 uses the closed form formula

(108)\[\frac{d^p}{d u^{p}}M_n(u_j) = \sum_{t=\textrm{max}\{0,j-k\}}^{\textrm{min}\{j-1,p\}} \binom{p}{t}(-1)^{t}M_{k}(u_j-t)~~.\label{dmnj2}\]

In DL_POLY_4 the stress tensor due to the reciprocal space, for an arbitrary \(p_\textrm{th}\) order multipolar electrostatic interaction is computed by the formula

(109)\[\begin{split}\begin{aligned} V\sigma_{\alpha \beta}^{\textrm{rec}} = \frac{1}{2V_o\epsilon_0\epsilon} \displaystyle \sum_{\mathbf{k}\ne 0} \frac{\textrm{exp}(-k^2/4\eta^2)}{k^2} &\left\{ \left |S(\mathbf{k})\right|^2 \left[\delta_{\alpha \beta}-2 \left( \frac{k^2/4\eta^2 + 1}{k^2}\right) k_{\alpha} k_{\beta}\right] \right. \nonumber \\ & \left. +2S(\mathbf{k}) S_i^{\beta}(-\mathbf{k}) \frac{k_{\alpha}}{k_{\beta}} \right\}~~, \end{aligned}\end{split}\]

where

\[\mathcal{J}_i^{\textbf{ \ell}}(\mathbf{k}) = \mathcal{M}_i^{\textbf{ \ell}}{\partial}_i^{\textbf{ \ell}}{\textrm{e}}^{\imath \mathbf{k}\cdot \mathbf{r_i}}~~,\]

\[S_i^{\beta}(-\mathbf{k}) = \sum_{\textbf{ \ell}= \mathbf{0}}^{p} \ell_{\beta} \sum_{i=1}^{N} \mathcal{J}_i^{\textbf{ \ell}}(-\mathbf{k})~~,\]

and \(\textbf{ \ell}= (\ell_1,\ell_2,\ell_3)\) .

In DL_POLY_4 the routine ewald_spme_mforces_d computes the reciprocal space interactions explicitly for simulations with multipoles of maximum order 2. The routine ewald_spme_mforces handles the general version with multipoles up to order 4.

The DL_POLY_4 subroutines required to calculate the contributions from the reciprocal space, in addition to the routines used for the point charges, are:

1. bspgen_mpoles, in spme_container evaluates equation (106) or (107) to compute the B-splines.

2. limit_erfr_deriv in mpoles_container which computes the limit of the derivatives of the kernel for the self-interaction term. limit_erfr_deriv is called in ewald_spme_mforces.

[1]

Unlike the other elements of the force field, the electrostatic forces are NOT specified in the input FIELD file, but by setting appropriate directives in the CONTROL file. See Section The CONTROL File.

[2]

Strictly speaking, the real space sum ranges over all periodic images of the simulation cell, but in the DL_POLY_4 implementation, the parameters are chosen to restrict the sum to the simulation cell and its nearest neighbours, i.e. the minimum images of the cell contents.