Nucleus energy
Denote the atomic point charge within the unit cell as
(1)\[\begin{equation}
\rho ^{\text{atom}}(\vb{r})=-\sum_{\ell} Z_{\ell }\delta (\vb{r}-\vb*{\tau }_{\ell })
\end{equation}\]
The Coulombic potential generated from this charge is
(2)\[\begin{equation}
\begin{split}
V_{\text{ext}}(\vb{r}) =& \rho ^{\text{atom}} \star \frac{1}{r}
= \int_{\Omega + \vb{R} } \dd{\vb{r}'} \frac{1}{\norm{\vb{r} - \vb{r}'}} \left( -\sum_{\ell} Z_{\ell }\delta (\vb{r}'-\vb*{\tau }_{\ell }) \right) \\
=& - \sum_{\vb{R}} \sum_{\ell} Z_{\ell } \frac{1}{\norm{\vb{r} - \vb*{\tau }_{\ell } - \vb{R}}}
\end{split}
\end{equation}\]
The nucleus energy is then
(3)\[\begin{equation}
E_{\text{nuc}}
= \frac{1}{2} \braket{V_{\text{ext}}}{\rho ^{\text{atom}}}
=- \frac{1}{2}\sum_{\vb{R}} \sum_{\ell'} \sum_{\ell} Z_{\ell '}Z_{\ell } \frac{1}{\norm{\vb*{\tau }_{\ell' } - \vb*{\tau }_{\ell } - \vb{R}}}
\end{equation}\]
where the \(1 / 2\) factor prevents double counting.
There are two problem with the above formula:
In practical calculation the summation over \(\vb{R}\) is truncated to a finite sized Bravais lattice, which need to be very big for the energy to converge
An uniform negative background charge \(\rho^{-}=\vb{Z}_{\text{tot}} / \Omega\) is required for the energy to be convergent, which is hard to apply in real space.
Reciprocal space representation via the Yukawa kernel
The Fourier space representation of the electrostatic potential can then be determined by the Poisson equation in Fourier space:
(4)\[\begin{equation}
\nabla ^{2} V(\vb{r}) = -4\pi \rho (\vb{r}) \Rightarrow
-\norm{\vb{G}}^{2} \tilde{V} (\vb{G}) = -4\pi \tilde{\rho} (\vb{G}) \Rightarrow
\tilde{V} (\vb{G}) = \frac{4\pi}{\norm{\vb{G}}^{2}} \tilde{\rho} (\vb{G}).
\end{equation}\]
Note that at \(\vb{G}=\vb{0}\) we have a singularity, so this ill-defined.
To fix this, we introduce the Yukawa kernel \(\nu_{\alpha }(r)=\frac{e^{-\alpha r}}{r}\), which represents a screened Coulomb interaction.
It’s Fourier representation is given by
(5)\[\begin{equation}
\tilde{\nu} ^{\alpha }(\vb{G})= \int_{\mathbb{R}^{3}} \dd{\vb{r}} \nu_{\alpha }(\vb{r}) e^{-\text{i} \vb{G} \cdot \vb{r}} = \dfrac{4\pi}{\norm{\vb{G}}^2 + \alpha^2}.
\end{equation}\]
This is equivalent to treat Fourier transform as the limit of the Laplace transform.
Now for any charge distribution \(\rho\), we can define the Fourier transform of Coulombic potentials generated from \(\rho\) as a limit, since \(\lim_{\alpha \to 0} -\frac{1}{4\pi} \nu _{\alpha }=\frac{1}{r}\):
(6)\[\begin{equation}
\tilde{V} (\vb{G})=\lim_{\alpha \to 0} [\widetilde{-\frac{1}{4\pi}\nu_{\alpha } \star -4\pi \rho }] (\vb{G}) =\lim_{\alpha \to 0} \tilde{\nu}_{\alpha } (\vb{G}) \tilde{\rho} (\vb{G})
= \lim_{\alpha \to 0} \frac{4\pi \tilde{\rho} (\vb{G})}{\norm{\vb{G}}^{2}+ \alpha ^{2}}.
\end{equation}\]
Note that this formula implies that the Coulombic potential generated by a neutral charge distribution has zero average, since \(\tilde{\rho}(\vb{0})=0\) implies \(\tilde{V}(\vb{0})=0\).
\(V_{\text{ext}}\) can now be represented in the reciprocal space via the Yukawa kernel:
(7)\[\begin{equation}
\tilde{V}_{\text{ext}} (\vb{G}) =\lim_{\alpha \to 0} \tilde{\nu}_{\alpha } (\vb{G}) \tilde{\rho}^{\text{atom}} (\vb{G})
= -\lim_{\alpha \to 0} \sum_{\ell} \frac{4\pi Z_{\ell } e^{-\text{i} \vb{G} ^{\top} \vb*{\tau }_{\ell}}}{\norm{\vb{G}}^{2}+ \alpha ^{2}}
= -\sum_{\ell} \frac{ 4\pi Z_{\ell } e^{-\text{i} \vb{G} ^{\top} \vb*{\tau }_{\ell}}}{\norm{\vb{G}}^{2}}.
\end{equation}\]
Now since
(8)\[\begin{equation}
\tilde{\rho} ^{\text{atom}}(\vb{G}) = - \sum_{\ell } Z_{\ell } e^{-\text{i} \vb{G}^{\top} \vb*{\tau }_{\ell }}
\end{equation}\]
by Parseval’s theorem we have
(9)\[\begin{equation}
E_{\text{nuc}} = \sum_{\vb{G}\neq \vb{0}} \tilde{V} _{\text{ext}}(\vb{G})^{*} \tilde{\rho} ^{\text{atom}}(\vb{G})
= \sum_{\vb{G}\neq \vb{0}} \sum_{\ell ,\ell '} \frac{1}{\norm{\vb{G}}^{2}} e^{-\text{i} \vb{G}^{\top} (\vb*{\tau }_{\ell } - \vb*{\tau }_{\ell '})}
\end{equation}\]
where the exclusion of the \(\vb{G}=\vb{0}\) term is equivalent to applying the uniform background charge. However, it still converges slowly in terms of the g grid size.
Ewald summation
The idea of Ewald summation is to do summation in both the real and the reciprocal space, where the Coulombic kernel is broken up into a short-range term that converges fast in real space, and a long-range term that converges fast in the reciprocal space. Intuitively, why this works is that the long-range terms will by definition have lower frequency.
Specifically, the Coulombic kernel is decomposed as follows
(10)\[\begin{equation}
\frac{1}{r} = \underbrace{\text{erfc}(\eta r) \frac{1}{r}}_{\nu _{\text{sr}}} + \underbrace{\text{erf}(\eta r) \frac{1}{r}}_{\nu _{\text{lr}}}
\end{equation}\]
where \(\nu_{\text{sr}}, \nu _{\text{lr}}\) are the short-range and the long-range kernel respectively.
Long range potential
Obviously \(\nu_{\text{sr}}\) decays exponentially in the real space. We need to show that the long-range kernel decays fast in the reciprocal space.
The long range potential is given by
(11)\[\begin{equation}
V_{\text{lr}}(\vb{r}) = \rho ^{\text{atom}} \star \nu _{\text{lr}} = - \sum_{\ell } \frac{Z_{\ell } \text{erf}(\eta \norm{ \vb{r} - \vb*{\tau }_{\ell }})}{\norm{\vb{r} - \vb*{\tau }_{\ell }}}
\end{equation}\]
The long-range part will be evaluated in the reciprocal space. First, we derive the Fourier transform of the long-range kernel calculated analytically using the same trick as Yukawa, which treats Fourier transform as the limit of Laplace transform:
(12)\[\begin{equation}
\begin{aligned}
\tilde{\nu }^\alpha_{\text{lr}} (\vb{G})= \int_{\mathbb{R}^{3}} \dd{\vb{r}} \nu _{\text{lr}} (\vb{r}) e^{-\alpha r} e^{-\text{i} \vb{G} \cdot \vb{r}}
=& \dfrac{2 \pi }{\text{i} G } \left[ -\frac{e^{(\frac{\alpha+ \text{i} G}{2\eta })^2 }}{\text{i} G + \alpha } + \frac{e^{(\frac{\alpha- \text{i} G}{2\eta })^2 }}{-\text{i} G + \alpha } \right].
\end{aligned}
\end{equation}\]
where \(G=\norm{\vb{G}}\), and in the last equality we have used square completion trick with the Gaussian integrals. For \(\vb{G}\neq 0\), we can take \(\alpha \to 0\) and get
(13)\[\begin{equation}
\tilde{\nu }_{\text{lr}}(\vb{G}) = \lim_{\alpha \to 0} \tilde{\nu }_{\text{lr}}^\alpha(\vb{G})
= \dfrac{2 \pi }{\text{i} G } \left[ \frac{e^{(\frac{ \text{i} G}{2\eta })^2 }}{\text{i} G } - \frac{e^{(\frac{ -\text{i} G}{2\eta })^2 }}{-\text{i} G } \right]
= \dfrac{4 \pi }{ G^2 } \exp {\frac{-G^2}{4\eta^2 } }.
\end{equation}\]
and the long-range potential for is then given by
(14)\[\begin{equation}
\tilde{V} _{\text{lr}}(\vb{G}) = \widetilde{[\nu_{\text{lr}} \star \rho^{\text{atom}}]}(\vb{G}) = \tilde{\nu} _{\text{lr}}(\vb{G})^{*} \tilde{\rho} ^{\text{atom}}(\vb{G})
= \dfrac{4 \pi }{ G^2 } \exp {\frac{-G^2}{4\eta^2 } } e^{-\text{i} \vb{G} ^\top \vb*{\tau }_{\ell }}
\end{equation}\]
which decays exponentially instead of polynomially.
Interpretation as screening by a Gaussian charge
Although we have shown that this decomposition of the Coulomb kernel yields two fast converging terms, we haven’t handle the neutral charge requirement. The perspective on kernel decomposition is not very helpful in this aspect. Instead, we need to interpret the Ewald summation as screening by a Gaussian charge.
Specifically, we will show that Ewald summation split the neutralized atomic charge distribution \(\rho^{\text{atom}}+ \rho ^{-}\) into two part: (a) the atomic charge distribution neutralized by a Gaussian charge distribution \(\rho^{\text{tot}}\); (b) the Gaussian charge distribution neutralized by the uniform background charge \(\rho^{-}\). The Coulombic potential generated by these two charge distributions corresponds to the short-range and the long-range part:
(15)\[\begin{equation}
\begin{split}
V_{\text{nuc}}
=& \frac{1}{r} \star [\rho ^{\text{atom}} + \rho ^{-}] \\
=& \frac{1}{r} \star [(\rho ^{\text{atom}}-\rho ^{\text{tot}}) + (\rho ^{\text{tot}} + \rho ^{-})] \\
=& (\nu_{\text{sr}} \star \rho ^{\text{atom}} - C) + (\nu_{\text{lr}}\star \rho ^{\text{atom}} + \frac{1}{r} \star \rho^{-}).
\end{split}
\end{equation}\]
Long-range part is the Coulombic potential generated via a Gaussian charge distribution
Consider the positive Gaussian charge distribution
(16)\[\begin{equation}
\rho ^{g}(\vb{r}) = -(\pi /\eta^{2} )^{-\frac{3}{2}} \exp (- \eta^{2} \norm{\vb{r}}^{2})
\end{equation}\]
whose Fourier space representation is
(17)\[\begin{equation}
\tilde{\rho } ^{g}(\vb{G}) = - \frac{1}{\Omega } \exp (-\frac{\norm{\vb{G}}^{2}}{4 \eta ^{2} } ).
\end{equation}\]
Now consider the adding such Gaussian charge distribution at each atom’s site in the crystal, the total charge distribution is:
(18)\[\begin{equation} \label{eqn:n-tot-real}
\rho ^{\text{tot}}(\vb{r})
= \sum_{\ell } \sum_{\vb{n}}^{\infty } Z_{\ell } \rho ^{g}(\vb{r}- \vb*{\tau }_{\ell } - \vb{R}_{\vb{n}})
= - (\pi /\eta^{2} )^{-\frac{3}{2}} \sum_{\ell } \sum_{\vb{n}}^{\infty } Z_{\ell } \exp (-\eta ^{2} \norm{\vb{r}- \vb*{\tau }_{\ell } - \vb{R}_{\vb{n}}}^{2})
\end{equation}\]
whose Fourier space representation is
(19)\[\begin{equation} \label{eqn:n-tot-fourier}
\tilde{\rho }^{\text{tot}}(\vb{G})
= -\frac{1}{\Omega }\sum_{\ell } Z_{\ell } e^{-\text{i} \vb{G}\cdot \vb*{\tau }_{\ell }} \exp (-\frac{\norm{\vb{G}}^{2}}{4 \eta ^{2} } )
\end{equation}\]
From the previous section, the Coulombic potential generated by \(\rho ^{\text{tot}}\) is
(20)\[\begin{equation}
\tilde{V}^{\text{tot}}(\vb{G}) = \left\{\begin{array}{cc}
\frac{4\pi}{\norm{\vb{G}}^{2}} \tilde{\rho }^{\text{tot}}(\vb{G}), & \vb{G}\neq \vb{0}, \\
0, & \vb{G}= \vb{0}, \\
\end{array}\right.
\end{equation}\]
therefore for the case of \(\vb{G}\neq \vb{0}\)
(21)\[\begin{equation}
\hat{V}^{\text{tot}}(\vb{G})
=- \frac{1}{\Omega } \sum_{\vb{G}\neq 0} \sum_{\ell } Z_{\ell } \frac{4\pi}{\norm{\vb{G}}^{2}} \exp (-\frac{\norm{\vb{G}}^{2}}{4 \eta ^{2} } ) e^{-\text{i}\vb{G}^{\top} \vb*{\tau }_{\ell }} = \tilde{V} _{\text{lr}}(\vb{G}).
\end{equation}\]
where removing the \(\vb{G}=\vb{0}\) term is equivalent to screening by uniform charge. Therefore
(22)\[\begin{equation}
\nu _{\text{lr}} \star \rho ^{\text{atom}} + \frac{1}{r} \star \rho ^{-} = V_{\text{lr}} = V^{\text{tot}} = \frac{1}{r} \star [\rho^{\text{tot}} + \rho^{-}]
\end{equation}\]
which says the long-range potential is the Coulombic potential generated via a neutralized Gaussian charge distribution \(\rho^{\text{tot}} + \rho ^-\).
Short-range part
For simplicity, consider the case where \(\rho^{\text{atom}}(\vb{r})=\delta (\vb{r})\), i.e. a system where there’s a unit positive charge at the origin. Then
(23)\[\begin{equation}
\nabla ^{2} \left[ \nu_{\text{lr}} \star \rho ^{\text{atom}} \right] = \nabla ^{2} \left[\text{erf}(\eta \norm{\vb{r}}) \frac{1}{ \norm{\vb{r}}}\right]
= -4\pi [-\rho^{g}(\vb{r})].
\end{equation}\]
Plug this into Poisson equation we get
(24)\[\begin{equation}
\nabla ^{2} [\nu _{\text{sr}}\star \rho^{\text{atom}}](\vb{r})
=\nabla ^{2} [(\frac{1}{r} - \nu _{\text{lr}})\star \rho^{\text{atom}}](\vb{r})
=\nabla ^{2} [\frac{1}{r} \star \rho^{\text{atom}}](\vb{r}) - \nabla ^{2} [\nu _{\text{lr}} \star \rho^{\text{atom}}](\vb{r})
= -4\pi [\rho ^{\text{atom}}-\rho ^{g }](\vb{r}) .
\end{equation}\]
It is easy to see that in the general case where \(\rho ^{\text{atom}}(\vb{r})=-\sum_{\ell} Z_{\ell }\delta (\vb{r}-\vb*{\tau }_{\ell })\) we have
(25)\[\begin{equation}
\nabla ^{2} [\nu _{\text{sr}}\star \rho^{\text{atom}}](\vb{r})
= -4\pi [\rho ^{\text{atom}}-\rho ^{\text{tot} }](\vb{r}) .
\end{equation}\]
Now this implies that
(26)\[\begin{equation} \label{eqn:ewald-sr}
\nu_{\text{sr}}\star \rho^{\text{atom}} - C = \frac{1}{r} \star [\rho^{\text{atom}}-\rho^{\text{tot}}]
\end{equation}\]
where the constant \(C\) is the average value of \(G_{\text{sr}}\star \rho^{\text{atom}}\), which is there due to the zero average condition of the potential. The value of \(C\) can be calculated:
(27)\[\begin{equation}
C = \frac{1}{\Omega} \int_{\Omega }\dd{\vb{r}} [\nu _{\text{sr}}\star \rho ^{\text{atom}}](\vb{r})
= -\frac{\pi }{\Omega \eta ^{2}}.
\end{equation}\]
We see that this construction converts the hard-to-apply charge neutral correct in real space to just a zero average condition of the potential.
Now we can calculate the short-range part of the nucleus energy via a real space summation:
(28)\[\begin{equation}
\begin{split}
E_{\text{nuc}}^{\text{sr}} =& \frac{1}{2} \int_{\Omega } \dd{\vb{r}} [\nu _{\text{sr}}\star \rho ^{\text{atom}}-C](\vb{r}) \rho ^{\text{atom}}(\vb{r}). \\
\approx& \frac{1}{2} \sum_{\ell, \ell'} \sum_{\vb{n}}^{L} Z_{\ell }Z_{\ell '} \frac{\text{erfc}(\eta \norm{\vb*{\tau }_{\ell }- \vb*{\tau }_{\ell' } - \vb{R}_{\vb{n}}})}{\norm{\vb*{\tau }_{\ell } - \vb*{\tau }_{\ell' } - \vb{R}_{\vb{n}}}} - \frac{\pi Z_{\text{tot}}^{2}}{2\Omega \eta ^{2}}.
\end{split}
\end{equation}\]
where \(L\) is the Ewald lattice. Furthermore, terms with \(\ell =\ell '\) are excluded when \(\vb{n}=\vb{0}\) to remove self-interaction.
Long-range part
The long-range contribution is
(29)\[\begin{equation}
E_{\text{nuc}}^{\text{lr}} = \frac{1}{2} \int_{\Omega } \dd{\vb{r}} V_{\text{lr}}(\vb{r})
\rho ^{\text{atom}}(\vb{r}) = \frac{1}{2} \sum_{\vb{G}\neq \vb{0}} \tilde{V}_{\text{lr}} (\vb{G})^* \tilde{\rho}^{\text{atom}}(\vb{G}).
\end{equation}\]
It decays fast in the reciprocal space. So with a reasonably small reciprocal mesh \(L'\), we can calculate the long-range contribution to energy:
(30)\[\begin{equation}
\begin{aligned}
E^{\text{lr}}_{\text{nuc}}
\approx& \frac{1}{2\Omega } \sum_{\vb{G}\neq \vb{0}}^{L'} \sum_{\ell } \sum_{\ell' } Z_{\ell } Z_{\ell'} \frac{4\pi}{\norm{\vb{G}}^{2}} \exp (-\frac{\norm{\vb{G}}^{2}}{4 \eta ^{2} } ) e^{\text{i}\vb{G} \cdot (\vb*{\tau }_{\ell '}-\vb*{\tau }_{\ell })} \\
=& \frac{2\pi }{\Omega } \sum_{\vb{G}\neq \vb{0}}^{L'} \frac{1}{\norm{\vb{G}}^{2}} \exp (-\frac{\norm{\vb{G}}^{2}}{4 \eta ^{2} } ) \left[\sum_{\ell } Z_{\ell } e^{\text{i}\vb{G} \cdot \vb*{\tau }_{\ell }} \right]^{2}.
\end{aligned}
\end{equation}\]
We also need to remove the self-interaction for the long-range part like in the short-range part. Specifically,
since
(31)\[\begin{equation}
\nu _{\text{lr}} \star \rho ^{\text{atom}} + \frac{1}{r} \star \rho ^{-} = V_{\text{lr}},
\end{equation}\]
we need to remove the interaction between \(\nu_{\text{lr}} \star \rho^{\text{atom}}\) and the portion of \(\rho^{\text{atom}}\) within \(\Omega\):
(32)\[\begin{equation}
E^{\text{self}}_{\text{nuc}}
= \frac{1}{2} \int_{\Omega}\dd{\vb{r}} [\nu _{\text{lr}} \star \rho ^{\text{atom}}](\vb{r}) \rho ^{\text{atom}}(\vb{r})
= \sum_{\ell } Z_{\ell }^{2} \eta / \sqrt{\pi}.
\end{equation}\]
Summary
The nucleus repulsion energy can now be evaluated using Ewald summation:
(33)\[\begin{equation}
\begin{split}
E_{\text{nuc}}' =
\frac{1}{2} \int_{\Omega } \dd{\vb{r}} V_{\text{nuc}}'(\vb{r}) \rho^{\text{atom}}(\vb{r}) - E^{\text{self}}_{\text{nuc}}
= E^{\text{sr}}_{\text{nuc}} + E^{\text{lr}}_{\text{nuc}} - E^{\text{self}}_{\text{nuc}}
\end{split}
\end{equation}\]
where \(V_{\text{nuc}}',E_{\text{nuc}}'\) means we have incorporated self-interaction correction,
and we have the following decomposition of the nucleus potential
(34)\[\begin{equation} \label{eqn:ewald-pot}
\begin{split}
V_{\text{nuc}}(\vb{r}) =& (\frac{1}{r} \star [\rho^{\text{atom}}+\rho^{-}])(\vb{r}) = [G_{\text{sr}}\star \rho^{\text{atom}} -C + V^{\text{tot}}](\vb{r}) \\
=&
- \sum_{\ell} \sum_{\vb{n}}^{L} Z_{\ell } \frac{\text{erfc}(\eta \norm{\vb{r}- \vb*{\tau }_{\ell } - \vb{R}_{\vb{n}}})}{\norm{\vb{r} - \vb*{\tau }_{\ell} - \vb{R}_{\vb{n}}}}
+ \frac{\pi Z_{\text{tot}}}{\Omega \eta ^{2}}
- \frac{1}{\Omega } \sum_{\vb{G}\neq 0} \sum_{\ell } Z_{\ell } \frac{4\pi}{\norm{\vb{G}}^{2}} \exp (-\frac{\norm{\vb{G}}^{2}}{4 \eta ^{2} } ) e^{\text{i}\vb{G} \cdot (\vb{r}-\vb*{\tau }_{\ell })}.
\end{split}
\end{equation}\]
