Hautman Klein Ewald (HKE)

The method of Hautman and Klein is an adaptation of the Ewald method for systems which are periodic in two dimensions only [32]. (DL_POLY_2 assumes this periodicity is in the XY plane.)

The HKE method gives the following formula for the electrostatic energy of a system of $N$ (nonbonded) ions that is overall charge neutral^2.4:

$$U_{c} = \frac{1}{4\epsilon_{0}A}\sum_{n=0}^{n_{max}}a_{n}\sum_{i,j}^{N} q... ...erline{0}$}} f_{n}(g;\alpha) g^{2n-1}exp(i\mbox{$\underline{g}$}\cdot s_{ij}) +$$

$$\frac{1}{8\pi\epsilon_{0}}\sum_{i\ne j}^{N}q_{i}q_{j}\sum_{\mbox{... ...{max}} a_{n}z_{ij}^{2n} \frac{h_{n}(s_{ij,L};\alpha)}{s_{ij,L}^{2n+1}}\right )+$$

$$\frac{1}{8\pi\epsilon_{0}}\sum_{i}^{N}q_{i}^{2}\sum_{\mbox{$\unde... ...-h_{0}(L;\alpha))}{L}-\frac{\alpha}{\epsilon_{0}\pi^{3/2}}\sum_{i}^{N}q_{i}^{2}$$ (2.147)

In this formula $A$ is the system area (in the XY plane), $\mbox{$\underline{L}$}$ is a 2D lattice vector representing the 2D periodicity of the system, $s_{ij}$ is the in-plane (XY) component of the interparticle distance $r_{ij}$ and $\mbox{$\underline{g}$}$ is a reciprocal lattice vector. Thus

$$\mbox{$\underline{L}$}=\ell_{1}\mbox{$\underline{a}$}+\ell_{2}\mbox{$\underline{b}$},$$ (2.148)

where $\ell_{1},\ell_{2}$ are integers and vectors $\mbox{$\underline{a}$}$ and $\mbox{$\underline{b}$}$ are the lattice basis vectors. The reciprocal lattice vectors are:

$$\mbox{$\underline{g}$}=n_{1} \mbox{$\underline{u}$} + n_{2} \mbox{$\underline{v}$}$$ (2.149)

where $n_{1},n_{2}$ are integers $\mbox{$\underline{u}$},\mbox{$\underline{v}$}$ are reciprocal space vectors (defined in terms of the vectors $\mbox{$\underline{a}$}$ and $\mbox{$\underline{b}$}$):

$$\mbox{$\underline{u}$} = 2\pi(b_{y},-b_{x})^{\dagger}/(a_{x}b_{y}-a_{y}b_{x})$$

$$\mbox{$\underline{v}$} = 2\pi(-a_{y},a_{x})^{\dagger}/(a_{x}b_{y}-a_{y}b_{x}).$$ (2.150)

The functions $h_{n}(s;\alpha)$ and $f_{n}(s;\alpha)$ are the HKE convergence functions, in real and reciprocal space respectively. (C.f. the complementary error and gaussian functions of the original Ewald method.) However they occur to higher orders here, as indicated by the sum over subscript $n$, which corresponds to terms in a Taylor expansion of $r^{-1}$ in $s$, the in-plane distance [32]. Usually this sum is truncated at $n_{max}=1$, but in DL_POLY_2 can go as high as $n_{max}=3$. In the HKE method the convergence functions are defined as follows:

$$h_{n}(s;\alpha)/s^{2n+1}=\frac{1}{a_{n}(2n)!}\nabla^{2n}(h_{0}(s;\alpha)/s)$$ (2.151)

with

$$h_{0}(s;\alpha)=erf(\alpha s)$$ (2.152)

and

$$f_{n}(g;\alpha)=\frac{1}{a_{n}(2n)!}f_{0}(g;\alpha)$$ (2.153)

with

$$f_{0}(g;\alpha)=erfc(g/2\alpha).$$ (2.154)

In DL_POLY_2 the $h_{n}(s;\alpha)/s^{2n+1}$ functions are derived by a recursion algorithm, while the $f_{n}(g;\alpha)$ functions are obtained by direct evaluation. The coefficients $a_{n}$ are given by

$$a_{n}=(-1)^{n}(2n)!/(2^{2n}(n!)^{2}).$$ (2.155)

As pointed out by Hautman and Klein, the equation (2.147) allows separation of the $z_{ij}^{2n}$ components via the binomial expansion, which greatly simplifies the double sum over atoms in reciprocal space. Thus the reciprocal space part of equation (2.147) becomes

$$U_{recip}=\frac{1}{4\epsilon_{0}A}\sum_{n=0}^{n_{max}} a_{n}... ...p}(\mbox{$\underline{g}$})Z_{2n-p}^{*}(\mbox{$\underline{g}$})$$ (2.156)

with $C_{p}^{2n}$ a binomial coefficient and

$$Z_{p}(\mbox{$\underline{g}$})=\sum_{j=1}^{N}q_{j}z_{j}^{p}exp(i\mbox{$\underline{g}$}\cdot \mbox{$\underline{s_{j}}$})$$ (2.157)

