The Ewald Coulombic energy

First, assume we have a collection of charged particles in a cubic box with side length $L$, with periodic boundary conditions. The collection is assumed neutral; there is an equal number of positive and negative charges. The total Coulombic energy in this system is given by:

$$\mathscr{U}_{\rm Coul} = \frac{1}{2}\sum_{i=1}^N q_i \phi(r_i)$$ (226)

$\phi(r_i)$ is the electrostatic potential at position $r_i$:

$$\phi\left(r_i\right) = {\sum_{j=1}^{N}}^\prime{\sum_{{\bf n}\in\mathbb{Z}^3}} \frac{q_j}{\left\vert{\bf r}_{ij} + {\bf n}L\right\vert}$$ (227)

${\bf n}$ is a three dimensional integer vector. The prime on the first summation indicates that we do not admit the term for which $j=i$ if ${\bf n} = (0,0,0)$. That is, we allow each particle to interact with its periodic images, but not with itself.

To evaluate $\mathscr{U}$ efficiently, we break it into two parts:

• A short-ranged potential treated with a simple cutoff;
• A long-ranged potential which is periodic and slowly varying, which can therefore be represented to an acceptable level of accuracy by a finite Fourier series.

How can we do this? The idea of Ewald is to do two things: first, screen each point charge using a diffuse cloud of opposite charge around each point charge, and then compensate for these screening charges using a smoothly varying, periodic charge density. The screening charge is constructed to make the electrostatic potential due to a charge at position $r_j$ decay rapidly to near zero at a prescribed distance. These interactions are treated in real space. The compensating charge density, which is the sum of all screening densities except with opposite charges, is treated using a Fourier series.

The standard choice for a screening potential is Gaussian:

$$\rho_{\rm s}(r) = -q_i(\alpha/\pi)^{3/2} e^{-\alpha r^2}$$ (228)

So for each charge, we add such a screening potential. Now, to evaluate $\mathscr{U}_{\rm Coul}$, we have to evaluate the potential of a charge density that compensates for the screening charge densities at each particle. This is done in Fourier space.

The potential of a given charge distribution is given by Poisson's equation:

$$-\nabla^2\phi(r) = 4\pi\rho(r)$$ (229)

Now, the compensating charge distribution, denoted $\rho_1$, can be written:

