Ewald summation

Ewald summation, named after Paul Peter Ewald, is a method for computing long-range interactions (e.g. electrostatic interactions) in periodic systems. It was first developed as the method for calculating the electrostatic energies of ionic crystals, and is now commonly used for calculating long-range interactions in computational chemistry. Ewald summation is a special case of the Poisson summation formula, replacing the summation of interaction energies in real space with an equivalent summation in Fourier space. In this method, the long-range interaction is divided into two parts: a short-range contribution, and a long-range contribution which does not have a singularity. The short-range contribution is calculated in real space, whereas the long-range contribution is calculated using a Fourier transform. The advantage of this method is the rapid convergence of the energy compared with that of a direct summation. This means that the method has high accuracy and reasonable speed when computing long-range interactions, and it is thus the de facto standard method for calculating long-range interactions in periodic systems. The method requires charge neutrality of the molecular system to accurately calculate the total Coulombic interaction. A study of the truncation errors introduced in the energy and force calculations of disordered point-charge systems is provided by Kolafa and Perram.

Ewald summation rewrites the interaction potential as the sum of two terms, $${\displaystyle \varphi (\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \varphi _{sr}(\mathbf {r} )+\varphi _{\ell r}(\mathbf {r} ),}$$ where ${\displaystyle \varphi _{sr}(\mathbf {r} )}$ represents the short-range term whose sum quickly converges in real space and ${\displaystyle \varphi _{\ell r}(\mathbf {r} )}$ represents the long-range term whose sum quickly converges in Fourier (reciprocal) space. The long-ranged part should be finite for all arguments (most notably r = 0) but may have any convenient mathematical form, most typically a Gaussian distribution. The method assumes that the short-range part can be summed easily; hence, the problem becomes the summation of the long-range term. Due to the use of the Fourier sum, the method implicitly assumes that the system under study is infinitely periodic (a sensible assumption for the interiors of crystals). One repeating unit of this hypothetical periodic system is called a unit cell. One such cell is chosen as the "central cell" for reference and the remaining cells are called images.

The long-range interaction energy is the sum of interaction energies between the charges of a central unit cell and all the charges of the lattice. Hence, it can be represented as a double integral over two charge density fields representing the fields of the unit cell and the crystal lattice $${\displaystyle E_{\ell r}=\iint d\mathbf {r} \,d\mathbf {r} ^{\prime }\,\rho _{\text{TOT}}(\mathbf {r} )\rho _{uc}(\mathbf {r} ^{\prime })\ \varphi _{\ell r}(\mathbf {r} -\mathbf {r} ^{\prime })}$$ where the unit-cell charge density field ${\displaystyle \rho _{uc}(\mathbf {r} )}$ is a sum over the positions ${\displaystyle \mathbf {r} _{k}}$ of the charges ${\displaystyle q_{k}}$ in the central unit cell $${\displaystyle \rho _{uc}(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{\mathrm {charges} \ k}q_{k}\delta (\mathbf {r} -\mathbf {r} _{k})}$$ and the total charge density field ${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )}$ is the same sum over the unit-cell charges ${\displaystyle q_{k}}$ and their periodic images $${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{n_{1},n_{2},n_{3}}\sum _{\mathrm {charges} \ k}q_{k}\delta (\mathbf {r} -\mathbf {r} _{k}-n_{1}\mathbf {a} _{1}-n_{2}\mathbf {a} _{2}-n_{3}\mathbf {a} _{3})}$$

Here, ${\displaystyle \delta (\mathbf {x} )}$ is the Dirac delta function, ${\displaystyle \mathbf {a} _{1}}$, ${\displaystyle \mathbf {a} _{2}}$ and ${\displaystyle \mathbf {a} _{3}}$ are the lattice vectors and ${\displaystyle n_{1}}$, ${\displaystyle n_{2}}$ and ${\displaystyle n_{3}}$ range over all integers. The total field ${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )}$ can be represented as a convolution of ${\displaystyle \rho _{uc}(\mathbf {r} )}$ with a lattice function ${\displaystyle L(\mathbf {r} )}$ $${\displaystyle L(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{n_{1},n_{2},n_{3}}\delta (\mathbf {r} -n_{1}\mathbf {a} _{1}-n_{2}\mathbf {a} _{2}-n_{3}\mathbf {a} _{3})}$$

Since this is a convolution, the Fourier transformation of ${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )}$ is a product $${\displaystyle {\tilde {\rho }}_{\text{TOT}}(\mathbf {k} )={\tilde {L}}(\mathbf {k} ){\tilde {\rho }}_{uc}(\mathbf {k} )}$$ where the Fourier transform of the lattice function is another sum over delta functions $${\displaystyle {\tilde {L}}(\mathbf {k} )={\frac {\left(2\pi \right)^{3}}{\Omega }}\sum _{m_{1},m_{2},m_{3}}\delta (\mathbf {k} -m_{1}\mathbf {b} _{1}-m_{2}\mathbf {b} _{2}-m_{3}\mathbf {b} _{3})}$$ where the reciprocal space vectors are defined ${\displaystyle \mathbf {b} _{1}\ {\stackrel {\mathrm {def} }{=}}\ 2\pi {\frac {\mathbf {a} _{2}\times \mathbf {a} _{3}}{\Omega }}}$ (and cyclic permutations) where ${\displaystyle \Omega \ {\stackrel {\mathrm {def} }{=}}\ \mathbf {a} _{1}\cdot \left(\mathbf {a} _{2}\times \mathbf {a} _{3}\right)}$ is the volume of the central unit cell (if it is geometrically a parallelepiped, which is often but not necessarily the case). Note that both ${\displaystyle L(\mathbf {r} )}$ and ${\displaystyle {\tilde {L}}(\mathbf {k} )}$ are real, even functions.

For brevity, define an effective single-particle potential $${\displaystyle v(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \int d\mathbf {r} ^{\prime }\,\rho _{uc}(\mathbf {r} ^{\prime })\ \varphi _{\ell r}(\mathbf {r} -\mathbf {r} ^{\prime })}$$

Since this is also a convolution, the Fourier transformation of the same equation is a product $${\displaystyle {\tilde {V}}(\mathbf {k} )\ {\stackrel {\mathrm {def} }{=}}\ {\tilde {\rho }}_{uc}(\mathbf {k} ){\tilde {\Phi }}(\mathbf {k} )}$$ where the Fourier transform is defined $${\displaystyle {\tilde {V}}(\mathbf {k} )=\int d\mathbf {r} \ v(\mathbf {r} )\ e^{-i\mathbf {k} \cdot \mathbf {r} }}$$

The energy can now be written as a single field integral $${\displaystyle E_{\ell r}=\int d\mathbf {r} \ \rho _{\text{TOT}}(\mathbf {r} )\ v(\mathbf {r} )}$$