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)$$ (310)

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

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

${\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$ decays 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}$$ (312)

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)$$ (313)

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 ... ...t{\bf r} - \left({\bf r}_j + {\bf n}L\right)\right\vert\right]$$ (314)

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.

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

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

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})$$ (316)

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

(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$$ (318)

We can use Eq. 314 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}$$ (319)

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}}$$ (320)

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

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

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)$$ (322)

$$= \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}$$ (323)

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

where

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

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 Poission'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)$$ (326)

yielding

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

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

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

At $r=0$, we have

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

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$$ (330)

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}... ...}r\right) \equiv \frac{q_i}{r}{\rm erfc}\left(\sqrt{a}r\right)$$ (331)

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{a}r_{ij}\right)$$ (332)

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}$$ (333)

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

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

where

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

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. 7.3.3.