本文是「第一性原理的微观计算模拟」系列的第 8 篇（共 11 篇），内容整理自同名书稿第 3 章，公式编号与原书一致。文中引用的文献依据原书参考文献清单整理并列于文末，按原书清单顺序从 1 开始编号。\nocite{*}
Slater 于 1937 年提出，可以将体系分成两部分，即以原子核为球心的球形邻域（I区）和各个球形领域之间的空隙区（II区）\cite{slater1937wave}。I区受离子势影响强烈，因此电子波函数变化比较剧烈，近似于原子轨道，而II区受离子势影响较小，电子波函数相应变化比较平缓。基于这种考虑，Slater 提出，可以将体系的本征波函数在I区中用局域化轨道波函数展开，在II区中则用平面波展开，然后通过在边界处函数值连续的条件将这两部分结合起来。这种用混合基函数展开本征波函数的方法称为缀加平面波（augmented plane wave，APW）方法。
APW 方法的理论基础及公式推导
本节讨论 APW 方法的矩阵元构成，其数学推导比较繁难。在下面的讨论中，我们采取了 Hartree 原子单位制，以避免 $\hbar$、$m_{\mathrm{e}}$ 等常数频繁出现。为了进一步简化模型，I区中的势函数取球对称形式，II区中的势场取常函数，即
$$ V(u)=\begin{cases}V(u)=-Z/u,&u\leqslant r_s\\V_0,&u\gt r_s\end{cases} \tag{3.522} $$
式中：$r_s$ 为单胞内位于 $\tau_s$ 处的第 $s$ 个原子的I区半径。方程（3.522）所描述的势函数被称为松糕势（muffin-tin potential，MT）。图 3.13 给出了由两种原子组成的二维带心正方格子的 MT 势示意图，其中 $Z_1=7.0q$，$Z_2=5.6q$，$q$ 为单位正电荷。I区半径分别为 $r_1$、$r_2$。两种原子相对于晶胞原点的坐标分别为 $\tau_1=0$ 和 $\tau_2$。相应地，在I区内的基函数 $\chi_{\mathrm{I}}$ 表示为类氢波函数的线性组合：
$$ \chi_{\mathrm{I}}(\boldsymbol{r})=\sum_{s}\sum_{l=0}\sum_{m=-l}^{l}A_{lm}^{s}R_l(u)\mathrm{Y}_l^{m}(\theta,\phi) \tag{3.523} $$
角向部分为球谐函数，而径向部分满足
$$ \left\{-\frac{1}{2}\frac{\mathrm{d}^{2}}{\mathrm{d}u^{2}}+\frac{l(l+1)}{2u^{2}}+V(u)\right\}uR_l(u)=E'uR_l(u) \tag{3.524} $$
式中：$E'$ 为任意参数。我们选择用 $\boldsymbol{u}$ 来描述径向部分，因为对于给定的原点，空间中 $\boldsymbol{r}$ 终点处的点相对于不同的原子核有不同的距离和方位。显然 $\boldsymbol{r}$ 和 $\boldsymbol{u}$ 之间有如下关系：
$$ \boldsymbol{r}=\boldsymbol{u}+\boldsymbol{\tau}_s \tag{3.525} $$
II区中因为势函数恒等于 $V_0$，所以基函数 $\chi_{\mathrm{II}}$ 可用平面波展开：
$$ \chi_{\mathrm{II}}=\frac{1}{\sqrt{\varOmega}}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}}=\frac{1}{\sqrt{\varOmega}}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{\tau}_s}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{u}}=\frac{4\pi}{\sqrt{\varOmega}}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{\tau}_s}\sum_{l=0}\sum_{m=-l}^{l}i^{l}\mathrm{j}_l(ku)\mathrm{Y}_l^{m*}(\hat{\boldsymbol{k}})\mathrm{Y}_l^{m}(\hat{\boldsymbol{u}}) \tag{3.526} $$
图 3.13 由两种原子组成的二维带心正方格子的 MT 势
在周期性体系之内，体系的本征波函数 $\psi_{\boldsymbol{k}}(\boldsymbol{r})$ 为 Bloch 波函数，也即 $\psi_{\boldsymbol{k}}(\boldsymbol{r})$ 可展开为
$$ \psi_{\boldsymbol{k}}(\boldsymbol{r})=\sum_{i=1}^{M}c_{i,\boldsymbol{k}}\chi(\boldsymbol{r},\boldsymbol{k}+\boldsymbol{G}_i) \tag{3.527} $$
式中：$\boldsymbol{G}$ 是倒格矢。在球面处要求 $\chi_{\mathrm{I}}$ 与 $\chi_{\mathrm{II}}$ 函数值连续。因此由式（3.523）和式（3.526）可知，系数 $A_{lm}$ 满足
$$ A_{lm}^{s,\boldsymbol{K}_i}=\frac{4\pi}{\sqrt{\varOmega}}\mathrm{e}^{\mathrm{i}\boldsymbol{K}_i\cdot\boldsymbol{\tau}_s}i^{l}\mathrm{Y}_l^{m*}(\hat{\boldsymbol{K}}_i)\frac{\mathrm{j}_l(K_ir_s)}{R_l(r_s)} \tag{3.528} $$
为了简化公式，设 ${ \boldsymbol{K}_i=\boldsymbol{k}+\boldsymbol{G}_i}$。至此，可以写出 APW 方法中的基函数 $\chi(\boldsymbol{K}_i,\boldsymbol{r})$：
$$ \chi(\boldsymbol{K}_i,\boldsymbol{r})=\begin{cases}\displaystyle\sum_{s}\frac{4\pi}{\sqrt{\varOmega}}\mathrm{e}^{\mathrm{i}\boldsymbol{K}_i\cdot\boldsymbol{\tau}_s}\sum_{l=0}^{\infty}\sum_{m=-l}^{l}i^{l}\mathrm{Y}_l^{m*}(\hat{\boldsymbol{K}}_i)\mathrm{Y}_l^{m}(\hat{\boldsymbol{u}})\mathrm{j}_l(K_ir_s)\frac{R_l(u)}{R_l(r_s)},&u\leqslant r_s\\\dfrac{1}{\sqrt{\varOmega}}\mathrm{e}^{\mathrm{i}\boldsymbol{K}_i\cdot\boldsymbol{r}},&u\gt r_s\end{cases} \tag{3.529} $$
图 3.14 给出了 $\chi(\boldsymbol{r})$ 的一个例子。图中竖直虚线表示I区与II区的边界，大、小实心原点分别表示两种原子，其电荷电量分别为 $7.0q$ 和 $5.6q$。可以看到，在I区与II区边界处，$V(\boldsymbol{r})$ 和 $\chi(\boldsymbol{r})$ 虽然连续，但是并不平滑。
![图 3.14 沿图 3.13 中[110]方向（沿图中 $\tau_2$ 方向）的 MT 势函数 $V(\boldsymbol{r})$ 及基函数 $\chi(\boldsymbol{r})$](assets/fig3-14.jpg)
图 3.14 沿图 3.13 中[110]方向（沿图中 $\tau_2$ 方向）的 MT 势函数 $V(\boldsymbol{r})$ 及基函数 $\chi(\boldsymbol{r})$
（a）MT 势函数 $V(\boldsymbol{r})$；（b）基函数 $\chi(\boldsymbol{r})$
方程（3.529）保证了基函数在全空间内的连续性，但是并不能保证在边界处的平滑性，即在各个I区与II区的边界处，$\chi(\boldsymbol{K}_i,\boldsymbol{u})$ 的左导数不等于右导数。这种导函数的不连续性使得动量算符在每一个边界处均对体系总能 $E$ 产生一项额外的贡献 $T_{\mathrm{S}}$，第 $s$ 个边界处的 $T_{\mathrm{S}}^{s}$ 在哈密顿矩阵中的矩阵元为
$$ T_{{\mathrm S},ij}^{s} =\frac12\oint_{|\boldsymbol u|=r_s}\!\chi_i^*(\boldsymbol r) \left[\partial_u\chi_{j,\mathrm I}(\boldsymbol r) -\partial_u\chi_{j,\mathrm{II}}(\boldsymbol r)\right]\mathrm dS \tag{3.530} $$
式中：$\boldsymbol{G}_{ij}=\boldsymbol{K}_j-\boldsymbol{K}_i$；积分域 $\varOmega_\varepsilon$ 是指与球形边界 $s$ 同心，且内、外径分别为 $r_s-\varepsilon$ 和 $r_s+\varepsilon$ 的球壳，$\varepsilon\to0$。式（3.530）可以计算如下：
$$ \begin{aligned} -\frac12\chi_i^*\nabla^2\chi_j &=-\frac12\nabla\cdot(\chi_i^*\nabla\chi_j) +\frac12\nabla\chi_i^*\cdot\nabla\chi_j,\\ T_{{\mathrm S},ij}^{s} &=\frac12\oint_{|\boldsymbol u|=r_s}\!\chi_i^* \bigl(\partial_u\chi_{j,\mathrm I}-\partial_u\chi_{j,\mathrm{II}}\bigr)\,\mathrm dS. \end{aligned} \tag{3.531} $$
在式（3.531）里，第一项利用格林公式变成了沿边界的面积分，并且考虑到 $\varOmega_\varepsilon$ 内、外表面的法线方向彼此相反，而第二项在 $\varepsilon\to0$ 的情况下可忽略不计。
计入 $T_{\mathrm{S}}$ 之后，可以写出体系的总能 $E$ 的表达式：
$$ E=\frac{\sum_{ij}c_i^*\left(H_{ij}+\sum_sT_{{\mathrm S},ij}^{s}\right)c_j} {\sum_{ij}c_i^*S_{ij}c_j},\qquad H_{ij}=\int_{\varOmega}\chi_i^*\hat H\chi_j\,\mathrm d^3\boldsymbol r, \quad S_{ij}=\int_{\varOmega}\chi_i^*\chi_j\,\mathrm d^3\boldsymbol r \tag{3.532} $$
式中：$\hat{H}=-\boldsymbol{\nabla}^{2}/2+V(\boldsymbol{r})$。下文的 $\Delta_{ij}$ 与此处 $S_{ij}$ 均表示基函数重叠矩阵，即 $\Delta_{ij}=S_{ij}$。将式（3.527）及式（3.531）代入式（3.532），并适当移项，得
$$ \sum_{ij}c_i^*\left(H_{ij}-ES_{ij}+\sum_s T_{{\mathrm S},ij}^{s}\right)c_j=0 \tag{3.533} $$
可以证明，式（3.533）确实是总能的变分表达式，即 $E$ 的极值可在式（3.533）对 $\psi$ 求变分极值时得到。证明的过程比较烦琐，具体请参看文献\cite{loucks1967augmented,schlosser1963composite}。因此，为了求得本征函数，需要满足
$$ \frac{\partial E}{\partial c_i^{*}}=0,\quad i=1,2,\cdots,M \tag{3.534} $$
将式（3.533）代入式（3.534），得线性方程组
$$ \sum_{j}^{M}\left(H_{ij}-E\Delta_{ij}+\sum_{s}T_{\mathrm{S},ij}^{s}\right)c_j=0,\quad i=1,2,\cdots,M \tag{3.535} $$
式中：$T_{\mathrm{S},ij}^{s}$ 由式（3.531）给出，且
$$ H_{ij}=\int_{\varOmega}\chi^{*}(\boldsymbol{K}_i,\boldsymbol{r})\hat{H}\chi(\boldsymbol{K}_j,\boldsymbol{r})\,\mathrm{d}\boldsymbol{r} \tag{3.536} $$
$$ \Delta_{ij}=\int_{\varOmega}\chi^{*}(\boldsymbol{K}_i,\boldsymbol{r})\chi(\boldsymbol{K}_j,\boldsymbol{r})\,\mathrm{d}\boldsymbol{r} \tag{3.537} $$
因此，APW 方法归结为求解久期方程
$$ \det\left|H_{ij}-E\Delta_{ij}+\sum^{s}T_{\mathrm{S},ij}^{s}\right|=0 \tag{3.538} $$
因为体系分为I区和II区，$\chi$ 在两个区域中各不相同。所以 $H_{ij}-E\Delta_{ij}$ 也可分别写为两个区域的贡献：$H_{ij}^{\mathrm{I}}-E\Delta_{ij}^{\mathrm{I}}$ 和 $H_{ij}^{\mathrm{II}}-E\Delta_{ij}^{\mathrm{II}}$。从原则上讲，至此可以计算哈密顿矩阵的矩阵元。但是考虑以下两个因素可以极大地简化计算过程。
首先，根据式（3.524），可得
$$ H_{ij}^{\mathrm{I}}-E\Delta_{ij}^{\mathrm{I}}=(E'-E)\Delta_{ij}^{\mathrm{I}} \tag{3.539} $$
因此，若取 $E'=E$，则 $H_{ij}-E\Delta_{ij}$ 在I区内的贡献为零，这无疑可以极大地简化矩阵元的计算，但是会导致 $\chi$ 成为体系能量 $E$ 的隐函数，因此无法通过常规的矩阵对角化求解本征值和本征函数。另一方面，这种 $E'=E$ 的强制选择将使得 APW 方法在求解过程中内秉地调节试探波函数，直至对于给定的势函数找到最优解为止，一般认为，这一特点正是 APW 方法取得普遍成功的原因。
其次，II区内的 $H_{ij}-E\Delta_{ij}$ 计算比较困难，但是因为 $\chi_{\mathrm{II}}$ 是平面波，所以在全空间内的积分满足正交条件。因此，将 $H_{ij}^{\mathrm{II}}-E\Delta_{ij}^{\mathrm{II}}$ 写为如下形式：
$$ H_{ij}^{\mathrm{II}}-E\Delta_{ij}^{\mathrm{II}}=(H_{ij}^{\mathrm{II}}-E\Delta_{ij}^{\mathrm{II}})_{\varOmega}-\sum_{s}(H_{ij}^{\mathrm{II}}-E\Delta_{ij}^{\mathrm{II}})_{\text{sphere-}s} \tag{3.540} $$
也即首先计算平面波在全空间的积分（意味着I区为空区，即势函数 $V\equiv0$），然后减去I区中的贡献。式（3.540）右端第一项很容易求出：
$$ (H_{ij}^{\mathrm{II}}-E\Delta_{ij}^{\mathrm{II}})_{\varOmega}=\frac{1}{\varOmega}\int_{\varOmega}\mathrm{e}^{-\mathrm{i}\boldsymbol{K}_i\cdot\boldsymbol{r}}\left(-\frac{\boldsymbol{\nabla}^{2}}{2}-E\right)\mathrm{e}^{\mathrm{i}\boldsymbol{K}_j\cdot\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r}=\left(\frac{|\,\boldsymbol{K}_j\,|^{2}}{2}-E\right)\delta_{ij} \tag{3.541} $$
而第 $s$ 个球内的贡献可计算如下：
$$ \begin{aligned} (H_{ij}^{\mathrm{II}}-ES_{ij}^{\mathrm{II}})_{\mathrm{sphere}\,s} ={}&\left(\frac{|\boldsymbol K_j|^2}{2}-E\right) \frac{e^{\mathrm i\boldsymbol G_{ij}\cdot\boldsymbol\tau_s}}{\varOmega} \int_{|\boldsymbol u|\lt r_s}e^{\mathrm i\boldsymbol G_{ij}\cdot\boldsymbol u}\,\mathrm d^3\boldsymbol u\\ ={}&\left(\frac{|\boldsymbol K_j|^2}{2}-E\right) \frac{4\pi r_s^2e^{\mathrm i\boldsymbol G_{ij}\cdot\boldsymbol\tau_s}}{\varOmega} \frac{j_1(|\boldsymbol G_{ij}|r_s)}{|\boldsymbol G_{ij}|}, \end{aligned} \tag{3.542} $$
当 $|\boldsymbol G_{ij}|=0$ 时，上式中的 $j_1(|\boldsymbol G_{ij}|r_s)/|\boldsymbol G_{ij}|$ 按连续极限取值 $r_s/3$。
因此，由式（3.539）至式（3.542）可得
$$ H_{ij}^{\mathrm{II}}-ES_{ij}^{\mathrm{II}} =\left(\frac{|\boldsymbol K_j|^2}{2}-E\right) \left[\delta_{ij}-\sum_s\frac{4\pi r_s^2}{\varOmega} e^{\mathrm i\boldsymbol G_{ij}\cdot\boldsymbol\tau_s} \frac{j_1(|\boldsymbol G_{ij}|r_s)}{|\boldsymbol G_{ij}|}\right] \tag{3.543} $$
矩阵元中面积分的贡献 $T_{\mathrm{S},ij}^{s}$ 分为内表面积分和外表面积分两项，这两项也是体系总能 $E$ 的隐函数。两项积分的求解过程类似，所以这里将它们放在一起讨论。因为在表面上 $\chi$ 连续，所以方程（3.531）中的 $\chi^{*}(\boldsymbol{K}_i,\boldsymbol{u})$ 为
$$ \chi^{*}(\boldsymbol{K}_i,\boldsymbol{u})\,|_{\text{sphere-}s}=\frac{4\pi\mathrm{e}^{-\mathrm{i}\boldsymbol{K}_i\cdot\boldsymbol{\tau}_s}}{\sqrt{\varOmega}}\sum_{l}\sum_{m=-l}^{l}(-i)^{l}\mathrm{j}_l(r_s\,|\,\boldsymbol{K}_i\,|)\mathrm{Y}_l^{m}(\hat{\boldsymbol{K}}_i)\mathrm{Y}_l^{m*}(\hat{\boldsymbol{u}}) \tag{3.544} $$
而
$$ \begin{aligned} \frac{\partial}{\partial u}[\chi_{\mathrm{II}}(\boldsymbol{K}_j,\boldsymbol{u})-\chi_{\mathrm{I}}(\boldsymbol{K}_j,\boldsymbol{u})]={}&\frac{4\pi\mathrm{e}^{\mathrm{i}\boldsymbol{K}_j\cdot\boldsymbol{\tau}_s}}{\sqrt{\varOmega}}\sum_{l}\sum_{m=-l}^{l}i^{l}\mathrm{j}_l(r_s\,|\,\boldsymbol{K}_j\,|)\mathrm{Y}_l^{m*}(\hat{\boldsymbol{K}}_j)\mathrm{Y}_l^{m}(\hat{\boldsymbol{u}})\\ &\times\left[\frac{|\,\boldsymbol{K}_j\,|\,\mathrm{j}_l'(|\,\boldsymbol{K}_j\,|\,u)}{\mathrm{j}_l(|\,\boldsymbol{K}_j\,|\,r_s)}-\frac{R_l'(u)}{R_l(r_s)}\right] \end{aligned} \tag{3.545} $$
其中 $\mathrm{j}_l'(x)=\mathrm{d}\mathrm{j}_l(x)/\mathrm{d}x$，而 $x=|\,\boldsymbol{K}_j\,|\,u$。面积分的积分元
$$ \mathrm{d}S=r_s^{2}\sin\theta\,\mathrm{d}\theta\,\mathrm{d}\phi=r_s^{2}\,\mathrm{d}\hat{\boldsymbol{u}} \tag{3.546} $$
此外还有球谐函数的正交关系
$$ \int\mathrm{Y}_l^{m*}(\hat{\boldsymbol{u}})\mathrm{Y}_{l'}^{m'}(\hat{\boldsymbol{u}})\,\mathrm{d}\hat{\boldsymbol{u}}=\delta_{ll'}\delta_{mm'} \tag{3.547} $$
$$ \sum_{m=-l}^{l}\mathrm{Y}_l^{m*}(\hat{\boldsymbol{K}}_i)\mathrm{Y}_l^{m}(\hat{\boldsymbol{K}}_j)=\frac{2l+1}{4\pi}\mathrm{P}_l(\cos\theta_{\boldsymbol{K}_i\boldsymbol{K}_j}) \tag{3.548} $$
式中：$\theta_{\boldsymbol{K}_i\boldsymbol{K}_j}$ 表示 $\boldsymbol{K}_i$ 和 $\boldsymbol{K}_j$ 间的夹角。
将式（3.544）至式（3.547）代入 $T_{\mathrm{S},ij}^{s}$ 的表达式（3.531），可得
$$ T_{{\mathrm S},ij}^{s} =\frac12\oint_{|\boldsymbol u|=r_s}\!\chi_i^* \left(\partial_u\chi_{j,\mathrm I}-\partial_u\chi_{j,\mathrm{II}}\right)\mathrm dS \tag{3.549} $$
至此，我们得出了 APW 方法久期方程的矩阵元 $M_{ij}$：
$$ M_{ij}(E)=H_{ij}^{\mathrm{II}}-ES_{ij}^{\mathrm{II}} +\sum_s\left[H_{ij}^{\mathrm I,s}-ES_{ij}^{\mathrm I,s}+T_{{\mathrm S},ij}^{s}\right] \quad\text{（各球内通道分别求和）} \tag{3.550} $$
球贝塞尔函数满足加法定理\cite{grosso2000solid}
$$ \sum_{l=0}^{\infty}(2l+1)P_l(\hat{\boldsymbol K}_i\!\cdot\!\hat{\boldsymbol K}_j) j_l(K_ir_s)j_l(K_jr_s) =j_0(|\boldsymbol K_i-\boldsymbol K_j|r_s) \tag{3.551} $$
对分段连续的 APW 基函数，采用包含球面跳跃项的变分二次型可直接得到厄米矩阵\cite{slater1937wave,loucks1967augmented}：
$$ M_{ij}(E)=\frac12\int_{\varOmega}\nabla\chi_i^*\cdot\nabla\chi_j\,\mathrm d^3\boldsymbol r +\int_{\varOmega}\chi_i^*V\chi_j\,\mathrm d^3\boldsymbol r -E\int_{\varOmega}\chi_i^*\chi_j\,\mathrm d^3\boldsymbol r, \qquad M_{ji}(E)=M_{ij}(E)^* \tag{3.552} $$
图 3.15 用 APW 方法求解本征值示意图
前面的讨论已经指出，因为 $M_{ij}$ 是待求的 $E$ 的隐函数，所以求解方程（3.538）的通常做法是给定倒空间中的一点 $\boldsymbol{k}$，然后改变 $E$，对于每一个 $E$ 的取值，通过式（3.524）、式（3.550）求得每一个 $M_{ij}$，直到找到 $\det|\,M_{ij}\,|=0$ 的 $E_n$，此即体系第 $n$ 条能带的本征值，继而可得相应的本征波函数 $\psi_{n,\boldsymbol{k}}(\boldsymbol{r})$。按上述方法找到的 $N$ 个解即构成 $\boldsymbol{k}$ 处的一套本征能级。当 $\boldsymbol{k}$ 遍历第一布里渊区时，也可由 APW 方法得出体系完整的能带结构，图 3.15 为用 APW 方法求解本征值示意图。可见，APW 方法每次只能更新一条能带，且对于倒空间内的所有 $\boldsymbol{k}$ 点都需重复进行上述步骤，因此其效率是很低的。在对其进行线性化处理之后，APW 方法的计算效率有很大提升。
APW 方法的线性化处理
APW 方法在计算上的主要困难是由面积分 $T_{\mathrm{S}}$ 中的 $R_l'(r_s)/R(r_s)$ 项造成的，这一项一般写为对数导数 $\partial\ln R_l(u)/\partial u\,|_{u=r_s}$。对数项的存在使得 APW 方法的基函数成为关于能量的非线性函数，这意味着基函数依赖于待求能量 $E$，从而导致了求解上的困难以及大计算量。为了解决这个困难，Andersen 在 1975 年提出对 APW 方法中I区基函数的线性化处理方法——线性缀加平面波（linearized APW，LAPW）方法\cite{andersen1975linear}。
LAPW 方法中，径向函数改由下式确定：
$$ R_l(u,E)=R_l(u,E_\nu)+(E-E_\nu)\dot{R}_l(u,E_\nu)+\cdots \tag{3.553} $$
式中：$\dot{R}_l(u,E_\nu)=\partial R_l(u,E_\nu)/\partial E$。将 $R_l(u)$ 的归一化条件表示为
$$ \int_0^{r_s}u^2|R_l(u;E_\nu)|^2\,\mathrm du=1 \tag{3.554} $$
因此对式（3.554）关于 $E$ 求导，即
$$ \int_0^{r_s}u^2R_l^*(u;E_\nu)\dot R_l(u;E_\nu)\,\mathrm du=0 \quad\text{（实径向函数及固定归一化）} \tag{3.555} $$
即 $R_l(u)$ 和 $\dot{R}_l(u)$ 彼此正交。这个性质表明，可以将相对于给定 $E_\nu$ 所求得的 $R_l(u;E_\nu)$ 以及 $\dot{R}_l(u;E_\nu)$ 同时作为一组基函数，用来展开波函数。这样，LAPW 方法的基函数表示为\cite{marcus1967variational}
$$ \chi(\boldsymbol{K}_i,\boldsymbol{r})=\begin{cases}\displaystyle\sum_{l=0}\sum_{m=-l}^{l}(A_{lm}^{s,\boldsymbol{K}_i}R_l^{s}(u)+B_{lm}^{s,\boldsymbol{K}_i}\dot{R}_l^{s}(u))\mathrm{Y}_l^{m}(\hat{\boldsymbol{u}}),&u\leqslant r_s\\\dfrac{1}{\sqrt{\varOmega}}\mathrm{e}^{\mathrm{i}\boldsymbol{K}_i\cdot\boldsymbol{r}},&u\gt r_s\end{cases} \tag{3.556} $$
至此，基函数不再依赖于 $E$。这里需要强调，虽然 LAPW 方法多用了一组 $\dot{R}_l(u)$，但是这并不意味着基组的数目增大了一倍，因为函数空间由 $\chi$ 张开，所以哈密顿矩阵的维数由 $\chi$ 的数目决定。容易证明，$\dot{R}_l(u,E_\nu)$ 和 $R_l(u,E_\nu)$ 满足如下关系：
$$ (\hat{H}-E_\nu)\dot{R}_l(u)=R_l(u) \tag{3.557} $$
这个关系在后面的计算中将起到重要的作用。
与原始的 APW 方法相同，将II区中的 $\chi(\boldsymbol{K}_i,\boldsymbol{r})$ 按照式（3.526）展开，同时在I区与II区边界处要求函数值连续以及一阶导数值连续，以此来确定I区中的基函数 $\chi_{\mathrm{I}}$。Anderson 所得到的公式形式上比较复杂，所以在这里我们介绍 Koelling 和 Arbman 同样在 1975 年独立得出的结果\cite{koelling1975use}。求解根据两个边界连续性条件得到的关于 $A$、$B$ 的二元一次方程组，可得
$$ \begin{aligned} A_{lm}^{s,\boldsymbol K_i}&=C_{lm}^{s,\boldsymbol K_i} \frac{K_i j_l'(K_ir_s)\dot R_l(r_s)-j_l(K_ir_s)\dot R_l'(r_s)}{D_l},\\ C_{lm}^{s,\boldsymbol K_i}&=\frac{4\pi i^l e^{\mathrm i\boldsymbol K_i\cdot\boldsymbol\tau_s}}{\sqrt\varOmega} Y_l^{m*}(\widehat{\boldsymbol K}_i),\\ D_l&=R_l'(r_s)\dot R_l(r_s)-R_l(r_s)\dot R_l'(r_s). \end{aligned} \tag{3.558} $$
$$ B_{lm}^{s,\boldsymbol K_i}=C_{lm}^{s,\boldsymbol K_i} \frac{j_l(K_ir_s)R_l'(r_s)-K_i j_l'(K_ir_s)R_l(r_s)}{D_l} \tag{3.559} $$
利用关系式
$$ r_s^2\left[R_l'(r_s)\dot R_l(r_s)-R_l(r_s)\dot R_l'(r_s)\right] =2\int_0^{r_s}u^2|R_l(u;E_\nu)|^2\,\mathrm du=2 \tag{3.560} $$
可以进一步简化 $A_{lm}^{s,\boldsymbol{K}_i}$ 和 $B_{lm}^{s,\boldsymbol{K}_i}$。这样，最后得到 LAPW 的基函数
$$ \chi(\boldsymbol K_i,\boldsymbol r)= \begin{cases} \displaystyle\sum_{l=0}^{\infty}\sum_{m=-l}^{l} \frac{2\pi r_s^2i^l e^{\mathrm i\boldsymbol K_i\cdot\boldsymbol\tau_s}}{\sqrt\varOmega} Y_l^{m*}(\widehat{\boldsymbol K}_i)Y_l^m(\hat{\boldsymbol u}) \left[a_l^{s,\boldsymbol K_i}R_l^s(u)+b_l^{s,\boldsymbol K_i}\dot R_l^s(u)\right], &\boldsymbol r=\boldsymbol\tau_s+\boldsymbol u,\ u\le r_s,\\ \displaystyle\varOmega^{-1/2}e^{\mathrm i\boldsymbol K_i\cdot\boldsymbol r},&\boldsymbol r\text{ 在球外}, \end{cases} \tag{3.561} $$
式中
$$ a_l^{s,\boldsymbol K_i}=K_i j_l'(K_ir_s)\dot R_l(r_s)-j_l(K_ir_s)\dot R_l'(r_s) \tag{3.562} $$
$$ b_l^{s,\boldsymbol K_i}=j_l(K_ir_s)R_l'(r_s)-K_i j_l'(K_ir_s)R_l(r_s) \tag{3.563} $$
与原始的 APW 方法计算过程相似，分别计算I区和II区的贡献，则得 LAPW 方法的哈密顿矩阵元 $H_{ij}$ 为
$$ H_{ij}=\int_{\varOmega}\chi_i^*(\boldsymbol r) \left[-\frac12\nabla^2+V(\boldsymbol r)\right]\chi_j(\boldsymbol r)\,\mathrm d^3\boldsymbol r, \qquad\chi_i=\chi(\boldsymbol K_i,\boldsymbol r) \tag{3.564} $$
重叠矩阵元 $\Delta_{ij}$ 为
$$ S_{ij}=\int_{\varOmega}\chi_i^*(\boldsymbol r)\chi_j(\boldsymbol r)\,\mathrm d^3\boldsymbol r \tag{3.565} $$
式中的 $\dot R_l$ 并不自动归一；其模应由径向积分直接计算：
$$ \langle\dot R_l|\dot R_l\rangle_s =\int_0^{r_s}u^2|\dot R_l(u;E_\nu)|^2\,\mathrm du\ge0 \tag{3.566} $$
若基函数在球面连续且一阶可微，分部积分给出显式厄米的等价形式：
$$ H_{ij}=\frac12\int_{\varOmega}\nabla\chi_i^*\cdot\nabla\chi_j\,\mathrm d^3\boldsymbol r +\int_{\varOmega}\chi_i^*V\chi_j\,\mathrm d^3\boldsymbol r \quad\text{（连续可微、周期边界）} \tag{3.567} $$
于是有
$$ H_{ji}=H_{ij}^{*},\qquad S_{ji}=S_{ij}^{*} \tag{3.568} $$
这样，LAPW 方法的久期方程即为标准的广义本征值问题：
$$ \det|\,H_{ij}-E\Delta_{ij}\,|=0 \tag{3.569} $$
与原始 APW 方法的久期方程（式（3.538））相比，LAPW 方法的久期方程中表面积分不见了，而在 $H_{ij}$ 中多了 $\dot{R}_l(u;E_\nu)$ 的贡献。
此外，我们在这里不加证明地给出 Anderson 得到的I区中的基函数：
$$ \chi_{\mathrm{I}}(\boldsymbol{K}_i,\boldsymbol{r})=\sum_{s}\frac{4\pi r_s^{2}}{\sqrt{\varOmega}}\mathrm{e}^{\mathrm{i}\boldsymbol{K}_i\cdot\boldsymbol{\tau}_s}\sum_{l=0}\sum_{m=-l}^{l}i^{l}\mathrm{j}_l(K_ir_s)\mathrm{Y}_l^{m*}(\hat{\boldsymbol{K}}_i)\mathrm{Y}_l^{m}(\hat{\boldsymbol{u}})\frac{\varPhi_l^{\nu}(\tilde{D}_{\boldsymbol{K}_i},u)}{\varPhi_l^{\nu}(\tilde{D}_{l,\boldsymbol{K}_i},r_s)} \tag{3.570} $$
式中
$$ \varPhi_l^{\nu}(\tilde{D}_{\boldsymbol{K}_i},u)=R_l(u;E_\nu)+\omega(\tilde{D}_{l,\boldsymbol{K}_i})\dot{R}_l(u;E_\nu) \tag{3.571} $$
而 $D=\dfrac{x\,\partial\ln\varphi(x)}{\partial(x)}$。式（3.571）中出现的变量具体如下：
$$ \tilde{D}_{l,\boldsymbol{K}_i}=\frac{K_ir_s}{\mathrm{j}_l(K_ir_s)}\left.\frac{\partial\mathrm{j}_l(x)}{\partial x}\right|_{x=K_ir_s} \tag{3.572} $$
$$ \begin{gathered} \omega(\tilde{D}_{l,\boldsymbol{K}_i})=-\frac{R_l(r_s;E_\nu)}{\dot{R}_l(r_s;E_\nu)}\frac{\tilde{D}_{\boldsymbol{K}_i}-D_l^{\nu}}{\tilde{D}_{l,\boldsymbol{K}_i}-D_l^{\nu}}\\ D_l^{\nu}=r_sR_l'(r_s)/R_l(r_s),\quad D_l^{\nu}=r_s\dot{R}_l'(r_s;E_\nu)/\dot{R}_l(r_s;E_\nu) \end{gathered} $$
不难证明，Koelling 和 Arbman 得到的I区基函数与 Anderson 的结果等价。
关于势函数的讨论
从 3.5.1 节的讨论可知，MT 势的合适选取（或称构建）决定了 APW 方法的计算结果。需要指出，I区中的 MT 势不能简单理解为离子势，它是包含离子-电子相互作用、电子-电子库仑相互作用以及电子-电子多体作用的等效势函数。因此，由正确的 $V(\boldsymbol{r})$ 得到的结果与由第 3 章中介绍的 Hartree-Fock 方法及密度泛函理论得到的结果相同，或至少非常接近。早期的工作中通常先利用 Hartree-Fock 方法对各原子特定的电子组态进行自洽场计算，得到孤立的原子轨道以及原子电荷分布。例如第 $s$ 个原子的电荷为
$$ \rho_0^{s}(r)=\sum_{\mathrm{occ}}|\,\phi_{lm}^{s}(\boldsymbol{r})\,|^{2} \tag{3.573} $$
其中求和遍历该原子所有占据轨道，且 $\rho(r)$ 满足球对称条件。然后利用泊松方程
$$ \boldsymbol{\nabla}^{2}V_{\mathrm{H}}^{s}(r)=-4\pi\rho_0^{s}(r) \tag{3.574} $$
求出电子-电子库仑势，从而得到第 $s$ 个原子的库仑势
$$ V_{\mathrm{Coul}}^{s}=-\frac{Z_s}{r}+V_{\mathrm{H}}^{s}(r) \tag{3.575} $$
孤立原子组成体系时，需要额外考虑能级的对齐，这反映在邻近原子对 $V_{\mathrm{Coul}}^{s}$ 的修正上。Ern 与 Switendick、Scop 均指出，该修正可以由 Madelung 常数 $\alpha$ 确定\cite{ern1965electronic,scop1965band}，即
$$ V_{\mathrm{M}}=-4\alpha/a_0 \tag{3.576} $$
然后从阳离子的 $V_{\mathrm{Coul}}$ 中扣除 $V_{\mathrm{M}}$，而在阴离子的 $V_{\mathrm{Coul}}$ 中加入 $V_{\mathrm{M}}$。
此外，还应该考虑电子-电子的交换势：
$$ V_{\mathrm{x}}^{s}(r)=-6[3\rho^{s}(r)/(8\pi)]^{1/3} \tag{3.577} $$
式中：$\rho^{s}(r)=\rho_0^{s}(r)+\displaystyle\sum_{j}\rho_0^{j}(r_{sj})$，即计入了邻近原子电荷密度的贡献。
因此，实际计算中，MT 势应为
$$ V(u)=\begin{cases}V_{\mathrm{Coul}}^{s}(u)\pm V_{\mathrm{M}}(u)+V_{\mathrm{x}}^{s}(u)-V_0,&u\leqslant r_s\\0,&u\gt r_s\end{cases} \tag{3.578} $$
式中：$V_0$ 是一个可调参数，目的是将II区变为 $V(\boldsymbol{r})\equiv0$ 的自由空间，可以根据实验值确定，也可通过对II区的 $V(\boldsymbol{r})$ 求平均值 $\bar{V}_{\mathrm{II}}$ 求得。关于 MT 势的确定还有其他一些处理办法，这里不再详述，请参看文献\cite{loucks1967augmented}。
密度泛函理论提出之后，可以严格地通过电荷的空间分布确定体系的多体作用势，这使得我们可以通过自洽场计算，不断地更新电荷密度而对体系进行求解。更为精确的计算可以通过取消对 MT 势的形状假设， 即允许球内势含非球对称分量、间隙区势随位置变化，并在两区连续地表示体系的有效势。这种计算方法称为全势线性缀加平面波（full-potential LAPW，FLAPW）方法， 它是常用的高精度全电子计算方法之一。