$$\rho_1(r) = \sum_{j=1}^{N}\sum_{{\bf n}\in\mathbb{Z}^3} q_j (\alp... ...[-\alpha\left\vert{\bf r} - \left({\bf r}_j + {\bf n}L\right)\right\vert\right]$$ (230)

Notice that the sum over $j$ includes the self-interaction when we include the potential due to this charge density in the calcuation of the total Coulombic energy (i.e., we have omitted the prime on the outer summation over particle indices).

Now, consider the Fourier transform of Poisson's equation:

$$k^2 \tilde{\phi}(k) = 4\pi\tilde{\rho}(k)$$ (231)

The Fourier transform of $\rho_1(r)$ is given by

$$\rho_1(k) = \int_V d{\bf r} e^{-i{\bf k}\cdot{\bf r}}\rho({\bf r})$$ (232)

$$= \sum_{j=1}^{N}q_j e^{-i{\bf k}\cdot{\bf r}_j}e^{-k^2/4\alpha}$$ (233)

(The math required to show this involves noting that the the Fourier transform of a Gaussian is another Gaussian, and that the integral over all space of a normalized Gaussian is unity.) The ${\bf k}$-vectors are given by

$${\bf k} = \frac{2\pi}{L} {\bf l} {\bf l} \in \mathbb{Z}^3$$ (234)

We can use Eq. 229 to solve for $\tilde{\phi}(k)$:

$$\tilde{\phi}(k) = \frac{4\pi}{k^2} \sum_{j=1}^{N}q_j e^{-i{\bf k}\cdot{\bf r}_j}e^{-k^2/4\alpha}$$ (235)

Note that this solution is not defined for $k = 0$. In fact, we have to assume that $\tilde\phi(0) = 0$, which is consistent with the notion that our system and all its periodic images is embedded in a medium of infinite dielectric constant (a perfect conductor; the “tinfoil” boundary condition).

Fourier inverting $\tilde\phi(k)$ gives

$$\phi_i(r) = \frac{1}{V}\sum_{{\bf k}\ne 0}\tilde\phi(k)e^{i{\bf k}\cdot{\bf r}}$$ (236)

which, when we substitute for $\tilde{\phi}(k)$ from Eq. 235 yields

$$\phi_i(r) = \sum_{{\bf k}\ne 0}\sum_{j=1}^{N} \frac{4\pi q_j}{Vk^2} e^{i{\bf k}\cdot\left({\bf r}-{\bf r}_j\right)}e^{-k^2/4\alpha}$$ (237)

So, the total Coulombic energy due to the compensating charge distribution is

$$\mathscr{U}_1 = \frac{1}{2}\sum_i q_i\phi_1(r_i)$$ (238)

$$= \frac{1}{2}\sum_{{\bf k}\ne 0}\sum_{j=1}^{N} \frac{4\pi q_i q_j}{Vk^2} e^{i{\bf k}\cdot\left({\bf r}-{\bf r}_j\right)}e^{-k^2/4\alpha}$$ (239)

$$= \frac{1}{2V}\sum_{{\bf k}\ne 0}\frac{4\pi}{k^2}\left\vert\rho(k)\right\vert^2e^{-k^2/4\alpha}$$ (240)

where

$$\rho(k) = \sum_{i=1}^{N}q_ie^{i{\bf k}\cdot{\bf r}}$$ (241)

Notice that this does indeed include a spurious self-self interaction, because the point charge at ${\bf r}_i$ interacts with the compensating charge cloud also at ${\bf r}_i$. This self-interaction is the potential at the center of a Gaussian charge distribution. First, we solve Poisson's equation for the potential due to a Gaussian charge distribution (details in F&S):

$$-\frac{1}{r}\frac{\partial^2r\phi_{\rm Gauss}}{\partial r^2} = 4\pi\rho_{\rm Gauss}(r)$$ (242)

yielding

$$\phi_{\rm Gauss}(r) = \frac{q_i}{r}{\rm erf}(\sqrt{a}r)$$ (243)

where ${\rm erf}$ is the error function:

$${\rm erf}(x) = \frac{2}{\sqrt{\pi}}\int_0^x e^{-r^2}dr$$ (244)

At $r=0$, we have

$$\phi_{\rm self} = \phi_{\rm Gauss}(0) = 2\left(\frac{\alpha}{\pi}\right)^{1/2}q_i$$ (245)

So the total self-interaction energy becomes

$$\mathscr{U}_{\rm self} = \left(\frac{\alpha}{\pi}\right)^{1/2}\sum_{i=1}^{N}q_i^2$$ (246)

which must be subtracted from the total Coulombic energy.

Finally, the real-space contribution of the point charge at ${\bf r}_i$ is the screened potential:

$$\phi_{\rm short}(r) = \frac{q_i}{r} - \frac{q_i}{r}{\rm erf}\left(\sqrt{a}r\right) \equiv \frac{q_i}{r}{\rm erfc}\left(\sqrt{\alpha}r\right)$$ (247)

where ${\rm erfc}$ is the complementary error function. The total real-space Coulombic potential energy is therefore

$$\mathscr{U}_{\rm short} = \frac{1}{2}\sum_{i\ne j}^{N} \frac{q_i q_j}{r_{ij}} {\rm erfc}\left(\sqrt{\alpha}r_{ij}\right)$$ (248)

Putting it all together:

$$\mathscr{U}_{\rm Coul} = \frac{1}{2V}\sum_{{\bf k}\ne 0}\frac{4\pi}{k^2}\left\vert\rho(k)\right\vert^2e^{-k^2/4\alpha}$$ (249)

$$-\left(\frac{\alpha}{\pi}\right)^{1/2}\sum_{i=1}^{N}q_i^2$$ (250)

$$+\frac{1}{2}\sum_{i\ne j}^{N} \frac{q_i q_j}{r_{ij}} {\rm erfc}\left(\sqrt{\alpha}r_{ij}\right)$$ (251)

where

$$\rho(k) = \sum_{i=1}^{N}q_ie^{i{\bf k}\cdot{\bf r}}.$$ (252)

Now, the arbitrariness left to us at this point is in a choice for the parameter $\alpha$. Clearly, very small alphas make the Gaussians tighter and therefore the compensating charge distribution less smoothly varying. This means a Fourier series representation of $\mathscr {U}_1$ with a given number of terms is more accurate for larger $\alpha$. We'll evaluate choice of $\alpha$ in Sec. 6.3.