本文是「第一性原理的微观计算模拟」系列的第 6 篇（共 11 篇），内容整理自同名书稿第 3 章，公式编号与原书一致。文中引用的文献依据原书参考文献清单整理并列于文末，按原书清单顺序从 1 开始编号。\nocite{*}
平面波-赝势框架下体系的总能
通过 3.3 节中的讨论可知，使用赝势可以显著减少描述本征轨道所需的平面波数量，从而有效地提高计算效率。此外，利用平面波可以轻松地计算动能项，这一点相对于局域轨道基（如 GTO（高斯型轨道）、STO（斯莱特型轨道）等）或混合基（如 LAPW（线性缀加平面波）等）是一个很大的优势。在可以获得高质量赝势的前提下，利用平面波-赝势框架进行 DFT 计算已经成为材料计算和模拟领域的重要研究方向。目前广泛应用的计算软件，如 VASP、CASTEP、PWSCF 等，都采用了这一类方法。在本节中，我们将详细推导这类方法的总能量表达式。
总能表达式的推导
Ihm、Zunger 和 Cohen 在 1979 年提出了倒空间下利用赝势计算体系总能的公式，奠定了目前计算物理中应用广泛的平面波-赝势方法的理论基础\cite{ihm1979momentum}。为简单起见，设体系为单质，由 $N_{\mathrm{cell}}$ 个单胞组成，每个单胞中有一个原子。该体系在 DFT 的框架下体系总能可以表示为
$$ E_{\mathrm{tot}}=T_{\mathrm s}[\rho]+E_{\mathrm H}[\rho]+E_{\mathrm{xc}}[\rho]+E_{\mathrm{Ie}}+E_{\mathrm{II}} \tag{3.378} $$
式中 $E_{\mathrm{xc}}=\int\rho(\boldsymbol r)\varepsilon_{\mathrm{xc}}(\boldsymbol r)\,\mathrm d\boldsymbol r$，$\varepsilon_{\mathrm{xc}}$ 是每电子交换关联能。如果采用赝势 $V_{\mathrm{ps},l}(\boldsymbol{r}-\boldsymbol{R}_{k,I})$ 代表位于 $\boldsymbol{R}_{k,I}$ 终点处的原子核（$k$ 代表单胞序号）对角动量为 $l$ 的电子波函数的作用，则 $E_{\mathrm{Ie}}[\rho]$（I 代表该单胞中的原子，e 表示电子）可以表示为
$$ E_{\mathrm{Ie}}[\rho]=\sum_{i,l,k,I}\int\psi_i^{*}(\boldsymbol{r})V_{\mathrm{ps},l}(\boldsymbol{r}-\boldsymbol{R}_{k,I})\hat{P}_l\psi_i(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r} \tag{3.379} $$
式中：$\hat{P}_l$ 表示将波函数投影到角动量为 $l$ 的电子波函数上的投影算符。根据式（3.378）导出的单电子薛定谔方程为
$$ \left[-\frac12\nabla^2+V_{\mathrm H}(\boldsymbol r)+V_{\mathrm{xc}}(\boldsymbol r) +\sum_{k,I,l}V_{\mathrm{ps},l}(\boldsymbol r-\boldsymbol R_{k,I})\hat P_l\right] \psi_n(\boldsymbol r)=\varepsilon_n\psi_n(\boldsymbol r) \tag{3.380} $$
式中：$\mu_{\mathrm{xc}}(\boldsymbol{r})$ 为交换关联势，${\mu_{\mathrm{xc}}(\boldsymbol r)=\delta E_{\mathrm{xc}}[\rho]/\delta\rho(\boldsymbol r)}$。$u_{\mathrm{xc}}(\boldsymbol{r})$ 具体的函数形式有很多种，Ihm 等人采用了 $X_\alpha$ 方法中的结果，即方程（3.143），则
$$ E_{\mathrm{xc}}(\boldsymbol{r})=\int\mu_{\mathrm{xc}}(\boldsymbol{r})\,\mathrm{d}\rho(\boldsymbol{r})=\frac{3}{4}\mu_{\mathrm{xc}}(\boldsymbol{r})\rho(\boldsymbol{r}) \tag{3.381} $$
当然也可以采取其他形式，如 3.2 节中介绍的各种交换关联势。
采用平面波 $\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}}$ 作为基函数，对于给定的 $\boldsymbol{k}_i$，$\psi_n(\boldsymbol{r})$ 可以展开为
$$ \psi_n(\boldsymbol{r})=\sum_{\boldsymbol{G}}c_{n,\boldsymbol{k}_i}(\boldsymbol{G})\mathrm{e}^{\mathrm{i}(\boldsymbol{k}_i+\boldsymbol{G})\cdot\boldsymbol{r}} \tag{3.382} $$
则式（3.378）和式（3.380）中的各项均需在倒空间内展开，也即需进行傅里叶变换。首先考虑 Hartree 势（有时也称为库仑势）
$$ V_{\mathrm{H}}(\boldsymbol{r})=\int\frac{\rho(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}' \tag{3.383} $$
其傅里叶变换为
$$ \rho(\boldsymbol G)=\frac1\varOmega\int_\varOmega\rho(\boldsymbol r)e^{-\mathrm i\boldsymbol G\cdot\boldsymbol r}\,\mathrm d\boldsymbol r,\qquad V_{\mathrm H}(\boldsymbol G)=\frac{4\pi\rho(\boldsymbol G)}{|\boldsymbol G|^2} \quad(\boldsymbol G\ne0) \tag{3.384} $$
式（3.384）也可以通过泊松方程或者函数卷积的傅里叶变换得到。相应的 Hartree 能为
$$ E_{\mathrm H}=\frac12\iint\frac{\rho(\boldsymbol r)\rho(\boldsymbol r')}{|\boldsymbol r-\boldsymbol r'|}\,\mathrm d\boldsymbol r\,\mathrm d\boldsymbol r' =\frac\varOmega2\sum_{\boldsymbol G\ne0}V_{\mathrm H}(\boldsymbol G)\rho(-\boldsymbol G) \tag{3.385} $$
式中：$\varOmega$ 是体系的总体积。与其相似，交换关联能 $E_{\mathrm{xc}}$ 在倒空间的表达式为
$$ E_{\mathrm{xc}}=\varOmega\sum_{\boldsymbol{G}}\rho(\boldsymbol{G})\varepsilon_{\mathrm{xc}}(-\boldsymbol{G}) \tag{3.386} $$
式中 $\varepsilon_{\mathrm{xc}}(\boldsymbol G)$ 是每电子交换关联能的归一化傅里叶系数；乘积在倒空间中按指标相反的两个系数求和。
动能项比较简单，计算如下：
$$ \begin{aligned} T&=\frac{1}{N_k}\sum_{n,i}\int\psi_n^{*}(\boldsymbol{r})\frac{-\boldsymbol{\nabla}^{2}}{2}\psi_n(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r}\\ &=\sum_{n,i,\boldsymbol{G},\boldsymbol{G}'}\int c_{n,\boldsymbol{k}_i}^{*}(\boldsymbol{G}')\mathrm{e}^{-\mathrm{i}(\boldsymbol{k}_i+\boldsymbol{G}')\cdot\boldsymbol{r}}\frac{-\boldsymbol{\nabla}^{2}}{2}c_{n,\boldsymbol{k}_i}(\boldsymbol{G})\mathrm{e}^{\mathrm{i}(\boldsymbol{k}_i+\boldsymbol{G})\cdot\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r}\\ &=\sum_{n,i,\boldsymbol{G},\boldsymbol{G}'}\frac{|\,\boldsymbol{k}_i+\boldsymbol{G}\,|^{2}}{2}\int c_{n,\boldsymbol{k}_i}^{*}(\boldsymbol{G}')c_{n,\boldsymbol{k}_i}(\boldsymbol{G})\mathrm{e}^{\mathrm{i}(\boldsymbol{G}-\boldsymbol{G}')\cdot\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r}\\ &=\frac{\varOmega}{2}\sum_{n,i,\boldsymbol{G},\boldsymbol{G}'}|\,\boldsymbol{k}_i+\boldsymbol{G}\,|^{2}\delta_{\boldsymbol{G}\boldsymbol{G}'}\,|\,c_{n,\boldsymbol{k}_i}(\boldsymbol{G})\,|^{2}\\ &=\frac{\varOmega}{2}\sum_{n,i,\boldsymbol{G}}|\,\boldsymbol{k}_i+\boldsymbol{G}\,|^{2}\,|\,c_{n,\boldsymbol{k}_i}(\boldsymbol{G})\,|^{2} \end{aligned} \tag{3.387} $$
式中：$N_k$ 是第一布里渊区中 $\boldsymbol{k}$ 点个数。
比较困难的是原子核-电子相互作用能 $E_{\mathrm{Ie}}$。如前所述，原子核的势场由赝势 $V_{\mathrm{ps},l}(\boldsymbol{r}-\boldsymbol{R}_{k,I})$ 表示，则
$$ E_{\mathrm{Ie}}=\sum_{i,k,l}\int\psi_i^{*}(\boldsymbol{r})V_{\mathrm{ps},l}(\boldsymbol{r}-\boldsymbol{R}_{k,I})\hat{P}_l\psi_i(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r} \tag{3.388} $$
该积分是在全空间中进行的。但是考虑到原子赝势 $V_{\mathrm{ps},l}$ 拥有平移对称性，因此对于每个 $\boldsymbol{R}_k$ 都可以做变量代换：$\boldsymbol{r}'_k=\boldsymbol{r}-\boldsymbol{R}_k$。因此，全空间积分可以表示为单胞 $k$ 的积分且对 $k$ 求和。在各个原胞中，$\boldsymbol{r}'$ 的下标可以忽略。因此，式（3.388）可写为
$$ \begin{aligned} E_{\mathrm{Ie}}&=\sum_{n,k,l}\int\psi_n^{*}(\boldsymbol{r})V_{\mathrm{ps},l}(\boldsymbol{r}-\boldsymbol{R}_k)\hat{P}_l\psi_n(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r}\\ &=\frac{\varOmega}{N_k}\sum_{n,i,\boldsymbol{G},\boldsymbol{G}',l}c_{n,\boldsymbol{k}_i}^{*}(\boldsymbol{G})c_{n,\boldsymbol{k}_i}(\boldsymbol{G}')\times\frac{1}{N_{\mathrm{cell}}\varOmega_{\mathrm{cell}}}\sum_{k,I}\int\mathrm{e}^{-\mathrm{i}(\boldsymbol{k}_i+\boldsymbol{G})(\boldsymbol{r}'+\boldsymbol{R}_{k,I})}\\ &\quad\times V_{\mathrm{ps},l}(\boldsymbol{r}')\hat{P}_l\mathrm{e}^{\mathrm{i}(\boldsymbol{k}_i+\boldsymbol{G}')(\boldsymbol{r}'+\boldsymbol{R}_{k,I})}\,\mathrm{d}\boldsymbol{r}'\\ &=\frac{\varOmega}{N_k}\sum_{n,i,\boldsymbol{G},\boldsymbol{G}',l}c_{n,\boldsymbol{k}_i}^{*}(\boldsymbol{G})c_{n,\boldsymbol{k}_i}(\boldsymbol{G}')\sum_{I}\mathrm{e}^{(\boldsymbol{G}'-\boldsymbol{G})\cdot\boldsymbol{R}_I}\cdot\frac{1}{\varOmega_{\mathrm{cell}}}\int\mathrm{e}^{-\mathrm{i}(\boldsymbol{k}_i+\boldsymbol{G})\boldsymbol{r}}V_{\mathrm{ps},l}(\boldsymbol{r})\hat{P}_l\mathrm{e}^{\mathrm{i}(\boldsymbol{k}_i+\boldsymbol{G}')\cdot\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r}\\ &=\frac{\varOmega}{N_k}\sum_{n,i,\boldsymbol{G},\boldsymbol{G}',l}c_{n,\boldsymbol{k}_i}^{*}(\boldsymbol{G})c_{n,\boldsymbol{k}_i}(\boldsymbol{G}')S(\boldsymbol{G}'-\boldsymbol{G})V_{\mathrm{ps},l,\boldsymbol{k}_i+\boldsymbol{G},\boldsymbol{k}_i+\boldsymbol{G}'} \end{aligned} \tag{3.389} $$
式中：$S(\boldsymbol{G}'-\boldsymbol{G})$ 称为结构因子（注意对原子 $I$ 的求和只限制在一个单胞中）；$V_{\mathrm{ps},l,\boldsymbol{k}_i+\boldsymbol{G},\boldsymbol{k}_i+\boldsymbol{G}'}$ 称为位势因子。位势因子的数学推导比较烦琐，需要将平面波用球谐函数和球贝塞尔函数展开。下面给出具体推导过程。
考虑恒等式
$$ \mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}}=4\pi\sum_{L}i^{l}\mathrm{j}_l(kr)\mathrm{Y}_L^{*}(\hat{\boldsymbol{k}})\mathrm{Y}_L(\hat{\boldsymbol{r}}) \tag{3.390} $$
式中：$\mathrm{j}_l(kr)$ 是球贝塞尔函数；$\mathrm{Y}_L(\hat{\boldsymbol{r}})$ 是球谐函数（见附录 A.3 节）；$L=(l,m)$，$l$ 是角动量量子数，$m$ 是角动量 $z$ 分量量子数；$\hat{\boldsymbol{r}}$ 代表单位矢量。此外，因为式（3.390）左端项与球谐函数中的变量 $\phi$ 无关，所以又有以下恒等式
$$ \mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}}=\sum_{l}(2l+1)i^{l}\mathrm{j}_l(kr)\mathrm{P}_l(\cos\theta) \tag{3.391} $$
由方程（3.389）可知，$V_{\mathrm{ps},l,\boldsymbol{k}_i+\boldsymbol{G},\boldsymbol{k}_i+\boldsymbol{G}'}$ 涉及三个矢量，分别为 $\boldsymbol{r}$、$\boldsymbol{k}_i+\boldsymbol{G}$ 和 $\boldsymbol{k}_i+\boldsymbol{G}'$。设 $\boldsymbol{k}_i+\boldsymbol{G}'$ 为 $z$ 轴，$\boldsymbol{k}_i+\boldsymbol{G}$ 与其夹角为 $\gamma$，$\boldsymbol{r}$ 与其夹角为 $\theta$，则利用式（3.391）展开 $\mathrm{e}^{\mathrm{i}(\boldsymbol{k}_i+\boldsymbol{G}')\cdot\boldsymbol{r}}$，而用式（3.390）展开 $\mathrm{e}^{\mathrm{i}(\boldsymbol{k}_i+\boldsymbol{G})\cdot\boldsymbol{r}}$，可得
$$ \begin{aligned} V_{\mathrm{ps},l,\boldsymbol{k}_i+\boldsymbol{G},\boldsymbol{k}_i+\boldsymbol{G}'}={}&\int\mathrm{e}^{-\mathrm{i}(\boldsymbol{k}_i+\boldsymbol{G}\cdot\boldsymbol{r})}V_{\mathrm{ps},l}\hat{P}_l\sum_{l'}(2l'+1)\mathrm{i}^{l'}\mathrm{j}_{l'}(|\,\boldsymbol{k}_i+\boldsymbol{G}'\,|\,r)\mathrm{P}_{l'}(\cos\theta)\,\mathrm{d}\boldsymbol{r}\\ ={}&4\pi(2l+1)\mathrm{i}^{l}\sum_{L'}(-\mathrm{i})^{l'}\mathrm{j}_{l'}(|\,\boldsymbol{k}_i+\boldsymbol{G}\,|\,r)\mathrm{Y}_{l'}^{-m'}(\gamma,\phi)\mathrm{Y}_{l'}^{-m'}(\theta,\phi)V_{\mathrm{ps},l}(\boldsymbol{r})\\ &\times\mathrm{j}_l(|\,\boldsymbol{k}_i+\boldsymbol{G}'\,|\,r)\,\mathrm{d}\boldsymbol{r}\\ ={}&4\pi\,(2l+1)\mathrm{i}^{l}\sum_{L'}\int(-\mathrm{i})^{l'}\mathrm{j}_{l'}(|\,\boldsymbol{k}_i+\boldsymbol{G}\,|\,r)\mathrm{j}_l(|\,\boldsymbol{k}_i+\boldsymbol{G}'\,|\,r)V_{\mathrm{ps},l}(r)\\ &\cdot\mathrm{Y}_{l'}^{m'}(\gamma,\phi)r^{2}\,\mathrm{d}r\cdot\left(\frac{4\pi}{2l+1}\right)^{1/2}\int\mathrm{Y}_{l'}^{-m'}(\theta,\phi)\mathrm{Y}_l^{0}(\theta,\phi)\sin\theta\,\mathrm{d}\theta\,\mathrm{d}\phi\\ ={}&4\pi(2l+1)\mathrm{i}^{l}\sum_{L'}\int(-\mathrm{i})^{l'}\mathrm{j}_{l'}(|\,\boldsymbol{k}_i+\boldsymbol{G}\,|\,r)\mathrm{j}_l(|\,\boldsymbol{k}_i+\boldsymbol{G}'\,|\,r)V_{\mathrm{ps},l}(r)\\ &\cdot\sqrt{\frac{(2l'+1)(l'+m')!}{(2l+1)(l'-m')!}}\mathrm{P}_{l'}^{-m'}(\cos\gamma)\times\mathrm{e}^{-\mathrm{i}m'\phi}r^{2}\,\mathrm{d}r\,\delta_{ll',m'0}\\ ={}&4\pi(2l+1)\int\mathrm{j}_l(|\,\boldsymbol{k}_i+\boldsymbol{G}\,|\,r)\mathrm{j}_l(|\,\boldsymbol{k}_i+\boldsymbol{G}'\,|\,r)V_{\mathrm{ps},l}(r)\mathrm{P}_l(\cos\gamma)r^{2}\,\mathrm{d}r \end{aligned} \tag{3.392} $$
对于原子赝势，我们采取第 2 章中介绍过的半局域分部形式。其中球对称的局域部分为 $V_{\mathrm{ps}}^{\mathrm{loc}}(\boldsymbol{r}-\boldsymbol{R}_{k,I})$，则其对总能的贡献有很简单的形式：
$$ \sum_{i,k,I}\langle\psi_i|V_{\mathrm{ps}}^{\mathrm{loc}}(\boldsymbol r-\boldsymbol R_{k,I})|\psi_i\rangle =\varOmega\sum_{\boldsymbol G}S(\boldsymbol G)V_{\mathrm{ps}}^{\mathrm{loc}}(\boldsymbol G)\rho(-\boldsymbol G) \tag{3.393} $$
而非局域部分 $\delta V_{\mathrm{ps},l}^{\mathrm{nl}}$ 则利用方程（3.392）计算，仅需要将该方程中的 $V_{\mathrm{ps},l}$ 替换为 $\delta V_{\mathrm{ps},l}^{\mathrm{nl}}$ 即可。
因此，若取基函数为平面波，则在倒空间中体系的总能可以表示为
$$ \begin{aligned} E_{\mathrm{tot}}={}&\varOmega\left\{\frac{1}{2N_k}\sum_{n,i,\boldsymbol{G}}|\,c_{n,\boldsymbol{k}_i}(\boldsymbol{G})\,|^{2}\,|\,\boldsymbol{k}_i+\boldsymbol{G}\,|^{2}+\sum_{\boldsymbol{G}}\left[\frac{1}{2}V_{\mathrm{H}}(\boldsymbol{G})+\varepsilon_{\mathrm{xc}}(\boldsymbol{G})+V_{\mathrm{ps}}^{\mathrm{loc}}(\boldsymbol{G})S(\boldsymbol{G})\right]\rho(-\boldsymbol{G})\right.\\ &\left.+\frac{1}{N_k}\sum_{n,i,l,\boldsymbol{G},\boldsymbol{G}'}c_{n,\boldsymbol{k}_i}^{*}(\boldsymbol{G})c_{n,\boldsymbol{k}_i}(\boldsymbol{G}')\delta V_{\mathrm{ps},l,\boldsymbol{k}_i+\boldsymbol{G},\boldsymbol{k}_i+\boldsymbol{G}'}^{\mathrm{nl}}\right\}+\frac{1}{2}\sum_{k,I,m,J}\frac{Z^{2}}{|\,\boldsymbol{R}_{k,I}-\boldsymbol{R}_{m,J}\,|} \end{aligned} \tag{3.394} $$
将式（3.394）推广到普遍情况（即每个单胞中有 $P_s$ 种原子，每种原子有 $\mathrm{N}_s$ 个）非常简单，只需要额外给结构因子 $\mathrm{S}$ 及赝势 $V_{\mathrm{ps},l}$ 加上上标 $s$ 并对其求和即可：
$$ \begin{aligned} E_{\mathrm{tot}}={}&\varOmega\left\{\frac{1}{2N_k}\sum_{n,i,\boldsymbol{G}}|\,c_{n,\boldsymbol{k}_i}(\boldsymbol{G})\,|^{2}\,|\,\boldsymbol{k}_i+\boldsymbol{G}\,|^{2}+\sum_{\boldsymbol{G}}\left[\frac{1}{2}V_{\mathrm{H}}(\boldsymbol{G})+\varepsilon_{\mathrm{xc}}(\boldsymbol{G})\right.\right.\\ &\left.\left.+\sum_{s=1}^{P_s}V_{\mathrm{ps}}^{\mathrm{loc},s}(\boldsymbol{G})S^{s}(\boldsymbol{G})\right]\rho(-\boldsymbol{G})+\frac{1}{N_k}\sum_{n,i,l,\boldsymbol{G},\boldsymbol{G}'}c_{n,\boldsymbol{k}_i}^{*}(\boldsymbol{G})c_{n,\boldsymbol{k}_i}(\boldsymbol{G}')\sum_{s=1}^{P_s}\delta V_{\mathrm{ps},l,\boldsymbol{k}_i+\boldsymbol{G},\boldsymbol{k}_i+\boldsymbol{G}'}^{\mathrm{nl},s}\right\}\\ &+\frac{1}{2}\sum_{k,I,m,J}\frac{Z_IZ_J}{|\,\boldsymbol{R}_{k,I}-\boldsymbol{R}_{m,J}\,|} \end{aligned} \tag{3.395} $$
其中结构因子为
$$ S^{s}(\boldsymbol{G}'-\boldsymbol{G})=\sum_{I=1}^{N_s}\mathrm{e}^{\mathrm{i}(\boldsymbol{G}'-\boldsymbol{G})\cdot\boldsymbol{R}_{I,s}} \tag{3.396} $$
式中求和遍历第 $s$ 种元素的所有原子，$\boldsymbol{R}_{I,s}$ 为每个该种原子在单胞中的位置。注意：与方程（3.389）中定义的 $S(\boldsymbol{G}'-\boldsymbol{G})$ 不同，式（3.396）中对原子位置的求和限制在一个单胞内，而取消了分母上的单胞数 $N_{\mathrm{cell}}$。这显然是合理的，因为对于任意 $n$，$\mathrm{e}^{\mathrm{i}n\boldsymbol{G}\cdot\boldsymbol{L}}\equiv1$。
在上面的所有计算中，非局域赝势 $\delta V_{\mathrm{ps},l}^{s}$ 需要计算 $\mathrm{N}\ (\mathrm{N}+1)/2$ 个积分，所以计算量比较大。为了解决这个问题，可以利用 3.3.3 节中介绍的 KB 方法将其改写为局域赝势。我们在这里给出最终结果（更详细的讨论可参阅文献\cite{kohanoff2006electronic}）：
$$ \delta V_{\mathrm{ps},l,\boldsymbol{k}_i+\boldsymbol{G},\boldsymbol{k}_i+\boldsymbol{G}'}^{\mathrm{nl},s}=\sum_{m=-l}^{l}\sum_{I=1}^{N_s}\beta_{lm}^{s}\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{R}_{I,s}}f_{lm}^{s*}(\boldsymbol{k}_i+\boldsymbol{G})f_{lm}^{s}(\boldsymbol{k}_i+\boldsymbol{G}')\mathrm{e}^{\mathrm{i}\boldsymbol{G}'\cdot\boldsymbol{R}_{I,s}} \tag{3.397} $$
式中
$$ \beta_{lm}^{s}=\left(\int r^{2}\,\mathrm{d}r\,|\,\varPhi_{\mathrm{ps}}^{lm,s}(r)\,|^{2}\,\delta V_{\mathrm{ps},l}^{\mathrm{nl},s}(r)\right)^{-1} \tag{3.398} $$
$$ f_{lm}^{s}(\boldsymbol K)=4\pi(-\mathrm i)^lY_{lm}^{*}(\widehat{\boldsymbol K}) \int_0^\infty r^2\,\mathrm dr\;R_{\mathrm{ps}}^{l,s}(r) \delta V_{\mathrm{ps}}^{l,s}(r)j_l(|\boldsymbol K|r),\qquad\boldsymbol K=\boldsymbol k_i+\boldsymbol G \tag{3.399} $$
按照式（3.397），计算量减小为 $N$ 个积分。分析式（3.397），可知其包含了结构因子。
综合上述讨论，我们得到平面波-赝势框架下的 Kohn-Sham 方程为
$$ \sum_{\boldsymbol{G}'}\left(\frac{1}{2}\,|\,\boldsymbol{k}_i+\boldsymbol{G}'\,|^{2}\delta_{\boldsymbol{G}\boldsymbol{G}'}+V_{\boldsymbol{G},\boldsymbol{G}'}^{i}\right)c_{n,\boldsymbol{k}_i}(\boldsymbol{G}')=\varepsilon_nc_{n,\boldsymbol{k}_i}(\boldsymbol{G}) \tag{3.400} $$
式中
$$ \begin{aligned} V_{\boldsymbol{G},\boldsymbol{G}'}^{i}={}&V_{\mathrm{H}}(\boldsymbol{G}'-\boldsymbol{G})+\mu_{\mathrm{xc}}(\boldsymbol{G}'-\boldsymbol{G})+\sum_{s=1}^{P_s}S^{s}(\boldsymbol{G}'-\boldsymbol{G})V_{\mathrm{ps}}^{\mathrm{loc},s}(\boldsymbol{G}'-\boldsymbol{G})\\ &+\sum_{s=1}^{P_s}\sum_{l}\delta V_{\mathrm{ps},l,\boldsymbol{k}_i+\boldsymbol{G},\boldsymbol{k}_i+\boldsymbol{G}'}^{\mathrm{nl},s} \end{aligned} \tag{3.401} $$
需要指出的是，式（3.384）中的 $V_{\mathrm{H}}(\boldsymbol{G})$、式（3.393）中的 $V_{\mathrm{ps}}^{\mathrm{loc}}(\boldsymbol{G})$ 均正比于 $1/|\,\boldsymbol{G}\,|^{2}$，因此 $V_{\mathrm{H}}(0)$ 和 $V_{\mathrm{ps}}^{\mathrm{loc}}(0)$ 分别发散。而式（3.378）中 $E_{\mathrm{II}}$ 同样包含这样的发散项。可以严格证明，虽然发散项各自发散，但是因为整个体系呈电中性，因此对这三项求和时发散项可彼此抵消，总能仍然是收敛的。具体的做法是在体系中加上遵循某种分布的负电荷 $\rho_{\mathrm{aux}}(\boldsymbol{r})$，其总量 $\displaystyle\int\rho_{\mathrm{aux}}(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r}$ 等于全体离子所带电荷。然后相应地对体系的总能进行修正，加上并减去 $\rho_{\mathrm{aux}}(\boldsymbol{r})$ 的自相互作用 $E_{\mathrm{aux}}$\cite{kohanoff2006electronic,marx2009initio}。附加电荷的分布是任意的，但是为了计算方便，一般取为以各离子为中心呈高斯分布的电荷的叠加。对这个问题的处理需要用到 Ewald 求和，因此首先对其加以介绍。
Ewald 求和
Ewald 求和方法最早是由 Ewald 提出的，用于计算周期性排列的点电荷势能，或者更准确地说，用于计算静电能\cite{ewald1921berechnung}。其核心思想是将库仑势分为长程势与短程势两部分，分别在倒空间与实空间内计算其各自对静电能的贡献，再对结果求和。Ewald 求和作为一种成熟的、被广泛应用的计算方法，其公式的数学推导有若干不同的方法\cite{gao1998large,lee2009ewald}。本节中我们选择 Stanford 大学 Lee 和 Cai 的推导方法\cite{lee2009ewald}，因为该方法涉及的物理图像比较清晰，过程也比较直观。
设一个单胞周期为 $\boldsymbol{L}$ 的体系，单胞中有 $N_{\mathrm{ion}}$ 个离子，携带电荷数为 $\{Z_I\}$，所处位置为 $\{\boldsymbol{r}_I\}$，则体系的静电能（原子单位制下）为
$$ E_{\mathrm{es}}=\frac{1}{2}\sum_{n}\sum_{I}\sum_{J}{}'\frac{Z_IZ_J}{|\,\boldsymbol{r}_I-\boldsymbol{r}_J+n\boldsymbol{L}\,|} \tag{3.402} $$
式中求和符号上的“$'$”表示不包括 $n=0$，$I=J$ 的一项。同时定义下面两个势函数：
$$ \varphi_I(\boldsymbol{r})=\frac{Z_I}{|\,\boldsymbol{r}-\boldsymbol{r}_I\,|} \tag{3.403} $$
$$ \varphi(\boldsymbol{r})=\sum_{n}\sum_{J}\frac{Z_J}{|\,\boldsymbol{r}-\boldsymbol{r}_J+n\boldsymbol{L}\,|} \tag{3.404} $$
式中：$\varphi_I(\boldsymbol{r})$ 为离子 $I$ 在空间中产生的静电势；$\varphi(\boldsymbol{r})$ 为所有离子及其全部映像在空间中产生的静电势。由此可以定义嵌入势
$$ \varphi_{[I]}(\boldsymbol{r})=\varphi(\boldsymbol{r})-\varphi_I(\boldsymbol{r})=\sum_{n}\sum_{J}{}'\frac{Z_J}{|\,\boldsymbol{r}-\boldsymbol{r}_J+n\boldsymbol{L}\,|} \tag{3.405} $$
式中：$\varphi_{[I]}(\boldsymbol{r})$ 表示当 $n=0$ 的单胞内 $\boldsymbol{r}_I$ 终点处没有离子 $I$ 时空间中所有其他离子所产生的静电势。该势函数的重要性在于可以借助它将 $E_{\mathrm{es}}$ 表示为如下非常简单的形式：
$$ E_{\mathrm{es}}=\frac{1}{2}\sum_{I}Z_I\varphi_{[I]}(\boldsymbol{r}_I) \tag{3.406} $$
上述讨论均基于点电荷模型。借助电荷密度分布函数 $\rho_I(\boldsymbol{r})$，可以将上述讨论扩展到一般情况。容易写出，由 $\rho_I(\boldsymbol{r})$ 产生的静电势为
$$ \varphi_I(\boldsymbol{r})=\int\frac{\rho_I(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}' \tag{3.407} $$
由此，遵循与此前一样的步骤，写出一般情况下的静电势能：
$$ E_{\mathrm{es}}=\frac{1}{2}\sum_{n}\sum_{I}\sum_{J}{}'\iint\frac{\rho_I(\boldsymbol{r})\rho_J(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'+n\boldsymbol{L}\,|}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}' \tag{3.408} $$
而相应地，有
$$ \varphi_{[I]}(\boldsymbol{r})=\sum_{n}\sum_{J}{}'\int\frac{\rho_J(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'+n\boldsymbol{L}\,|}\,\mathrm{d}\boldsymbol{r}' \tag{3.409} $$
为了计算静电能，在每个点电荷 $Z_I\delta(\boldsymbol{r}-\boldsymbol{r}_I)$ 上先叠加再减去一个相同的呈高斯分布的电荷 $\rho_\sigma^{\mathrm{G}}(\boldsymbol{r})$，从而将其分解为短程电荷 $\rho_I^{\mathrm{S}}(\boldsymbol{r})$ 和长程电荷 $\rho_I^{\mathrm{L}}(\boldsymbol{r})$：
$$ \rho_I^{\mathrm{S}}=Z_I\delta(\boldsymbol{r}-\boldsymbol{r}_I)-Z_I\rho_\sigma^{\mathrm{G}}(\boldsymbol{r}-\boldsymbol{r}_I) \tag{3.410} $$
$$ \rho_I^{\mathrm{L}}=Z_I\rho_\sigma^{\mathrm{G}}(\boldsymbol{r}-\boldsymbol{r}_I) \tag{3.411} $$
式中：$\rho_\sigma^{\mathrm{G}}(\boldsymbol{r}-\boldsymbol{r}_I)$ 代表中心在 $\boldsymbol{r}_I$ 终点处、展宽为 $\sigma$ 的高斯电荷分布，且有
$$ \rho_\sigma^{\mathrm{G}}(\boldsymbol{r}-\boldsymbol{r}_I)=\frac{1}{(2\pi\sigma^{2})^{3/2}}\mathrm{e}^{-|\boldsymbol{r}-\boldsymbol{r}_I|^{2}/(2\sigma^{2})} \tag{3.412} $$
相应地有短程势 $\varphi_I^{\mathrm{S}}(\boldsymbol{r})$ 和长程势 $\varphi_I^{\mathrm{L}}(\boldsymbol{r})$：
$$ \varphi_I^{\mathrm{S}}(\boldsymbol{r})=Z_I\int\frac{\delta(\boldsymbol{r}'-\boldsymbol{r}_I)-\rho_\sigma^{\mathrm{G}}(\boldsymbol{r}'-\boldsymbol{r}_I)}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}' \tag{3.413} $$
$$ \varphi_I^{\mathrm{L}}(\boldsymbol{r})=Z_I\int\frac{\rho_\sigma^{\mathrm{G}}(\boldsymbol{r}'-\boldsymbol{r}_I)}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}' \tag{3.414} $$
而
$$ \varphi_I(\boldsymbol{r})=\varphi_I^{\mathrm{S}}(\boldsymbol{r})+\varphi_I^{\mathrm{L}}(\boldsymbol{r}) \tag{3.415} $$
显然，静电能可以写为
$$ E_{\mathrm{es}}=\frac{1}{2}\sum_{I}Z_I\varphi_{[I]}^{\mathrm{S}}(\boldsymbol{r}_I)+\frac{1}{2}\sum_{I}Z_I\varphi_{[I]}^{\mathrm{L}}(\boldsymbol{r}_I) \tag{3.416} $$
式中：$\varphi_{[I]}^{\mathrm{S}}(\boldsymbol{r}_I)$、$\varphi_{[I]}^{\mathrm{L}}(\boldsymbol{r}_I)$ 由方程（3.409）给出，仅需将其中的 $\rho_J$ 相应替换为 $\rho_J^{\mathrm{S}}$ 和 $\rho_J^{\mathrm{L}}$ 即可。在实际计算中，常常将 $E_{\mathrm{es}}$ 写成如下形式：
$$ E_{\mathrm{es}}=\frac{1}{2}\sum_{I}Z_I\varphi_{[I]}^{\mathrm{S}}(\boldsymbol{r}_I)+\frac{1}{2}\sum_{I}Z_I\varphi^{\mathrm{L}}(\boldsymbol{r}_I)-\frac{1}{2}\sum_{I}Z_I\varphi_I^{\mathrm{L}}(\boldsymbol{r}_I) \tag{3.417} $$
式（3.417）右端第一项称为短程静电能，用 $E_{\mathrm{es}}^{\mathrm{S}}$ 表示；第二项称为长程静电能，用 $E_{\mathrm{es}}^{\mathrm{L}}$ 表示；第三项称为自能项，用 $E^{\mathrm{self}}$ 表示。$E_{\mathrm{es}}^{\mathrm{S}}$ 与 $E^{\mathrm{self}}$ 需要在实空间内求解，而 $E_{\mathrm{es}}^{\mathrm{L}}$ 需要在倒空间内求解。
静电势 $\varphi$ 与电荷 $\rho$ 之间的关系由泊松方程确定。原子单位制下，有
$$ \boldsymbol{\nabla}^{2}\varphi_\sigma^{\mathrm{G}}(\boldsymbol{r})=-4\pi\rho_\sigma^{\mathrm{G}}(\boldsymbol{r}) \tag{3.418} $$
附录 A.2 节中给出了球坐标系下 $\boldsymbol{\nabla}^{2}$ 的具体形式。考虑到静电势的球对称性，$\varphi_\sigma^{\mathrm{G}}(\boldsymbol{r})$ 仅由 $r=|\,\boldsymbol{r}\,|$ 决定。因此，$\varphi_\sigma^{\mathrm{G}}(\boldsymbol{r})$ 对 $\theta$ 和 $\phi$ 的导数均为零。因此方程（3.418）可写为
$$ \frac{1}{r}\frac{\partial^{2}}{\partial r^{2}}[r\varphi_\sigma^{\mathrm{G}}(\boldsymbol{r})]=-4\pi\rho_\sigma^{\mathrm{G}}(\boldsymbol{r}) \tag{3.419} $$
由此解得
$$ \varphi_\sigma^{\mathrm{G}}(\boldsymbol{r})=\frac{1}{|\,\boldsymbol{r}\,|}\mathrm{erf}\left(\frac{|\,\boldsymbol{r}\,|}{\sqrt{2}\sigma}\right) \tag{3.420} $$
其中 $\mathrm{erf}(z)$ 为误差函数，有
$$ \mathrm{erf}(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}\mathrm{e}^{-t^{2}}\,\mathrm{d}t \tag{3.421} $$
将式（3.421）代入式（3.413）和式（3.414），可得
$$ \varphi_I^{\mathrm{S}}(\boldsymbol{r})=\frac{Z_I}{|\,\boldsymbol{r}-\boldsymbol{r}_I\,|}\mathrm{erfc}\left(\frac{|\,\boldsymbol{r}-\boldsymbol{r}_I\,|}{\sqrt{2}\sigma}\right) \tag{3.422} $$
$$ \varphi_I^{\mathrm{L}}(\boldsymbol{r})=\frac{Z_I}{|\,\boldsymbol{r}-\boldsymbol{r}_I\,|}\mathrm{erf}\left(\frac{|\,\boldsymbol{r}-\boldsymbol{r}_I\,|}{\sqrt{2}\sigma}\right) \tag{3.423} $$
其中，$\mathrm{erfc}(z)=1-\mathrm{erf}(z)$。将式（3.423）和式（3.409）代入式（3.417），可得
$$ E_{\mathrm{es}}^{\mathrm{S}}=\frac{1}{2}\sum_{n}\sum_{I}\sum_{J}{}'\frac{Z_IZ_J}{|\,\boldsymbol{r}_I-\boldsymbol{r}_J+n\boldsymbol{L}\,|}\mathrm{erfc}\left(\frac{|\,\boldsymbol{r}_I-\boldsymbol{r}_J+n\boldsymbol{L}\,|}{\sqrt{2}\sigma}\right) \tag{3.424} $$
因为 $E^{\mathrm{self}}$ 中 $|\,\boldsymbol{r}_I-\boldsymbol{r}_J\,|=0$，所以利用 $z\to0$ 时 $\mathrm{erf}(z)$ 的极限
$$ \lim_{z\to0}\mathrm{erf}(z)=\frac{2}{\sqrt{\pi}}z $$
可得
$$ \varphi_I^{\mathrm{L}}(\boldsymbol{r}_I)=\frac{Z_I}{\sigma}\sqrt{\frac{2}{\pi}} \tag{3.425} $$
因此，将其代入式（3.417）中，可求得
$$ E^{\mathrm{self}}=\frac{1}{\sqrt{2\pi}\sigma}\sum_{I}Z_I^{2} \tag{3.426} $$
接下来讨论 $E_{\mathrm{es}}^{\mathrm{L}}$ 的计算过程。由式（3.414）可以看出，$\varphi_I^{\mathrm{L}}(\boldsymbol{r})$ 是一个长程且无奇点的势函数，因此 $E_{\mathrm{es}}^{\mathrm{L}}$ 无法直接在实空间内求解。Ewald 借助傅里叶变换及其逆变换，首先将 $\rho^{\mathrm{L}}(\boldsymbol{r})$ 变换为倒空间中的 $\rho^{\mathrm{L}}(\boldsymbol{G})$，然后通过倒空间下的泊松方程求解 $\varphi^{\mathrm{L}}(\boldsymbol{G})$，再变换到实空间中，得到了最后结果。
式（3.417）右端第二项的求和中没有扣除任何离子的贡献，因此 $\varphi^{\mathrm{L}}(\boldsymbol{r})$ 是由在空间中呈周期分布的所有 $\rho_I^{\mathrm{L}}(\boldsymbol{r})$ 叠加而生成的长程势。可以写出长程电荷总和：
$$ \rho^{\mathrm{L}}(\boldsymbol{r})=\sum_{n}\sum_{I}\rho_I^{\mathrm{L}}(\boldsymbol{r}+n\boldsymbol{L}) \tag{3.427} $$
显然，$\rho^{\mathrm{L}}(\boldsymbol{r})$ 与 $\varphi^{\mathrm{L}}(\boldsymbol{r})$ 均为周期性函数。因此，可以定义二者各自的傅里叶变换：
$$ \tilde{\rho}^{\mathrm{L}}(\boldsymbol{G})=\frac{1}{\varOmega_{\mathrm{cell}}}\int_{\varOmega_{\mathrm{cell}}}\rho^{\mathrm{L}}(\boldsymbol{r})\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r} \tag{3.428} $$
$$ \tilde{\varphi}^{\mathrm{L}}(\boldsymbol{G})=\frac{1}{\varOmega_{\mathrm{cell}}}\int_{\varOmega_{\mathrm{cell}}}\varphi^{\mathrm{L}}(\boldsymbol{r})\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r} \tag{3.429} $$
其中积分在单胞 $\varOmega_{\mathrm{cell}}$ 中进行。相应的傅里叶逆变换为
$$ \rho^{\mathrm{L}}(\boldsymbol{r})=\sum_{\boldsymbol{G}}\tilde{\rho}^{\mathrm{L}}(\boldsymbol{G})\mathrm{e}^{\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{r}} \tag{3.430} $$
$$ \varphi^{\mathrm L}(\boldsymbol r) =\sum_{\boldsymbol G}\tilde\varphi^{\mathrm L}(\boldsymbol G)e^{\mathrm i\boldsymbol G\cdot\boldsymbol r} \tag{3.431} $$
将 $\rho^{\mathrm{L}}(\boldsymbol{r})$ 的表达式（3.427）、式（3.412）代入式（3.428），可得
$$ \begin{aligned} \tilde{\rho}^{\mathrm{L}}(\boldsymbol{G})&=\frac{1}{\varOmega_{\mathrm{cell}}}\int_{\varOmega_{\mathrm{cell}}}\sum_{n}\sum_{J}Z_J\rho_\sigma^{\mathrm{G}}(\boldsymbol{r}-\boldsymbol{r}_J+n\boldsymbol{L})\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r}\\ &=\frac{1}{\varOmega_{\mathrm{cell}}}\sum_{J}Z_J\int_{\mathrm{allspace}}\rho_\sigma^{\mathrm{G}}(\boldsymbol{r}-\boldsymbol{r}_J)\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r}\\ &=\frac{1}{\varOmega_{\mathrm{cell}}}\sum_{J}Z_J\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{r}_J}\mathrm{e}^{-\sigma^{2}|\boldsymbol{G}|^{2}/2} \end{aligned} \tag{3.432} $$
倒空间中泊松方程为
$$ |\,\boldsymbol{G}\,|^{2}\tilde{\varphi}^{\mathrm{L}}(\boldsymbol{G})=4\pi\tilde{\rho}^{\mathrm{L}}(\boldsymbol{G}) \tag{3.433} $$
将式（3.432）代入式（3.433），有
$$ \tilde{\varphi}^{\mathrm{L}}(\boldsymbol{G})=\frac{4\pi}{\varOmega_{\mathrm{cell}}|\,\boldsymbol{G}\,|^{2}}\sum_{J}Z_J\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{r}_J}\mathrm{e}^{-\sigma^{2}|\boldsymbol{G}|^{2}/2} \tag{3.434} $$
然后利用傅里叶逆变换得
$$ \begin{aligned} \varphi^{\mathrm{L}}(\boldsymbol{r})&=\sum_{\boldsymbol{G}}\tilde{\varphi}^{\mathrm{L}}(\boldsymbol{G})\mathrm{e}^{\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{r}}\\ &=\frac{4\pi}{\varOmega_{\mathrm{cell}}}\lim_{\boldsymbol{G}_0\to\boldsymbol{0}}\frac{\displaystyle\sum_{J}Z_J}{|\,\boldsymbol{G}_0\,|^{2}}+\frac{4\pi}{\varOmega_{\mathrm{cell}}}\sum_{\boldsymbol{G}\neq\boldsymbol{0}}\frac{1}{|\,\boldsymbol{G}\,|^{2}}\sum_{J}Z_J\mathrm{e}^{-\sigma^{2}|\boldsymbol{G}|^{2}/2}\mathrm{e}^{\mathrm{i}\boldsymbol{G}\cdot(\boldsymbol{r}-\boldsymbol{r}_J)} \end{aligned} \tag{3.435} $$
式（3.435）右端的第一项因为单胞呈电中性，即 $\displaystyle\sum_{J}Z_J=0$，所以为 0，故 $\varphi^{\mathrm{L}}(\boldsymbol{r})$ 仅为其第二项 $\boldsymbol{G}\neq\boldsymbol{0}$ 的求和。将式（3.435）代入式（3.417），可得长程静电能为
$$ E_{\mathrm{es}}^{\mathrm{L}}=\frac{1}{2}\sum_{I}Z_I\varphi^{\mathrm{L}}(\boldsymbol{r}_I)=\frac{2\pi}{\varOmega_{\mathrm{cell}}}\sum_{\boldsymbol{G}\neq\boldsymbol{0}}\sum_{I}\sum_{J}\frac{Z_IZ_J}{|\,\boldsymbol{G}\,|^{2}}\mathrm{e}^{\mathrm{i}\boldsymbol{G}\cdot(\boldsymbol{r}_I-\boldsymbol{r}_J)}\mathrm{e}^{-\sigma^{2}|\boldsymbol{G}|^{2}/2} \tag{3.436} $$
将式（3.424）、式（3.426）、式（3.436）代入 $E_{\mathrm{es}}$ 的表达式（3.417），得到最后结果：
$$ \begin{aligned} E_{\mathrm{es}}={}&\frac{1}{2}\sum_{n}\sum_{I}\sum_{J}{}'\frac{Z_IZ_J}{|\,\boldsymbol{r}_I-\boldsymbol{r}_J+n\boldsymbol{L}\,|}\mathrm{erfc}\frac{|\,\boldsymbol{r}_I-\boldsymbol{r}_J+n\boldsymbol{L}\,|}{\sqrt{2}\sigma}\\ &+\frac{2\pi}{\varOmega_{\mathrm{cell}}}\sum_{\boldsymbol{G}\neq\boldsymbol{0}}\sum_{I}\sum_{J}\frac{Z_IZ_J}{|\,\boldsymbol{G}\,|^{2}}\mathrm{e}^{\mathrm{i}\boldsymbol{G}\cdot(\boldsymbol{r}_I-\boldsymbol{r}_J)}\mathrm{e}^{-\sigma^{2}|\boldsymbol{G}|^{2}/2}-\frac{1}{\sqrt{2\pi}\sigma}\sum_{I}Z_I^{2} \end{aligned} \tag{3.437} $$
实际应用中的总能表达式
由 3.4.3.2 节的结果，可以推导在实际应用中的总能 $E_{\mathrm{tot}}$ 的表达式。为了后面的计算方便，取附加正电荷分布 $\rho_{\mathrm{aux}}(\boldsymbol{r})$ 为
$$ \rho_{\mathrm{aux}}(\boldsymbol{r})=-\sum_{i}\frac{Z_i}{(\pi\sigma^{2})^{3/2}}\mathrm{e}^{-|\boldsymbol{r}-\boldsymbol{R}_i|^{2}/\sigma^{2}} \tag{3.438} $$
式中：$\boldsymbol{r}$ 遍历全空间；右端最前面有负号是因为采用的是原子单位制，$e=1$。相应的傅里叶变换为
$$ \rho_{\mathrm{aux}}(\boldsymbol{G})=-\frac{1}{\varOmega}\mathrm{e}^{-|\boldsymbol{G}|^{2}\sigma^{2}/4}\left[\sum_{s=1}^{P_s}Z_sS^{s}(\boldsymbol{G})\right] \tag{3.439} $$
而相互作用能 $E_{\mathrm{aux}}$ 为
$$ E_{\mathrm{aux}}=\frac{1}{2}\iint\frac{\rho_{\mathrm{aux}}(\boldsymbol{r})\rho_{\mathrm{aux}}(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}' \tag{3.440} $$
将式（3.440）加入总能表达式（3.378），再将其从中减去，则根据式（3.385）、式（3.393）、式（3.402）等，可以将方程（3.378）中的静电能部分表示为
$$ \begin{aligned} E_{\mathrm{es}}={}&E_{\mathrm{ee}}+E_{\mathrm{Ie}}^{\mathrm{loc}}+E_{\mathrm{II}}\\ ={}&\frac{1}{2}\iint\frac{\rho(\boldsymbol{r})\rho(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}'+\int\rho(\boldsymbol{r})\left(\sum_{n}\sum_{s=1}^{P_s}\sum_{I=1}^{N_s}V_{\mathrm{ps}}^{\mathrm{loc},s}\,|\,\boldsymbol{r}-\boldsymbol{R}_I+n\boldsymbol{L}\,|\right)\mathrm{d}\boldsymbol{r}\\ &+\frac{1}{2}\sum_{n}\sum_{I}\sum_{J}{}'\frac{Z_IZ_J}{|\,\boldsymbol{r}_I-\boldsymbol{r}_J+n\boldsymbol{L}\,|}-\frac{1}{2}\iint\frac{\rho_{\mathrm{aux}}(\boldsymbol{r})\rho_{\mathrm{aux}}(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}'\\ &+\frac{1}{2}\iint\frac{\rho_{\mathrm{aux}}(\boldsymbol{r})\rho_{\mathrm{aux}}(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}' \end{aligned} \tag{3.441} $$
将赝势的非局域部分单独处理。引入总电荷密度 $\rho_{\mathrm{T}}(\boldsymbol{r})=\rho(\boldsymbol{r})+\rho_{\mathrm{aux}}(\boldsymbol{r})$。显然，如果体系呈电中性，则有
$$ Q_{\mathrm{T}}=\int\rho_{\mathrm{T}}(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r}=0 \tag{3.442} $$
经过简单的计算，可将 $E_{\mathrm{es}}$ 重新表示为
$$ \begin{aligned} E_{\mathrm{es}}={}&\frac{1}{2}\iint\frac{\rho_{\mathrm{T}}(\boldsymbol{r})\rho_{\mathrm{T}}(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}'+\int\rho(\boldsymbol{r})\left(\sum_{n}\sum_{s=1}^{P_s}\sum_{I=1}^{N_s}V_{\mathrm{ps}}^{\mathrm{loc},s}\,|\,\boldsymbol{r}-\boldsymbol{R}_{I,s}+n\boldsymbol{L}\,|-\int\mathrm{d}\boldsymbol{r}'\frac{\rho_{\mathrm{aux}}(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\right)\mathrm{d}\boldsymbol{r}\\ &+\frac{1}{2}\left[\sum_{n}\sum_{I}\sum_{J}{}'\frac{Z_IZ_J}{|\,\boldsymbol{r}_I-\boldsymbol{r}_J+n\boldsymbol{L}\,|}-\iint\frac{\rho_{\mathrm{aux}}(\boldsymbol{r})\rho_{\mathrm{aux}}(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}'\right] \end{aligned} \tag{3.443} $$
根据 3.4.3.1 节中的讨论，在平面波基组的表象下，式（3.443）中，
$$ E_{\mathrm T}^{\mathrm H}=\frac{\varOmega}{2} \sum_{\boldsymbol G\ne\boldsymbol0}\frac{4\pi}{|\boldsymbol G|^2} \rho_{\mathrm T}(\boldsymbol G)\rho_{\mathrm T}(-\boldsymbol G) =2\pi\varOmega\sum_{\boldsymbol G\ne\boldsymbol0} \frac{|\rho_{\mathrm T}(\boldsymbol G)|^2}{|\boldsymbol G|^2} \tag{3.444} $$
其中 $\boldsymbol{G}=\boldsymbol{0}$ 的一项为发散项，但是 $\rho_{\mathrm{T}}(0)$ 等于 $Q_{\mathrm{T}}/\varOmega$，由式（3.442）可知该项为零，因此发散项消失。同理，有
$$ \iint\frac{\rho_{\mathrm{aux}}(\boldsymbol{r})\rho_{\mathrm{aux}}(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}'=4\pi\varOmega\sum_{\boldsymbol{G}\neq\boldsymbol{0}}\frac{\rho_{\mathrm{aux}}(\boldsymbol{G})\rho_{\mathrm{aux}}(-\boldsymbol{G})}{|\,\boldsymbol{G}\,|^{2}} \tag{3.445} $$
式（3.443）中 $\boldsymbol{G}=\boldsymbol{0}$ 的一项在后面单独处理，而包括 $V_{\mathrm{ps}}^{\mathrm{loc},s}$ 的项由方程（3.393）给出，有
$$ \int\rho(\boldsymbol{r})\left(\sum_{n}\sum_{s=1}^{P_s}\sum_{I=1}^{N_s}V_{\mathrm{ps}}^{\mathrm{loc},s}\,|\,\boldsymbol{r}-\boldsymbol{R}_I+n\boldsymbol{R}_{I,s}\,|\right)\mathrm{d}\boldsymbol{r}=\varOmega\sum_{|\boldsymbol{G}|}\sum_{s=1}^{N_s}S^{s}(\boldsymbol{G})V_{\mathrm{ps}}^{\mathrm{loc},s}(\boldsymbol{G})\rho(\boldsymbol{G}) \tag{3.446} $$
注意，与式（3.396）相同，此时第二次求和只在一个单胞内进行。因此，方程（3.443）中的 $E_{\mathrm{Ie}}^{\mathrm{loc}}$ 项为
$$ E_{\mathrm{Ie}}^{\mathrm{loc}} =\varOmega\sum_{\boldsymbol G}\left[ \sum_{s=1}^{P_s}S^s(\boldsymbol G)V_{\mathrm{ps}}^{\mathrm{loc},s}(\boldsymbol G) -\frac{4\pi\rho_{\mathrm{aux}}(\boldsymbol G)}{|\boldsymbol G|^2}\right]\rho(-\boldsymbol G), \quad \boldsymbol G=0\text{ 项按共同的零频约定取极限} \tag{3.447} $$
其中 $\boldsymbol{G}=\boldsymbol{0}$ 的一项要进行特殊的处理。
首先考虑方程（3.447）方括号中的第一项 $\displaystyle\sum_{s=1}^{N_s}S^{s}(\boldsymbol{G})V_{\mathrm{ps}}^{\mathrm{loc},s}(\boldsymbol{G})$。当 $\boldsymbol{G}=\boldsymbol{0}$ 时，由式（3.396）可知 $S^{s}(\boldsymbol{G})=N_s$。而 $V_{\mathrm{ps}}^{\mathrm{loc},s}(\boldsymbol{G})$ 可计算如下：
$$ \begin{aligned} V_{\mathrm{ps}}^{\mathrm{loc},s}(\boldsymbol{G})&=\frac{1}{\varOmega_{\mathrm{cell}}}\int_{\varOmega_{\mathrm{cell}}}V_{\mathrm{ps}}^{\mathrm{loc},s}\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r}\\ &=\frac{1}{\varOmega_{\mathrm{cell}}}\int_{|\boldsymbol{r}|\lt r_{\mathrm{c}}}\left(V_{\mathrm{ps}}^{\mathrm{loc},s}+\frac{Z_s}{r}\right)\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r}+\frac{1}{\varOmega_{\mathrm{cell}}}\int_{\varOmega_{\mathrm{cell}}}\left(-\frac{Z_s}{r}\right)\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r}\\ &=\frac{1}{\varOmega_{\mathrm{cell}}}\int_{|\boldsymbol{r}|\lt r_{\mathrm{c}}}\left(V_{\mathrm{ps}}^{\mathrm{loc},s}+\frac{Z_s}{r}\right)\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r}-\frac{4\pi Z_s}{\varOmega_{\mathrm{cell}}\,|\,\boldsymbol{G}\,|^{2}} \end{aligned} \tag{3.448} $$
式（3.448）使用了核吸引势在截断半径之外的渐近形式 $V_{\mathrm{ps}}^{\mathrm{loc},s}(r)=-Z_s/r$。考虑 $\boldsymbol{G}=\boldsymbol{0}$ 的项，式（3.448）中最后一行第一项记为 $\alpha^{s}$，即
$$ \alpha^{s}=\frac{1}{\varOmega_{\mathrm{cell}}}\int_{|\boldsymbol{r}|\lt r_{\mathrm{c}}}\left(V_{\mathrm{ps}}^{\mathrm{loc},s}+\frac{Z_s}{r}\right)\mathrm{d}\boldsymbol{r} \tag{3.449} $$
它的值是非零的有限值；第二项为发散项，暂且记为 $V_{\mathrm{c}}^{\mathrm{loc},s}(0)$。
其次考虑式（3.447）方括号中的第二项。因为 $\rho_{\mathrm{aux}}(\boldsymbol{G})$ 前面有因子 $4\pi/|\,\boldsymbol{G}\,|^{2}$，所以应该将 $\rho_{\mathrm{aux}}(\boldsymbol{G})$ 按照 $|\,\boldsymbol{G}\,|$ 展开到 $|\,\boldsymbol{G}\,|^{2}$ 项，再取 $|\,\boldsymbol{G}\,|\to0$ 的极限。由式（3.428）、式（3.438）可得
$$ \begin{aligned} \rho_{\mathrm{aux}}(\boldsymbol G) &=-\frac{e^{-\sigma^2|\boldsymbol G|^2/4}}{\varOmega_{\mathrm{cell}}} \sum_I Z_Ie^{-\mathrm i\boldsymbol G\cdot\boldsymbol R_I}\\ &=-\frac Q{\varOmega_{\mathrm{cell}}} +\frac{\mathrm i}{\varOmega_{\mathrm{cell}}}\boldsymbol G\cdot\sum_I Z_I\boldsymbol R_I +\frac{Q\sigma^2|\boldsymbol G|^2}{4\varOmega_{\mathrm{cell}}}\\ &\quad+\frac{1}{2\varOmega_{\mathrm{cell}}}\sum_I Z_I(\boldsymbol G\cdot\boldsymbol R_I)^2 +O(|\boldsymbol G|^3) \end{aligned} \tag{3.450} $$
式（3.450）表明一般晶胞还含位置相关的低阶矩；以下零频极限的简式仅在这些矩项被相应的 Ewald 项抵消或单独处理时适用。由式（3.449）及 $V_{\mathrm{c}}^{\mathrm{loc},s}(0)$ 可得，$E_{\mathrm{Ie}}^{\mathrm{loc}}$ 中 $|\,\boldsymbol{G}\,|=0$ 的一项（记为 $\bar{E}_{\mathrm{Ie}}^{\mathrm{loc}}$）为
$$ \bar E_{\mathrm{Ie}}^{\mathrm{loc}} =\varOmega\lim_{\boldsymbol G\to0}\left[ \sum_{s=1}^{P_s}S^s(\boldsymbol G)V_{\mathrm{ps}}^{\mathrm{loc},s}(\boldsymbol G) -\frac{4\pi\rho_{\mathrm{aux}}(\boldsymbol G)}{|\boldsymbol G|^2}\right]\rho(-\boldsymbol G) \tag{3.451} $$
由于辅助电荷与离子电荷相反且总量相等，两项的 $|\boldsymbol G|^{-2}$ 发散部分在式（3.451）的组合中抵消。有限部分依赖于所选零频边界条件。结合式（3.439）、式（3.447）及式（3.451），得到
$$ E_{\mathrm{Ie}}^{\mathrm{loc}} =\varOmega\sum_{\boldsymbol G\ne0}\left[ \sum_{s=1}^{P_s}S^s(\boldsymbol G)V_{\mathrm{ps}}^{\mathrm{loc},s}(\boldsymbol G) -\frac{4\pi\rho_{\mathrm{aux}}(\boldsymbol G)}{|\boldsymbol G|^2}\right]\rho(-\boldsymbol G) +\bar E_{\mathrm{Ie}}^{\mathrm{loc}} \tag{3.452} $$
方程（3.443）中的第二行记为 $E_{\mathrm{II}}^{\mathrm{mix}}$，其中方括号内第二项记为 $E_{\mathrm{aux}}$，其表达式在倒空间中可以写为与式（3.444）、式（3.445）类似的形式：
$$ E_{\mathrm{aux}}^{(\boldsymbol G\ne0)} =\frac\varOmega2\sum_{\boldsymbol G\ne0}\frac{4\pi}{|\boldsymbol G|^2} \rho_{\mathrm{aux}}(\boldsymbol G)\rho_{\mathrm{aux}}(-\boldsymbol G) =\frac\varOmega2\sum_{\boldsymbol G\ne0}\frac{4\pi}{|\boldsymbol G|^2} \left|\rho_{\mathrm{aux}}(\boldsymbol G)\right|^2 \tag{3.453} $$
根据式（3.453）、式（3.437），可得
$$ E_{\mathrm{II}}^{\mathrm{mix}}=E_{\mathrm{II}}-E_{\mathrm{aux}} \quad\text{（两项使用同一周期零频与背景约定）} \tag{3.454} $$
其中右端倒数第三、四项与最后两项分别为第二行方括号中第一项和第二项在 $\boldsymbol{G}\to\boldsymbol{0}$ 时展开到 $|\,\boldsymbol{G}\,|^{2}$ 项的极限（见式（3.450））。$\rho_{\mathrm{aux}}(0)$ 由式（3.438）给出，容易求得
$$ \rho_{\mathrm{aux}}(0)=-Q/\varOmega_{\mathrm{cell}},\quad\rho''_{\mathrm{aux}}(0)=Q\sigma^{2}/(2\varOmega_{\mathrm{cell}}) $$
对带电的各单项应统一规定 $\boldsymbol G=0$ 的背景与边界条件；只有在整体中性且采用同一约定后，静电能组合才是有限的。式（3.453）给出非零倒格矢的自能贡献，式（3.454）保留其余项与离子间能的精确代数关系。
合并各项后，静电能的等价恒等式为：
$$ \begin{aligned} E_{\mathrm{es}}={}&\frac12\iint \frac{\rho_{\mathrm T}(\boldsymbol r)\rho_{\mathrm T}(\boldsymbol r')}{|\boldsymbol r-\boldsymbol r'|} \,\mathrm d\boldsymbol r\,\mathrm d\boldsymbol r'\\ &+\int\rho(\boldsymbol r)\left[V_{\mathrm{ps}}^{\mathrm{loc}}(\boldsymbol r) -V_{\mathrm{aux}}(\boldsymbol r)\right]\mathrm d\boldsymbol r +E_{\mathrm{II}}-E_{\mathrm{aux}},\\ &V_{\mathrm{aux}}(\boldsymbol r)=\int \frac{\rho_{\mathrm{aux}}(\boldsymbol r')}{|\boldsymbol r-\boldsymbol r'|}\,\mathrm d\boldsymbol r'. \end{aligned} \tag{3.455} $$
辅助高斯电荷只用于拆分长程库仑项；物理总能不应依赖任意选定的高斯宽度 $\sigma$。
总能的其他部分，如动能、交换关联能、非局域赝势项等已经在 3.4.3.1 节中给出。从原则上讲，体系的总能可以在求解本征值的同时得到。但是此时的本征函数及电荷分布均未更新，所以与通常所说的 Kohn-Sham 总能有所不同。这说明，仍然需要首先得到体系的本征值和相应的本征方程，之后才能计算更新后体系的 Kohn-Sham 总能。
倒空间中总能表达式（3.394）、式（3.395）与式（3.455）看起来并不是很协调，最明显的区别在于式（3.394）和式（3.395）中用 $\sum V_{\mathrm{H}}(\boldsymbol{G})\rho(\boldsymbol{G})$ 项来表示电子相互作用，而式（3.455）则借助人为构建的 $\rho_{\mathrm{T}}(\boldsymbol{G})$ 进行描述。为使这三个公式达成一致，将 $\rho_{\mathrm{T}}(\boldsymbol{r})=\rho(\boldsymbol{r})+\rho_{\mathrm{aux}}(\boldsymbol{r})$ 代入式（3.455），同时定义 $\gamma_{\mathrm{Ewald}}$ 为
$$ \begin{aligned} \gamma_{\mathrm{Ewald}}={}&\frac12\sum_{n,I,J}{}'\frac{Z_IZ_J\,\operatorname{erfc}\!\left(|\boldsymbol R_I-\boldsymbol R_J+n\boldsymbol L|/(\sqrt2\sigma)\right)}{|\boldsymbol R_I-\boldsymbol R_J+n\boldsymbol L|}\\ &+\frac{2\pi}{\varOmega_{\mathrm{cell}}}\sum_{\boldsymbol G\ne0}\frac{e^{-\sigma^2|\boldsymbol G|^2/2}}{|\boldsymbol G|^2} \left|\sum_I Z_Ie^{-\mathrm i\boldsymbol G\cdot\boldsymbol R_I}\right|^2 -\frac{1}{\sqrt{2\pi}\sigma}\sum_I Z_I^2 -\frac{\pi\sigma^2}{\varOmega_{\mathrm{cell}}}\left(\sum_I Z_I\right)^2 \end{aligned} \tag{3.456} $$
经过简单的计算，就可以得到比较常见的平面波-赝势框架下的单胞总能表达式：
$$ \begin{aligned} E_{\mathrm{tot}}={}&\varOmega_{\mathrm{cell}}\left[\frac{1}{2N_k}\sum_{n,i,\boldsymbol{G}}|\,c_{n,\boldsymbol{k}_i}(\boldsymbol{G})\,|^{2}\,|\,\boldsymbol{k}_i+\boldsymbol{G}\,|^{2}+\sum_{\boldsymbol{G}\neq\boldsymbol{0}}\left(\frac{1}{2}V_{\mathrm{H}}(\boldsymbol{G})+\sum_{s=1}^{P_s}V_{\mathrm{ps}}^{\mathrm{loc},s}(\boldsymbol{G})S^{s}(\boldsymbol{G})\right)\rho(-\boldsymbol{G})\right.\\ &+\sum_{\boldsymbol{G}}\rho(\boldsymbol{G})\varepsilon_{\mathrm{xc}}(-\boldsymbol{G})+\frac{1}{N_k}\sum_{n,i,l,\boldsymbol{G},\boldsymbol{G}'}c_{n,\boldsymbol{k}_i}^{*}(\boldsymbol{G})c_{n,\boldsymbol{k}_i}(\boldsymbol{G}')\sum_{s=1}^{P_s}\delta V_{\mathrm{ps},l,\boldsymbol{k}_i+\boldsymbol{G},\boldsymbol{k}_i+\boldsymbol{G}'}^{s}\\ &\left.\vphantom{\frac11}\right]+\gamma_{\mathrm{Ewald}}+\sum_{s=1}^{P_s}(N_s\alpha^{s}Z_s) \end{aligned} \tag{3.457} $$
利用 3.4.1 节中介绍的特殊 $\boldsymbol{k}$ 点法，将 $\boldsymbol{k}$ 的取值限制在第一布里渊区的不可约区域内，则可得
$$ \begin{aligned} E_{\mathrm{tot}}={}&\frac{\varOmega_{\mathrm{cell}}}{N_k}\sum_{i}\omega_{\boldsymbol{k}_i}\left[\sum_{n,\boldsymbol{G},\boldsymbol{G}'}c_{n,\boldsymbol{k}_i}^{*}(\boldsymbol{G})\left(\frac{1}{2}\,|\,\boldsymbol{k}_i+\boldsymbol{G}\,|^{2}\delta_{\boldsymbol{G},\boldsymbol{G}'}\right.\right.\\ &\left.\left.+\sum_{l}\sum_{s=1}^{P_s}\delta V_{\mathrm{ps},l,\boldsymbol{k}_i+\boldsymbol{G},\boldsymbol{k}_i+\boldsymbol{G}'}^{s}\right)c_{n,\boldsymbol{k}_i}(\boldsymbol{G}')\right]+\varOmega_{\mathrm{cell}}\sum_{\boldsymbol{G}}\rho(\boldsymbol{G})\varepsilon_{\mathrm{xc}}(-\boldsymbol{G})\\ &+\varOmega_{\mathrm{cell}}\sum_{\boldsymbol{G}\neq\boldsymbol{0}}\left(\frac{1}{2}V_{\mathrm{H}}(\boldsymbol{G})+\sum_{s=1}^{P_s}V_{\mathrm{ps}}^{\mathrm{loc},s}(\boldsymbol{G})S^{s}(\boldsymbol{G})\right)\rho(-\boldsymbol{G})+\gamma_{\mathrm{Ewald}}+\sum_{s=1}^{P_s}(N_s\alpha^{s}Z_s) \end{aligned} \tag{3.458} $$
式中：$N_k$ 是第一布里渊区内 $\boldsymbol{k}$ 点的数目。
为了计算体系本征值，需要对平面波-赝势框架下的 Kohn-Sham 方程（式（3.400））中的势能项 $V_{\boldsymbol{G},\boldsymbol{G}'}^{i}$ 加以限制。令 $V_{\mathrm{H}}(0)$ 及 $V_{\mathrm{ps}}^{\mathrm{loc}}(0)$ 等于 0，因此，当 $\boldsymbol{G}=\boldsymbol{G}'$ 时，有
$$ V_{\boldsymbol{G},\boldsymbol{G}}^{i}=\mu_{\mathrm{xc}}(0)+\sum_{s=1}^{P_s}N_s\sum_{l}\delta V_{\mathrm{ps},l,\boldsymbol{k}_i+\boldsymbol{G},\boldsymbol{k}_i+\boldsymbol{G}}^{\mathrm{nl},s} \tag{3.459} $$
这种直接忽略发散项的做法相当于平移了势能零点，因此需要对最后的总能表达式做出修正\cite{ihm1979momentum}。而由式（3.458）可知，修正由最后的 $\alpha Z$ 项给出。
有必要指出，实际计算中，即使在平面波基组下，也并不是所有能量项都适合在动量空间中求解。式（3.457）和式（3.458）中的交换关联项，在考虑更复杂的 $\varepsilon_{\mathrm{xc}}$ 函数形式时（例如在 GGA 中），可能无法表示成如此简单的形式。更适合的方法是在实空间中计算：
$$ E_{\mathrm{xc}}=\frac{\varOmega}{N_r}\sum_{i=1}^{N_r}\varepsilon_{\mathrm{xc}}[\rho(\boldsymbol{r}_i)]\rho(\boldsymbol{r}_i) \tag{3.460} $$
此外，赝势非局域部分对总能的贡献 $E_{\mathrm{ps}}^{\mathrm{nl}}$ 也可以通过 KB 分部形式重新写出：
$$ E_{\mathrm{ps}}^{\mathrm{nl}}=\frac{\varOmega_{\mathrm{cell}}}{N_k}\sum_{i}\omega_{\boldsymbol{k}_i}\sum_{l=0}^{l_{\max}}\sum_{m=-l}^{l}\sum_{s=1}^{P_s}\sum_{I=1}^{N_s}\sum_{n}^{N_{\mathrm{stat}}}\beta_{lm}^{s}\,|\,F_{I,n}^{lm,s}(\boldsymbol{k}_i+\boldsymbol{G})\,|^{2} \tag{3.461} $$
其中 $\beta_{lm}^{s}$ 由式（3.398）给出，而 $F_{I,n}^{lm,s}$ 表示为
$$ F_{I,n}^{lm,s}(\boldsymbol{k}_i+\boldsymbol{G})=\sum_{\boldsymbol{G}}\mathrm{e}^{\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{R}_{I,s}}f_{lm}^{s}(\boldsymbol{k}_i+\boldsymbol{G})c_{n,\boldsymbol{k}_i}(\boldsymbol{G}) \tag{3.462} $$