The force on an ion is obtained by the usual differentiation, however in this case the z components have different expressions from the x and y.

$$-\frac{\partial U_{c}}{\partial u_{j}} = \frac{1}{4\epsilon_{0}A}\sum_{\mbox{$\underline{g}$}\ne \mbox{$\u... ...ine{g}$})\frac{\partial Z_{p}(\mbox{$\underline{g}$})}{\partial u_{j}} \right )$$

$$+\frac{q_{j}}{4\pi \epsilon_{0}}\sum_{n=0}^{n_{max}}\sum_{\mbox{$... ...}}\left ( z_{ij,L}^{2n} \frac{h_{n}(s_{ij,L};\alpha)}{s_{ij,L}^{2n+1}} \right )$$ (2.158)

where $u_{j}$ is one of $x_{j},y_{j},z_{j}$ and (noting for brevity that $x$ and $y$ derivatives are similar)

$$\frac{\partial Z_{p}(\mbox{$\underline{g}$})}{\partial x_{j}} = ig_{x}q_{j}z_{j}^{p} exp(i\mbox{$\underline{g}$}\cdot\mbox{$\underline{s_{j}}$})$$

$$\frac{\partial Z_{p}(\mbox{$\underline{g}$})}{\partial z_{j}} = pq_{j}z_{j}^{p-1} exp(i\mbox{$\underline{g}$}\cdot\mbox{$\underline{s_{j}}$})$$ (2.159)

and

$$\frac{\partial}{\partial x_{j}}\left ( z_{ij,L}^{2n} \frac{h_{n}(s_{ij,L};\alpha)}{s_{ij,L}^{2n+1}} \right ) = s_{ij,L}^{x}\frac{z_{ij,L}^{2n} }{s_{ij,L}} \frac{\partial}{\partial x_{j}}\left (\frac{h_{n}(s_{ij,L};\alpha)}{s_{ij,L}^{2n+1}} \right )$$

$$\frac{\partial}{\partial z_{j}}\left ( z_{ij,L}^{2n} \frac{h_{n}(s_{ij,L};\alpha)}{s_{ij,L}^{2n+1}} \right ) = 2n z_{ij,L}^{2n-1} \frac{h_{n}(s_{ij,L};\alpha)}{s_{ij,L}^{2n+1}}.$$ (2.160)

In DL_POLY_2 the partial derivatives of $h_{n}(s_{ij,L};\alpha)/s_{ij,L}^{2n+1}$ are calculated by a recursion algorithm. Note that when $n=0$ there is no derivative w.r.t. $z$.

The virial and stress tensor terms in real space may be calculated directly from the pair forces and interatomic distances in the usual way, and need not be discussed further. The calculation of the reciprocal space contributions (the terms involving the $f_{n}(g;\alpha)$ functions) are more difficult. Firstly however we note that the reciprocal space contributions to $\sigma_{xz},\sigma_{yz}$ and $\sigma_{zz}$ may be obtained directly from the force calculations thus:

$$\sigma^{recip}_{xz} = \sum_{j} z_{j}f_{j}^{x}$$

$$\sigma^{recip}_{yz} = \sum_{j} z_{j}f_{j}^{y}$$ (2.161)

$$\sigma^{recip}_{zz} = \sum_{j} z_{j}f_{j}^{z}$$

which renders the calculation of these components trivial. The remaining components are calculated from

$$\sigma^{recip}_{uv} = U_{recip}\delta_{uv}+\frac{1}{4\epsilon_{0}A} \sum_{n=0}^{n_{max}... ...underline{g}$}\ne \mbox{$\underline{0}$}} g_{u}g_{v}\frac{g^{2n-2}}{a_{n}(2n)!}$$

$$\left (\frac{(2n-1)f_{0}(g;\alpha)}{g}-\frac{1}{\alpha\surd{\pi}}exp(-g^{2}/4\alpha^{2})\right )$$ (2.162)

$$\sum_{p=0}^{2n}(-1)^{p}C_{p}^{2n}Z_{p}(\mbox{$\underline{g}$})Z_{2n-p}^{*}(\mbox{$\underline{g}$})$$

where $u,v$ are one or both of the components $x,y$. Note that, although it is possible to define these contributions to the stress tensor, it is not possible to calculate a pressure from them unless a finite, arbitrary boundary is imposed on the z direction (which is an assumption applied in DL_POLY_2 , but without implications of periodicity in the z-direction). The $x,y$ components define the surface tension however.

For bonded molecules, as with the standard 3D Ewald sum, it is necessary to extract contributions associated with the excluded atom pairs. In the DL_POLY_2 HKE implementation this amounts to an a posteriori subtraction of the corresponding coulomb terms.

In DL_POLY_2 the HKE method is handled by several subroutines: HKGEN constructs the $h_{n}(s;\alpha)$ convergence functions and their derivatives; HKEWALD1 calculates the reciprocal space terms; HKEWALD2 and HKEWALD3 calculate the real space terms and the bonded atom corrections respectively. HKEWALD4 calculates the primary interactions in the multiple timestep implementation.