next up previous
Next: Ewald Forces Up: Long-Range Interactions: The Ewald Previous: Long-Range Interactions: The Ewald

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:


\begin{displaymath}
\mathscr{U}_{\rm Coul} = \frac{1}{2}\sum_{i=1}^N q_i \phi(r_i)
\end{displaymath} (310)

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

\begin{displaymath}
\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}
\end{displaymath} (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:

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:

\begin{displaymath}
\rho_{\rm s}(r) = -q_i(\alpha/\pi)^{3/2} e^{-\alpha r^2}
\end{displaymath} (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:

\begin{displaymath}
-\nabla^2\phi(r) = 4\pi\rho(r)
\end{displaymath} (313)

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

\begin{displaymath}
\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]
\end{displaymath} (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:

\begin{displaymath}
k^2 \tilde{\phi}(k) = 4\pi\tilde{\rho}(k)
\end{displaymath} (315)

The Fourier transform of $\rho_1(r)$ is given by
$\displaystyle \rho_1(k)$ $\textstyle =$ $\displaystyle \int_V d{\bf r} e^{-i{\bf k}\cdot{\bf r}}\rho({\bf r})$ (316)
  $\textstyle =$ $\displaystyle \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
\begin{displaymath}
{\bf k} = \frac{2\pi}{L} {\bf l}     {\bf l} \in \mathbb{Z}^3
\end{displaymath} (318)

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

\begin{displaymath}
\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}
\end{displaymath} (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

\begin{displaymath}
\phi_i(r) = \frac{1}{V}\sum_{{\bf k}\ne 0}\tilde\phi(k)e^{i{\bf k}\cdot{\bf r}}
\end{displaymath} (320)

which, when we substitute for $\tilde{\phi}(k)$ from Eq. 320 yields
\begin{displaymath}
\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}
\end{displaymath} (321)

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

$\displaystyle \mathscr{U}_1$ $\textstyle =$ $\displaystyle \frac{1}{2}\sum_i q_i\phi_1(r_i)$ (322)
  $\textstyle =$ $\displaystyle \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)
  $\textstyle =$ $\displaystyle \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
\begin{displaymath}
\rho(k) = \sum_{i=1}^{N}q_ie^{i{\bf k}\cdot{\bf r}}
\end{displaymath} (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):

\begin{displaymath}
-\frac{1}{r}\frac{\partial^2r\phi_{\rm Gauss}}{\partial r^2} = 4\pi\rho_{\rm Gauss}(r)
\end{displaymath} (326)

yielding
\begin{displaymath}
\phi_{\rm Gauss}(r) = \frac{q_i}{r}{\rm erf}(\sqrt{a}r)
\end{displaymath} (327)

where ${\rm erf}$ is the error function:
\begin{displaymath}
{\rm erf}(x) = \frac{2}{\sqrt{\pi}}\int_0^x e^{-r^2}dr
\end{displaymath} (328)

At $r=0$, we have
\begin{displaymath}
\phi_{\rm self} = \phi_{\rm Gauss}(0) = 2\left(\frac{\alpha}{\pi}\right)^{1/2}q_i
\end{displaymath} (329)

So the total self-interaction energy becomes

\begin{displaymath}
\mathscr{U}_{\rm self} = \left(\frac{\alpha}{\pi}\right)^{1/2}\sum_{i=1}^{N}q_i^2
\end{displaymath} (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:

\begin{displaymath}
\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)
\end{displaymath} (331)

where ${\rm erfc}$ is the complementary error function. The total real-space Coulombic potential energy is therefore
\begin{displaymath}
\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)
\end{displaymath} (332)

Putting it all together:

$\displaystyle \mathscr{U}_{\rm Coul}$ $\textstyle =$ $\displaystyle \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)
    $\displaystyle -\left(\frac{\alpha}{\pi}\right)^{1/2}\sum_{i=1}^{N}q_i^2$ (334)
    $\displaystyle +\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
\begin{displaymath}
\rho(k) = \sum_{i=1}^{N}q_ie^{i{\bf k}\cdot{\bf r}}.
\end{displaymath} (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.


next up previous
Next: Ewald Forces Up: Long-Range Interactions: The Ewald Previous: Long-Range Interactions: The Ewald
cfa22@drexel.edu