本文是「第一性原理的微观计算模拟」系列的第 4 篇（共 11 篇），内容整理自同名书稿第 3 章，公式编号与原书一致。文中引用的文献依据原书参考文献清单整理并列于文末，按原书清单顺序从 1 开始编号。\nocite{*}
正交化平面波
利用平面波函数作为基函数有很多优点，但是其中有一个很显著的缺陷，即原子的内层电子波函数在靠近原子核的区域有很大的振荡，需要用数目很大的平面波基组展开这些波函数才能获得比较精确的结果。这种处理方法无疑极大地增加了计算量。考虑到原子的内层电子并不参与成键，在固体或分子体系中芯区轨道与自由原子状态相比几乎不变，而外层电子，或者更精确地说，处于价带或者导带中的电子才是我们的研究重点，因此可以将这两种轨道分开处理。
首先介绍正交化平面波（orthogonalized plane wave，OPW）方法\cite{phillips1959method}。
价带或者导带的 Bloch 波函数应与芯区轨道的波函数正交。如果已知各芯区轨道的波函数 $\varphi_j$（满足薛定谔方程 $H\varphi_j=\varepsilon_j\varphi_j$），则可以构建满足上述正交条件的波函数的普遍表达式：
$$ \psi_n(\boldsymbol{k},\boldsymbol{r})=\chi_n(\boldsymbol{k},\boldsymbol{r})-\sum_{j}\langle\varphi_j\,|\,\chi_n\rangle\varphi_j \tag{3.236} $$
不难验证，$\langle\psi_n\,|\,\varphi_j\rangle=0$。式（3.236）中并没有对 $\chi_n$ 予以明确定义，如果取其为平面波，则式（3.236）所表示的波称为正交化平面波。考虑一个孤立原子，将式（3.236）代入该原子的薛定谔方程，有
$$ \hat{H}\chi_n+\sum_{j}(\varepsilon_n-\varepsilon_j)\,|\,\varphi_j\rangle\langle\varphi_j\,|\,\chi_n=\varepsilon_n\chi_n \tag{3.237} $$
与原方程比较，式（3.237）中的哈密顿量多了一项 $\displaystyle\sum_{j}(\varepsilon_n-\varepsilon_j)\,|\,\varphi_j\rangle\langle\varphi_j\,|$。原子的哈密顿量表示为动能算符与库仑势之和，即 $\hat{H}=\hat{T}+\hat{V}$，其中 $\hat{V}=-(Z/r)\boldsymbol{I}$，也即裸核的库仑势，$\boldsymbol{I}$ 是单位矩阵。因此式（3.237）表示 $\chi_n$ 满足下列方程
$$ (\hat{T}+\hat{V}_{\mathrm{PK}})\chi_n=\varepsilon_n\chi_n \tag{3.238} $$
式中
$$ \hat{V}_{\mathrm{PK}}=-\frac{Z}{r}\boldsymbol{I}+\sum_{j}(\varepsilon_n-\varepsilon_j)\,|\,\varphi_j\rangle\langle\varphi_j\,| \tag{3.239} $$
这意味着可以将内层电子视为一个等效屏蔽势函数，而不对其进行精确的求解。这也是赝势（pseudopotential）理论最早的由来\cite{phillips1959method}。因为 $\varepsilon_n-\varepsilon_j$ 恒大于零，所以 $\hat{V}_{\mathrm{PK}}$ 比 $\hat{V}$ 弱，其所对应的波函数 $\chi_n$ 在原子核附近更为平滑，随着 $\boldsymbol{G}$ 的增大，其傅里叶变换相应的分量也减小得更快。
由正交化平面波方法的思想出发，我们可以构建每种原子的等效势算符 $\hat{V}_{\mathrm{ps}}(r)$，称为赝势。相应的薛定谔方程的解 $\psi_{\mathrm{ps}}$ 称为赝波函数。由式（3.239）可知，$\hat{V}_{\mathrm{ps}}$ 应与轨道的角动量量子数 $L=\{l,m\}$ 相关。为使计算简便，我们可以进一步设定对于每个角动量量子数为 $L$ 的轨道，$\hat{V}_{\mathrm{ps}}(r)$ 是一个球对称的函数。综合以上讨论，比照式（3.239），可将赝势写为
$$ \hat{V}_{\mathrm{ps}}(r)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}V_{\mathrm{ps}}^{l}(r)\,|\,lm\rangle\langle lm\,|=\sum_{l=0}^{\infty}V_{\mathrm{ps}}^{l}(r)\hat{P}_l \tag{3.240} $$
坐标表象下 $|\,lm\rangle$ 为球谐函数 $\mathrm{Y}_l^{m}(\theta,\phi)$，而投影算符 $\hat{P}_l$ 为
$$ \hat{P}_l=\sum_{m=-l}^{l}|\,lm\rangle\langle lm\,| \tag{3.241} $$
因此，将 $\hat{V}_{\mathrm{ps}}(r)$ 作用于波函数时，首先将其投影到不同的 $l$ 分量上，并对该分量作用相应的 $V_{\mathrm{ps}}^{l}$，最后对各分量的结果求和。方程（3.241）表明，赝势算符同时包含 $|\,lm\rangle$ 的变量 $(\theta,\phi)$ 及 $\langle lm\,|$ 的变量 $(\theta',\phi')$，但是径向部分仅与 $r$ 一个变量有关。因此，$\hat{V}_{\mathrm{ps}}$ 是一个非局域（non-local）算符，更准确地说，是一个半局域（semi-local）算符（角向部分为非局域的，径向部分为局域的），这一点可以表示如下：将 $\hat{V}_{\mathrm{ps}}(r)$ 作用于某函数 $f(r,\theta',\phi')$，可得
$$ (\hat V_{\mathrm{ps}}f)(r,\theta,\phi) =\sum_{l=0}^{\infty}\sum_{m=-l}^{l}Y_l^m(\theta,\phi)V_{\mathrm{ps}}^l(r) \int_0^{2\pi}\!\mathrm d\phi'\int_0^\pi\!\sin\theta'\,\mathrm d\theta'\, Y_l^{m*}(\theta',\phi')f(r,\theta',\phi') \tag{3.242} $$
模守恒赝势
在 3.3.1 节中，我们推导出了赝势的普遍表达式，但是尚未讨论如何得到 $V_{\mathrm{ps}}^{l}(r)$ 以及构建赝势时应该满足的条件或者性质。在本节以及随后的几节中我们将对此进行详细探讨。实际上，从 3.3.1 节的讨论中已经可以得出某些结论，例如：赝波函数与全电子波函数 $\varPsi_{\mathrm{ae}}$（或称真实波函数）拥有相同的本征能级 $\varepsilon_l$；赝势的建立需要给定的参考态，对同种原子如果选取其不同的电子组态，因为 $\varepsilon_j$ 不同，构造出来的赝势也会有差异。Hamann、Schlüter 和 Chiang 最早提出模守恒赝势（norm-conserving pseudopotential，NCPP）的概念以及一个优良的赝势应该满足的条件\cite{hamann1979norm}：
（1）赝波函数与作为其参考态的全电子波函数拥有相同的本征能级，即
$$ \tilde{\varepsilon}_l=\varepsilon_l \tag{3.243} $$
（2）在芯区截断半径 $r_{\mathrm{c}}$ 之外，赝波函数与全电子波函数完全重合，即
$$ \varPsi_{\mathrm{ps}}^{l}(r)=\varPsi_{\mathrm{ae}}^{l}(r),\quad r\gt r_{\mathrm{c}} \tag{3.244} $$
（3）在 $r_{\mathrm{c}}$ 终点处，赝波函数与全电子波函数的对数导数相等，即
$$ \frac{\mathrm{d}}{\mathrm{d}r}\ln\varPsi_{\mathrm{ps}}^{l}(r)\,|_{r=r_{\mathrm{c}}}=\frac{\mathrm{d}}{\mathrm{d}r}\ln\varPsi_{\mathrm{ae}}^{l}(r)\,|_{r=r_{\mathrm{c}}} \tag{3.245} $$
波函数的对数导数记为 $D_{\mathrm{ps}}^{l}$ 及 $D_{\mathrm{ae}}^{l}$。
（4）在 $r_{\mathrm{c}}$ 之内，赝波函数与全电子波函数对体积的积分相等，即
$$ \int_0^{r_{\mathrm c}}r^2|R_{\mathrm{ps}}^l(r)|^2\,\mathrm dr =\int_0^{r_{\mathrm c}}r^2|R_{\mathrm{ae}}^l(r)|^2\,\mathrm dr \tag{3.246} $$
（5）在 $r_{\mathrm{c}}$ 终点处，赝波函数与全电子波函数的对数导数相对于能量的一阶导数相等，即
$$ \frac{\partial D_{\mathrm{ps}}^{l}}{\partial\varepsilon}=\frac{\partial D_{\mathrm{ae}}^{l}}{\partial\varepsilon} $$
条件（1）、（2）表明，引入赝势不应对元素在芯区之外的电子结构产生干扰。条件（3）要求赝波函数与全电子波函数在截断半径处光滑连续。条件（4）即“模守恒条件”，它表明赝波函数在芯区内的电荷量正确。由于芯区外的势函数取决于芯区内的电荷总量，因此符合“模守恒条件”的赝势保证了在多原子体系内对原子间相互作用的描述是准确的。条件（5）与赝势的移植性有关。
通常，赝势是根据原子在孤立环境下的电子结构构建的，将其应用到相互作用的多原子体系中时，本征波函数及本征能级都会发生变化。如果赝势满足条件（5），则它同样可以反映这种变化，且在线性项上是正确的。上述讨论也可以根据散射理论进行：环境变化会导致全电子波函数的相移（phase shift）$\delta\eta_{\mathrm{ae}}^{l}$，而用赝势生成的赝波函数在相同的环境变化下也会产生相移 $\delta\eta_{\mathrm{ps}}^{l}$。$\delta\eta$ 是本征能级 $\varepsilon^{l}$ 的函数，将 $\delta\eta$ 关于 $\varepsilon^{l}$ 展开，则 $\delta\eta_{\mathrm{ae}}^{l}$ 与 $\delta\eta_{\mathrm{ps}}^{l}$ 的线性项相同。从上面的讨论中也可以看出，最后两点的联系非常紧密（均与本征能级 $\varepsilon^{l}$ 的变化有关）。更严格的数学推导指出，一种赝势若满足条件（4），则必然满足条件（5）。详细过程请参看文献\cite{martin2004electronic,hamann1979norm,shirley1989extended}。
构建赝势的普遍过程
对于一个体系，如果确定了体系的势函数，可以唯一地求解对应的波函数。而如果预知了体系的本征波函数，则可以反推相应的势函数。因此，构建赝势实际上是求解薛定谔方程的反问题。首先预设一个合适的赝波函数，注意满足前述的几个条件，然后通过反解薛定谔方程得到体系的赝势。一般设赝势具有球对称性，因此，薛定谔方程可以分离变量，角向和径向函数可以分别求解，即 $\varPsi_{\mathrm{ps}}^{l}(\boldsymbol{r})=R_{\mathrm{ps}}^{l}(r)\mathrm{Y}_{lm}(\theta,\phi)$，而 $R_{\mathrm{ps}}^{l}(r)$ 满足径向薛定谔方程（在原子单位制下，参见 2.1.9 节）
$$ \left[-\frac{1}{2}\frac{\mathrm{d}^{2}}{\mathrm{d}r^{2}}+\frac{l(l+1)}{2r^{2}}+V_{\mathrm{ps,scr}}^{l}(r)\right]rR_{\mathrm{ps}}^{l}(r)=\varepsilon_lrR_{\mathrm{ps}}^{l}(r) \tag{3.247} $$
由方程（3.247）立即可以解得屏蔽有效势 $V_{\mathrm{ps,scr}}^{l}(r)$，即
$$ V_{\mathrm{ps,scr}}^{l}(r)=\varepsilon_l-\frac{l(l+1)}{2r^{2}}+\frac{1}{2rR_{\mathrm{ps}}^{l}(r)}\frac{\mathrm{d}^{2}}{\mathrm{d}r^{2}}(rR_{\mathrm{ps}}^{l}(r)) \tag{3.248} $$
需要注意的是，屏蔽有效势 $V_{\mathrm{ps,scr}}^{l}(r)$ 并不是所求的原子赝势，因为其包含多体效应。因此需要再进行如下“去屏蔽”的操作：
$$ V_{\mathrm{ps}}^{l}(r)=V_{\mathrm{ps,scr}}^{l}(r)-\int\frac{\rho_{\mathrm{v}}(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}'-\mu_{\mathrm{xc}}[\rho_{\mathrm{v}}(\boldsymbol{r})] \tag{3.249} $$
式中，三维价电子数密度应由各占据赝轨道的模平方构成；径向概率分布另含球坐标体积因子，即
$$ \rho_{\mathrm v}(\boldsymbol r)=\sum_{nlm\sigma}f_{nlm\sigma} \left|R_{nl}(r)Y_{lm}(\hat{\boldsymbol r})\right|^2,\qquad \bar\rho_{\mathrm v}(r)=\frac1{4\pi}\int\rho_{\mathrm v}(r,\Omega)\,\mathrm d\Omega \tag{3.250} $$
这里径向分布中的 $r^2$ 来自体积元 $r^2\sin\theta\,\mathrm dr\,\mathrm d\theta\,\mathrm d\phi$，不能直接放进每单位体积的密度。方程（3.249）右端的第二项与第三项分别为 Hartree 势以及交换关联势，对此我们将在第 4 章详细讨论。
Troullier-Martins 赝势
Troullier 与 Martins 于 1991 年提出了一种模守恒赝势——TM 赝势\cite{troullier1991efficient}，它也是其后提出的很多模守恒赝势的模板。TM 赝势是对 Kerker 早期工作\cite{kerker1980non}的扩展。TM 赝势将赝波函数表示为
$$ R_{\mathrm{ps}}^{l}(r)=\begin{cases}R_{\mathrm{ae}}^{l}(r),&r\gt r_{\mathrm{c}}\\r^{l}\exp[p(r)],&r\leqslant r_{\mathrm{c}}\end{cases} \tag{3.251} $$
式中：$p(r)$ 是一个多项式，有
$$ p(r)=c_0+c_2r^{2}+c_4r^{4}+c_6r^{6}+c_8r^{8}+c_{10}r^{10}+c_{12}r^{12} \tag{3.252} $$
其中 $c_0,c_2,\cdots,c_{12}$ 是七个待定系数。根据方程（3.248），可得 $V_{\mathrm{ps,scr}}^{l}(r)$ 满足
$$ V_{\mathrm{ps,scr}}^{l}(r)=\begin{cases}V_{\mathrm{ae}}^{l}(r),&r\gt r_{\mathrm{c}}\\\varepsilon_l+\dfrac{l+1}{r}p'(r)+\dfrac{1}{2}p''(r)+\dfrac{1}{2}[p'(r)]^{2},&r\leqslant r_{\mathrm{c}}\end{cases} \tag{3.253} $$
可以看到，式（3.252）中所有奇数项的系数均为零，因为 Troullier 与 Martins 发现这种设定可以使所生成的赝势随倒格矢 $\boldsymbol{G}$ 的增加而更快地趋于零，即改善赝势的光滑性。此外，引入限制条件
$$ \left.\frac{\mathrm d^2V_{\mathrm{ps,scr}}^l(r)}{\mathrm dr^2}\right|_{r=0}=0 \tag{3.254} $$
也可以有效地提升赝势的光滑性。$[V]''$ 代表 $V$ 对于 $r$ 的二阶导数。因此，$p(r)$ 中的七个待定系数可以由如下几个条件确定\cite{troullier1991efficient,giannozzi2004notes}。
（1）由模守恒条件，有
$$ 2c_0+\ln\left[\int_{0}^{r_{\mathrm{c}}}r^{2(l+1)}\exp(2p(r)-2c_0)\,\mathrm{d}r\right]=\ln\left(\int_{0}^{r_{\mathrm{c}}}r^{2}\,|\,R_{\mathrm{ae}}^{l}(r)\,|^{2}\,\mathrm{d}r\right) \tag{3.255} $$
（2）在 $r_{\mathrm{c}}$ 终点处 $rR_{\mathrm{ps}}^{l}(r)$ 与 $rR_{\mathrm{ae}}^{l}(r)$ 连续，即
$$ p(r_{\mathrm{c}})=\ln\left(\frac{P(r_{\mathrm{c}})}{r_{\mathrm{c}}^{l+1}}\right) \tag{3.256} $$
式中 $P(r_{\mathrm{c}})=rR_{\mathrm{ae}}^{l}(r)$，下同。
（3）在 $r_{\mathrm{c}}$ 终点处 $rR_{\mathrm{ps}}^{l}(r)$ 与 $rR_{\mathrm{ae}}^{l}(r)$ 对 $r$ 的一阶导数连续，即
$$ \left.\frac{\mathrm{d}(rR_{\mathrm{ps}}^{l})}{\mathrm{d}r}\right|_{r=r_{\mathrm{c}}}=\left.\frac{\mathrm{d}(rR_{\mathrm{ae}}^{l})}{\mathrm{d}r}\right|_{r=r_{\mathrm{c}}} \tag{3.257} $$
由此并利用条件（2）中的连续性条件，推出
$$ p'(r_{\mathrm c})=\frac{P'(r_{\mathrm c})}{P(r_{\mathrm c})}-\frac{l+1}{r_{\mathrm c}} \tag{3.258} $$
（4）在 $r_{\mathrm{c}}$ 终点处 $rR_{\mathrm{ps}}^{l}(r)$ 与 $rR_{\mathrm{ae}}^{l}(r)$ 对 $r$ 的二阶导数连续，并利用方程（3.248），有
$$ p''(r_{\mathrm{c}})=2V_{\mathrm{ae}}-2\varepsilon_l-\frac{2(l+1)}{r_{\mathrm{c}}}p'(r_{\mathrm{c}})-[p'(r_{\mathrm{c}})]^{2} \tag{3.259} $$
（5）在 $r_{\mathrm{c}}$ 终点处 $rR_{\mathrm{ps}}^{l}(r)$ 与 $rR_{\mathrm{ae}}^{l}(r)$ 对 $r$ 的三阶导数连续，直接对式（3.259）求导，得
$$ p^{(3)}(r_{\mathrm{c}})=2V'_{\mathrm{ae}}(r_{\mathrm{c}})+\frac{2(l+1)}{r_{\mathrm{c}}^{2}}p'(r_{\mathrm{c}})-\frac{2(l+1)}{r_{\mathrm{c}}}p''(r_{\mathrm{c}})-2p'(r_{\mathrm{c}})p''(r_{\mathrm{c}}) \tag{3.260} $$
（6）在 $r_{\mathrm{c}}$ 终点处 $rR_{\mathrm{ps}}^{l}(r)$ 与 $rR_{\mathrm{ae}}^{l}(r)$ 对 $r$ 的四阶导数连续，直接对式（3.260）求导，得
$$ \begin{aligned} p^{(4)}(r_{\mathrm{c}})={}&2V''_{\mathrm{ae}}(r_{\mathrm{c}})-\frac{4(l+1)}{r_{\mathrm{c}}^{3}}p'(r_{\mathrm{c}})+\frac{4(l+1)}{r_{\mathrm{c}}^{2}}p''(r_{\mathrm{c}})-\frac{2(l+1)}{r_{\mathrm{c}}}p^{(3)}(r_{\mathrm{c}})\\ &-2[p''(r_{\mathrm{c}})]^{2}-2p'(r_{\mathrm{c}})p^{(3)}(r_{\mathrm{c}}) \end{aligned} \tag{3.261} $$
（7）根据方程（3.253）以及方程（3.254），可得
$$ c_2^{2}+c_4(2l+5)=0 \tag{3.262} $$
全电子波函数以及全电子势对 $r$ 的导数可利用有限差分得到。由此，可以确定赝波函数，再根据方程（3.253）求得 $V_{\mathrm{ps,scr}}^{l}(r)$。如果已知交换关联势 $\mu_{\mathrm{xc}}[\rho_{\mathrm{v}}(r)]$，则根据式（3.249）可确定赝势 $V_{\mathrm{ps}}^{l}(r)$。
自旋-轨道耦合的处理
如果考虑自旋-轨道耦合，那么好量子数将是 $j=l\pm1$。首先生成 $j=l+1/2$ 的赝势 $V_{\mathrm{ps}}^{l+1/2}$ 以及 $j=l-1/2$ 的赝势 $V_{\mathrm{ps}}^{l-1/2}$。由此得
$$ V_{\mathrm{ps}}^{l}=\frac{1}{2l+1}[(l+1)V_{\mathrm{ps}}^{l+1/2}+lV_{\mathrm{ps}}^{l-1/2}] \tag{3.263} $$
$$ \delta V_{\mathrm{so}}^{l}=\frac{2}{2l+1}(V_{\mathrm{ps}}^{l+1/2}-V_{\mathrm{ps}}^{l-1/2}) \tag{3.264} $$
这样，赝势方程（3.242）可以表示为
$$ V_{\mathrm{ps}}^{l}=\sum_{lm}[\,|\,\mathrm{Y}_{lm}\rangle V_{\mathrm{ps}}^{l}\langle\mathrm{Y}_{lm}\,|+|\,\mathrm{Y}_{lm}\rangle\delta V_{\mathrm{so}}^{l}\boldsymbol{L}\cdot\boldsymbol{S}\langle\mathrm{Y}_{lm}\,|\,] \tag{3.265} $$
式中：$\boldsymbol{L}$ 为轨道角动量；$\boldsymbol{S}$ 为自旋角动量。更具体的讨论请参看文献\cite{bachelet1982pseudopotentials}。
赝势的分部形式
局域赝势形式
从原则上讲，对所有的 $l$ 轨道都应该单独建立赝势 $V_{\mathrm{ps}}^{l}(r)$。但是在实际应用中，仅对少数几个 $l\lt l_{\max}$ 的轨道分别建立赝势，而对 $l\gt l_{\max}$ 的轨道则认为感受到的赝势相同，也即其赝势与 $l$ 无关。与 $l$ 无关的 $l\gt l_{\max}$ 的轨道的赝势称为局域赝势 $V_{\mathrm{ps}}^{\mathrm{loc}}$。因此，可以重新将方程（3.242）写为
$$ \hat{V}_{\mathrm{ps}}(r)=\sum_{l=0}^{\infty}V_{\mathrm{ps}}^{\mathrm{loc}}(r)\hat{P}_l+\sum_{l=0}^{l_{\max}}(V_{\mathrm{ps}}^{l}(r)-V_{\mathrm{ps}}^{\mathrm{loc}})\hat{P}_l=V_{\mathrm{ps}}^{\mathrm{loc}}(r)\boldsymbol{I}+\sum_{l=0}^{l_{\max}}\delta V_{\mathrm{ps}}^{l}(r)\hat{P}_l \tag{3.266} $$
式中：$\delta V_{\mathrm{ps}}^{l}(r)$ 为短程函数，仅局限于芯区范围内，而其局域部分 $V_{\mathrm{ps}}^{\mathrm{loc}}(r)=V_{\mathrm{ps}}^{l_{\max}+1}(r)$，其中 $l_{\max}$ 一般选取芯区电子占据态的角动量量子数最大值，或者单质状态下最高占据轨道的角动量量子数。局域通道可以在满足芯区外正确渐近行为的候选势中选择，还应检查高角动量散射性质及分离式赝势可能出现的幽灵态。如果给定一组标准正交基 $\{\varphi_i\}$，则由方程（3.266）可计算矩阵元：
$$ \begin{aligned} V_{\mathrm{ps},i,j}&=\langle\varphi_i\,|\,\hat{V}_{\mathrm{ps}}(\boldsymbol{r})\,|\,\varphi_j\rangle=\langle\varphi_i\,|\,V_{\mathrm{ps}}^{\mathrm{loc}}(r)\boldsymbol{I}+\sum_{l=0}^{l_{\max}}\delta V_{\mathrm{ps}}^{l}(r)\hat{P}_l\,|\,\varphi_j\rangle\\ &=\langle\varphi_i|V_{\mathrm{ps}}^{\mathrm{loc}}|\varphi_j\rangle+\sum_{l=0}^{l_{\max}}\delta V_{\mathrm{ps}}^{l}(i,j) \end{aligned} \tag{3.267} $$
其中局域部分在一般基组中也可能有非对角矩阵元；半局域角动量修正的矩阵元为
$$ \delta V_{\mathrm{ps}}^{l}(i,j)= \sum_{m=-l}^{l}\int_0^\infty\!r^2\,\mathrm dr\;\delta V_{\mathrm{ps}}^l(r) \left[\int\!\mathrm d\Omega\,\varphi_i^*(r,\Omega)Y_{lm}(\Omega)\right] \left[\int\!\mathrm d\Omega'\,Y_{lm}^*(\Omega')\varphi_j(r,\Omega')\right]. $$
$$ \langle\varphi_i|\delta\hat V_{\mathrm{ps}}|\varphi_j\rangle =\sum_{l=0}^{l_{\max}}\sum_{m=-l}^{l}\int_0^\infty r^2\,\mathrm dr\; \delta V_{\mathrm{ps}}^l(r) \langle\varphi_i(r,\cdot)|Y_{lm}\rangle_\Omega \langle Y_{lm}|\varphi_j(r,\cdot)\rangle_\Omega \tag{3.268} $$
式（3.268）的第二行用到了方程（3.242）（将其中的 $V_{\mathrm{ps}}^{l}(r)$ 替换为 $\delta V_{\mathrm{ps}}^{l}(r)$），即径向的积分仅对 $r$ 进行。
两种常用的基函数——平面波函数和原子轨道波函数，均可表示为 $R(r)\mathrm{Y}_{lm}(\theta,\phi)$，所以相应地式（3.268）的积分可以表示成角向积分以及径向积分的乘积。其中角向部分由球谐函数的卷积决定，而径向部分则为 $\displaystyle\int R_i^{*}(r)\delta V_{\mathrm{ps}}^{l}(r)R_j(r)r^{2}\,\mathrm{d}r$。
对于平面波，$R_j(r)$ 为 $j$ 阶的球贝塞尔函数，而对于原子轨道波函数，$R_j(r)$ 可取类氢原子的径向波函数。
方程（3.268）径向部分的计算量相当大，因为需要对每一对 $\varphi_i$ 和 $\varphi_j$ 进行积分。以平面波为例，如果有 $N$ 个基函数，在第一布里渊区内有 $M$ 个 $\boldsymbol{k}$ 采样点，则对于每一个 $l$，需要进行 $MN^{2}/2$ 次计算。为了克服这个困难，Kleinman 和 Bylander 提出将半局域的 $\delta V_{\mathrm{ps}}^{l}(r)$ 表示成非局域的投影算符，从而达到减小计算量的目的。这就是著名的 Kleinman-Bylander（KB）非局域赝势形式。
KB 非局域赝势形式
1982 年，Kleinman 和 Bylander 提出了一种普适的方法\cite{kleinman1982efficacious}——将半局域的赝势变换为非局域的形式，也即式（3.268）可以写为如下形式
$$ \delta V_{\mathrm{NL}}(\boldsymbol r,\boldsymbol r') =\sum_a F_a(\boldsymbol r)G_a^*(\boldsymbol r') \tag{3.269} $$
其中 $F_i$ 与 $G_i$ 分别依赖于 $\boldsymbol{r}$ 和 $\boldsymbol{r}'$，同时需要满足一个条件，即当作用在赝波函数 $\varPsi_{\mathrm{ps}}^{lm}$ 上时，$\delta V_{\mathrm{NL}}(\boldsymbol{r},\boldsymbol{r}')$ 与 $\delta V_{\mathrm{ps}}^{l}(r)$ 的结果相同（其中 $\delta V_{\mathrm{ps}}^{l}(r)$ 由 $\varPsi_{\mathrm{ps}}^{lm}$ 生成）。为此 Kleinman 和 Bylander 构建了非局域算符
$$ \delta V_{\mathrm{NL}}=\sum_{l=0}^{l_{\max}}\sum_{m=-l}^{l}\frac{|\,\delta V_{\mathrm{ps}}^{l}\varPsi_{\mathrm{ps}}^{lm}\rangle\langle\delta V_{\mathrm{ps}}^{l}\varPsi_{\mathrm{ps}}^{lm}\,|}{\langle\varPsi_{\mathrm{ps}}^{lm}\,|\,\delta V_{\mathrm{ps}}^{l}\,|\,\varPsi_{\mathrm{ps}}^{lm}\rangle} \tag{3.270} $$
不难证明，利用式（3.270）构建的非局域赝势算符满足条件
$$ \delta\hat V_{\mathrm{NL}}|\Psi_{\mathrm{ps}}^{l'm'}\rangle =\delta V_{\mathrm{ps}}^{l'}(r)|\Psi_{\mathrm{ps}}^{l'm'}\rangle \tag{3.271} $$
KB 非局域赝势形式的优势在于将其作用在两个基函数上时，可写为
$$ \langle\varphi_i\,|\,\delta V_{\mathrm{NL}}\,|\,\varphi_j\rangle=\sum_{l=0}^{l_{\max}}\sum_{m=-l}^{l}\langle\varphi_i\,|\,\varPsi_{\mathrm{ps}}^{lm}\delta V_{\mathrm{ps}}^{l}\rangle\frac{1}{\langle\varPsi_{\mathrm{ps}}^{lm}\,|\,\delta V_{\mathrm{ps}}^{l}\,|\,\varPsi_{\mathrm{ps}}^{lm}\rangle}\langle\delta V_{\mathrm{ps}}^{l}\varPsi_{\mathrm{ps}}^{lm}\,|\,\varphi_j\rangle \tag{3.272} $$
显然，这正是方程（3.269）的形式。两个基函数的积分是分开进行的，因此对于给定的 $l$，计算次数下降为 $NM$，对实际工作而言，这是一个极大的改进。
引入 $\chi_{\mathrm{ps}}^{lm}$，有
$$ |\,\chi_{\mathrm{ps}}^{lm}\rangle=|\,\delta V_{\mathrm{ps}}^{l}\varPsi_{\mathrm{ps}}^{lm}\rangle \tag{3.273} $$
则可将 KB 非局域赝势形式改写为更简洁的形式：
$$ \delta V_{\mathrm{NL}}=\sum_{l=0}^{l_{\max}}\sum_{m=-l}^{l}\frac{|\,\chi_{\mathrm{ps}}^{lm}\rangle\langle\chi_{\mathrm{ps}}^{lm}\,|}{\langle\chi_{\mathrm{ps}}^{lm}\,|\,\varPsi_{\mathrm{ps}}^{lm}\rangle} \tag{3.274} $$
式（3.274）表明，可以在构建模守恒赝势 $V_{\mathrm{ps}}^{l}$ 的同时得到 KB 非局域赝势形式（因为可以同时得到 $|\,\varPsi_{\mathrm{ps}}^{lm}\rangle$ 与 $\delta V_{\mathrm{ps}}^{l}$）。不难看到，新引入的 $\chi_{\mathrm{ps}}^{lm}$ 满足
$$ \chi_{\mathrm{ps}}^{lm}(\boldsymbol{r})=\left[\varepsilon_l-\left(-\frac{1}{2}\boldsymbol{\nabla}^{2}+V_{\mathrm{ps}}^{\mathrm{loc}}(r)\right)\right]\varPsi_{\mathrm{ps}}^{lm}(\boldsymbol{r}) \tag{3.275} $$
显然，$\chi_{\mathrm{ps}}^{lm}$ 只在 $r_{\mathrm{c}}$ 内不为零，因此是一个局域函数。而赝势波函数满足：
$$ \left(-\frac{1}{2}\boldsymbol{\nabla}^{2}+V_{\mathrm{ps}}^{\mathrm{loc}}(r)+\delta V_{\mathrm{ps}}^{l}\right)\varPsi_{\mathrm{ps}}^{lm}(\boldsymbol{r})=\varepsilon_l\varPsi_{\mathrm{ps}}^{lm}(\boldsymbol{r}) \tag{3.276} $$
KB 非局域赝势形式的重要性还在于它为赝势构建新方法的应用铺平了道路。如式（3.274）所示，$\delta V_{\mathrm{NL}}$ 的形式与 OPW 方法中的 $\hat{V}_{\mathrm{PK}}$（见式（3.239））非常相似。这表明赝势可以表示为投影算符，而不是必须采用 3.3.3 节一开始介绍的半局域形式。在 KB 方法中，对于给定的 $l$，$\delta V_{\mathrm{NL}}^{l}$ 只扩展为一项投影算符。如果将其扩展为多项投影算符的线性组合，会产生怎样的效果呢？Blöchl 和 Vanderbilt 各自独立地研究了这个问题\cite{blochl1990generalized,vanderbilt1990soft}，并由此为赝势家族增添了超软赝势和投影缀加平面波两大类成员。这两类赝势也被广泛地应用于当前流行的电子结构计算软件中。
超软赝势
3.3.2 节中引入了模守恒条件，用以保证所生成的赝势可以适用于不同环境。但是在特定情况下，这个限制条件可能会影响计算效率。考虑元素周期表中第二行元素的 $2p$ 轨道（如 O2p）或者第四周期过渡金属元素的 $3d$ 轨道（如 Cu3d）等，因为原子径向波函数的节点数为 $n-l-1$，所以这两类价电子轨道没有节点。如果保持模守恒条件，可以预见赝波函数 $\varPsi_{\mathrm{ps}}^{l}$ 与全电子波函数 $\varPsi_{\mathrm{ae}}^{l}$ 形状相似。因此，在平面波计算方法中，仍然需要大量的平面波来展开这类赝波函数，这样将使计算效率降低。
为了解决这个问题，Vanderbilt 提出，可以放弃模守恒条件，从而生成对应于更平滑的赝波函数 $\varPsi_{\mathrm{ps}}^{l}$ 的赝势，称为超软赝势（ultrasoft pseudopotential，USPP）\cite{vanderbilt1990soft}。为了再次满足 3.3.2 节中提出的条件（5），USPP 需要引入额外的补偿项。尽管这些补偿项的引入可能导致体系总能以及原子受力的表达式的变化，但是实际上，补偿项只需要在构建赝势时计算一次即可，而在后续计算中保持不变。此外，能量和力的附加项可以与正常项同时求解，因此利用 USPP 可以有效降低哈密顿矩阵的维数，从而提高计算效率。下面将具体介绍这一方法。
对于给定的 $l$ 和 $m$，选定 $s$ 个（通常为 1～3 个）参考能量 $\varepsilon_i^{lm}$，对于每一个 $\varepsilon_i^{lm}$，求解薛定谔方程\cite{hamann1989generalized}
$$ \left(-\frac12\boldsymbol\nabla^2+V_{\mathrm{ae}}\right)\varPsi_{\mathrm{ae},i}^{lm}=\varepsilon_i^{lm}\varPsi_{\mathrm{ae},i}^{lm} \tag{3.277} $$
然后按照 3.3.2 节中介绍的普遍过程构造 $\varPsi_{\mathrm{ps},i}^{lm}$，并通过方程（3.275）计算 $\chi_i^{lm}$：
$$ |\,\chi_i^{lm}\rangle=\left[\varepsilon_i^{lm}-\left(-\frac{1}{2}\boldsymbol{\nabla}^{2}+V_{\mathrm{ps}}^{\mathrm{loc}}(r)\right)\right]|\,\varPsi_{\mathrm{ps},i}^{lm}\rangle \tag{3.278} $$
至此，对于给定的 $(l,m)$，我们得到了两类赝波函数：$\{\varPsi_{\mathrm{ps},i}^{lm}\}$ 和 $\{\chi_i^{lm}\}$。由此可以构建矩阵 $\boldsymbol{B}_{ij}^{lm}$：
$$ \boldsymbol{B}_{ij}^{lm}=\langle\varPsi_{\mathrm{ps},i}^{lm}\,|\,\chi_j^{lm}\rangle \tag{3.279} $$
$\boldsymbol{B}_{ij}^{lm}$ 是一个 $s$ 阶方阵。因为 $|\,\varPsi_{\mathrm{ps},i}^{lm}\rangle$ 和 $|\,\chi_{\mathrm{ps},i}^{lm}\rangle$ 总是选取同样的 $m$，所以 $\boldsymbol{B}_{ij}^{lm}$ 上标中的 $m$ 可以省略\cite{kohanoff2006electronic}。但是为了简化公式，在此仍然将其保留。通过 $\boldsymbol{B}_{ij}^{lm}$，可以定义新的局域波函数：
$$ |\beta_i^{lm}\rangle=\sum_j(B^{-1})_{ji}^{lm}|\chi_j^{lm}\rangle \tag{3.280} $$
此外，还可以定义一个补偿量 $Q_{ij}^{lm}$：
$$ Q_{ij}^{lm}=\langle\varPsi_{\mathrm{ae},i}^{lm}\,|\,\varPsi_{\mathrm{ae},j}^{lm}\rangle_{R_{\mathrm{C}}}-\langle\varPsi_{\mathrm{ps},i}^{lm}\,|\,\varPsi_{\mathrm{ps},j}^{lm}\rangle_{R_{\mathrm{C}}} \tag{3.281} $$
式中：$\langle\cdots\rangle_{R_{\mathrm{C}}}$ 代表积分在半径为 $R_{\mathrm{C}}$ 的球体内进行。Vanderbilt 在文献\cite{vanderbilt1990soft}中证明，若 $Q_{ij}^{lm}=0$，则矩阵 $\boldsymbol{B}_{ij}^{lm}$ 是厄米矩阵，可将赝势的非局域形式写为
$$ \delta V_{\mathrm{NL}}^{l}=\sum_{m=-l}^{l}\sum_{ij}\boldsymbol{B}_{ij}^{lm}\,|\,\beta_i^{lm}\rangle\langle\beta_j^{lm}\,| \tag{3.282} $$
对比式（3.282）与 KB 非局域赝势形式的方程（3.274），可以看到前者正是后者的推广，因为式（3.282）同样将赝势的非局域部分 $\delta V_{\mathrm{NL}}$ 写为投影算符。与 KB 非局域赝势不同的是，$\delta V_{\mathrm{NL}}$ 在这里被表示为 $s$ 个投影算符的线性组合。
但是 $Q_{ij}^{lm}$ 并非必须为零不可。若 $Q_{ij}^{lm}\neq0$，则对薛定谔方程的求解由本征值问题转化为广义本征值问题。由此，需要定义交叠算符：
$$ \hat{\mathrm{S}}=\boldsymbol{I}+\sum_{lm}\sum_{ij}Q_{ij}^{lm}\,|\,\beta_i^{lm}\rangle\langle\beta_j^{lm}\,| \tag{3.283} $$
不难证明，交叠算符 $\hat{\mathrm{S}}$ 具有以下性质：
$$ \langle\varPsi_{\mathrm{ps},i}^{lm}\,|\,\hat{\mathrm{S}}\,|\,\varPsi_{\mathrm{ps},i}^{lm}\rangle_{R_{\mathrm{C}}}=\langle\varPsi_{\mathrm{ae},i}^{lm}\,|\,\varPsi_{\mathrm{ae},i}^{lm}\rangle_{R_{\mathrm{C}}} \tag{3.284} $$
此外，定义算符 $D_{ij}^{lm}$ 为
$$ D_{ij}^{lm}=B_{ij}^{lm}+\varepsilon_jQ_{ij}^{lm} \tag{3.285} $$
并借此定义超软赝势的非局域形式为
$$ \delta V_{\mathrm{NL}}^{\mathrm{US}}=\sum_{lm}\sum_{ij}D_{ij}^{lm}\,|\,\beta_i^{lm}\rangle\langle\beta_j^{lm}\,| \tag{3.286} $$
由式（3.278）至式（3.286）可得，赝波函数 $\varPsi_{\mathrm{ps},i}^{lm}$ 满足
$$ \left(-\frac{1}{2}\boldsymbol{\nabla}^{2}+V_{\mathrm{ps}}^{\mathrm{loc}}+\delta V_{\mathrm{NL}}^{\mathrm{US}}-\varepsilon_i\hat{\mathrm{S}}\right)|\,\varPsi_{\mathrm{ps},i}^{lm}\rangle=0 \tag{3.287} $$
将模守恒条件 $Q_{ij}^{lm}=0$ 取消，意味着赝波函数的限制条件只有“在芯区半径 $r_{\mathrm{c}}$ 终点处及以外与全电子波函数一致”。这种宽松的条件使得超软赝势可以选取非常大的 $r_{\mathrm{c}}$，从而有效改善芯区内的波函数的平滑度，这无疑会提高对元素周期表中第二行元素及第三行过渡金属元素的计算效率。
另一方面，在具体的计算中，因为模守恒条件被取消，所以价电子密度“缺失”的部分需要用 $Q_{ij}^{lm}$ 来补偿，即
$$ \rho_{\mathrm{v}}(\boldsymbol{r})=\sum_{n}^{\mathrm{occ}}\varphi_n^{*}(\boldsymbol{r})\varphi_n(\boldsymbol{r})+\sum_{lm}\sum_{ij}\rho_{ij}^{lm}Q_{ij}^{lm}(\boldsymbol{r}) \tag{3.288} $$
式中
$$ \rho_{ij}^{lm}=\sum_{n}^{\mathrm{occ}}\langle\varphi_n\,|\,\beta_j^{lm}\rangle\langle\beta_i^{lm}\,|\,\varphi_n\rangle \tag{3.289} $$
$$ Q_{ij}^{lm}(\boldsymbol{r})=\varPsi_{\mathrm{ae},i}^{lm*}(\boldsymbol{r})\varPsi_{\mathrm{ae},j}^{lm}(\boldsymbol{r})-\varPsi_{\mathrm{ps},i}^{lm*}(\boldsymbol{r})\varPsi_{\mathrm{ps},j}^{lm}(\boldsymbol{r}) \tag{3.290} $$
方程（3.288）中的 $\varphi_n$ 满足广义正交性条件
$$ \langle\varphi_m\,|\,\hat{\mathrm{S}}\,|\,\varphi_n\rangle=\delta_{mn} \tag{3.291} $$
而体系的总能为
$$ \begin{aligned} E_{\mathrm{tot}}={}&\sum_{n=1}^{\mathrm{occ}}\langle\varphi_n\,|\left(-\frac{1}{2}\boldsymbol{\nabla}^{2}+V_{\mathrm{ps}}^{\mathrm{loc}}+\sum_{lm}\sum_{ij}D_{ij}^{lm}\,|\,\beta_i^{lm}\rangle\langle\beta_j^{lm}\,|\right)|\,\varphi_n\rangle\\ &+E_{\mathrm{H}}[\rho_{\mathrm{v}}]+E_{\mathrm{xc}}[\rho_{\mathrm{v}}]+E_{\mathrm{II}} \end{aligned} \tag{3.292} $$
在广义正交性条件（式（3.291））下求 $E_{\mathrm{tot}}$ 的变分极值，所得的结果即为式（3.288）至式（3.292）中出现的本征波函数 $\varphi_n$。
为了简化最后的久期方程，引入包含全体系 Hartree 势和交换关联势的局域有效势 $V_{\mathrm{eff}}^{\mathrm{loc}}$，以及修正后的投影系数 $\tilde{D}_{ij}^{lm}$，即
$$ V_{\mathrm{eff}}^{\mathrm{loc}}(\boldsymbol r) =\sum_I V_{\mathrm{ps},I}^{\mathrm{loc}}(\boldsymbol r-\boldsymbol R_I) +V_{\mathrm H}(\boldsymbol r)+V_{\mathrm{xc}}(\boldsymbol r) \tag{3.293} $$
$$ \tilde{D}_{ij}^{lm}=D_{ij}^{lm}+\int\mathrm{d}\boldsymbol{r}\,(V_{\mathrm{H}}(\boldsymbol{r})+V_{\mathrm{xc}}(\boldsymbol{r}))Q_{ij}^{lm}(\boldsymbol{r}) \tag{3.294} $$
由此，可以将利用超软赝势的久期方程写为
$$ \left[-\frac12\nabla^2+V_{\mathrm{eff}}^{\mathrm{loc}}(\boldsymbol r) +\sum_I\sum_{ij}\tilde D_{ij}^{I}|\beta_i^I\rangle\langle\beta_j^I|\right] |\varphi_n\rangle=\varepsilon_n\hat S|\varphi_n\rangle \tag{3.295} $$
式中，$I$ 标记原子中心及其角动量通道；投影系数由式（3.286）的 $D_{ij}^{lm}$ 替换为式（3.294）的 $\tilde D_{ij}^{lm}$ 得到。
近年来，Blöchl 提出的投影缀加平面波（projector augmented-wave，PAW）方法\cite{blochl1994projector}引起了越来越多的关注。与传统的赝势方法相比，PAW 方法最大的优点是可以重新构建出因为赝势化而丢失的芯区电子的信息，而其构建过程并不比 USPP 的构建过程复杂。事实上，USPP 方程与 PAW 方程有着类似的推导过程。Kresse 与 Joubert 给出了二者之间联系的详细证明\cite{kresse1999from}。受篇幅所限，我们在这里不展开讨论。