本文是「第一性原理的微观计算模拟」系列的第 10 篇（共 11 篇），内容整理自同名书稿第 3 章，公式编号与原书一致。文中引用的文献依据原书参考文献清单整理并列于文末，按原书清单顺序从 1 开始编号。\nocite{*}
激发态与多体理论
基于 Kohn-Sham 方程的局域密度近似和广义梯度近似的交换关联泛函在材料模拟领域内的应用获得了巨大的成功。然而，与此同时，它们也存在很大局限性，最受诟病的就是无法给出准确的电子激发能，即严重低估了半导体和绝缘体材料的能隙。考虑一个单电子激发的过程，例如能带 $i$ 中的一个电子受能量为 $\hbar\omega$ 的光子激发成为自由电子，或一个入射的自由电子被能带 $j$ 捕获，同时释放能量为 $\hbar\omega$ 的光子，由于体系受到扰动（电子数发生变化），所以体系必然会做出响应：当受激电子在体系中运动时，周围电荷密度会因为响应而发生变化。此时的受激电子相当于携带着一部分正电荷在体系中运动，因此电子-电子相互作用是动态变化的，与时间相关。这种动态效果无法用 LDA 和 GGA 泛函正确地描述。为了获得一个简明的物理图像，一般视这个携带着响应的受激电子为一个准粒子，而利用多体理论进行求解。
格林函数理论与 Dyson 方程
格林函数在凝聚态理论领域有着非常重要的位置。原则上，知道一个体系的格林函数，就可以得到大多数我们感兴趣的性质，如态密度、电荷密度分布、本征能级等。在零温下单粒子格林函数定义为
$$ \begin{aligned} G(\boldsymbol{r}t,\boldsymbol{r}'t')&=-\mathrm{i}\langle N\,|\,T[\psi(\boldsymbol{r}t)\psi^{\dagger}(\boldsymbol{r}'t')]\,|\,N\rangle\\ &=\begin{cases}-\mathrm{i}\langle N\,|\,\psi(\boldsymbol{r}t)\psi^{\dagger}(\boldsymbol{r}'t')\,|\,N\rangle,&t\gt t'\\\mathrm{i}\langle N\,|\,\psi^{\dagger}(\boldsymbol{r}'t')\psi(\boldsymbol{r}t)\,|\,N\rangle,&t'\geqslant t\end{cases} \end{aligned} \tag{3.602} $$
式中：$|\,N\rangle$ 是 $N$ 电子体系的基态；变量 $\boldsymbol{r}$ 同时指代电子的位置 $\boldsymbol{r}$ 和自旋态 $\sigma$；$\psi(\boldsymbol{r}t)$ 是海森堡绘景下的场算符，例如，$\psi^{\dagger}(\boldsymbol{r}t)\,|\,N\rangle$ 表示在时刻 $t$ 将一个自旋为 $\sigma$ 的电子加在 $\boldsymbol{r}$ 的终点处，使得体系含有 $N+1$ 个电子；$T$ 是时序算符，它保证更新的时刻总在左侧。若 $t\gt t'$，式（3.602）给出的是在 $t'$ 时刻将一个电子加入体系 $\boldsymbol{r}'$ 之后，$t$ 时刻在 $\boldsymbol{r}$ 的终点处观测到的一个电子的概率振幅；若 $t'\geqslant t$，则式（3.602）给出的是在 $t$ 时刻将一个电子于 $\boldsymbol{r}$ 的终点处取走后，$t'$ 时刻在 $\boldsymbol{r}'$ 的终点处观测到的一个空位的概率振幅。显然，准粒子格林函数描述的正是我们在 3.7.1 节中介绍的单电子激发过程。
由式（3.602），我们还可以得到格林函数 $G$ 的谱表示。在海森堡绘景下场算符 $\psi(\boldsymbol{r}t)$ 可表示为
$$ \psi(\boldsymbol{r}t)=\mathrm{e}^{\mathrm{i}\hat{H}t}\psi(\boldsymbol{r})\mathrm{e}^{-\mathrm{i}\hat{H}t} \tag{3.603} $$
式中：$\hat{H}$ 为多体哈密顿算符。可以将一组 $N+1$ 或 $N-1$ 电子体系完备基 $|\,N\pm1\rangle$ 根据 $t$ 与 $t'$ 的时序关系插入式（3.602）。举例来说，当 $t\gt t'$ 时，有
$$ \begin{aligned} G(\boldsymbol r t,\boldsymbol r't') &=-\mathrm i\langle N|\psi(\boldsymbol r,t)\psi^\dagger(\boldsymbol r',t')|N\rangle\\ &=-\mathrm i\sum_s\langle N|\psi(\boldsymbol r)|N+1,s\rangle e^{-\mathrm i(E_{N+1,s}-E_N)(t-t')} \langle N+1,s|\psi^\dagger(\boldsymbol r')|N\rangle, \qquad t\gt t'. \end{aligned} \tag{3.604} $$
对于 $t'\gt t$ 时的空穴运动，可类似地得到
$$ G(\boldsymbol r t,\boldsymbol r't')=\mathrm i\sum_s \langle N|\psi^\dagger(\boldsymbol r')|N-1,s\rangle e^{-\mathrm i(E_N-E_{N-1,s})(t-t')} \langle N-1,s|\psi(\boldsymbol r)|N\rangle,\qquad t\lt t'. \tag{3.605} $$
再对其做傅里叶变换，得\cite{onida2002electronic}
$$ G(\boldsymbol{r},\boldsymbol{r}',\omega)=\sum_{s}\frac{f_s(\boldsymbol{r})f_s^{*}(\boldsymbol{r}')}{\omega-\varepsilon_s+\mathrm{i}\eta\,\mathrm{sgn}(\varepsilon_s-\mu)} \tag{3.606} $$
式中，$\eta$ 为正的无穷小量；$\mu$ 为化学势；$\mathrm{sgn}(\varepsilon_s-\mu)$ 代表 $\varepsilon_s-\mu$ 的符号，当 $\varepsilon_s\gt \mu$ 时，$\varepsilon_s=E_{N+1,s}-E_N$，当 $\varepsilon_s\lt \mu$ 时，$\varepsilon_s=E_N-E_{N-1,s}$；$\omega$ 是具有能量的量纲；另有
$$ f_s(\boldsymbol{r})=\begin{cases}\langle N\,|\,\psi(\boldsymbol{r})\,|\,N+1,s\rangle,&\varepsilon_s\gt \mu\\\langle N-1,s\,|\,\psi(\boldsymbol{r})\,|\,N\rangle,&\varepsilon_s\leqslant\mu\end{cases} \tag{3.607} $$
式（3.606）和式（3.607）就是所求的谱表示。还可以据此得到谱函数 $A(\boldsymbol{r},\boldsymbol{r}',\omega)$，即
$$ A(\boldsymbol r,\boldsymbol r',\omega) =\frac{\mathrm i}{2\pi}\left[G^{\mathrm R}(\boldsymbol r,\boldsymbol r',\omega) -G^{\mathrm A}(\boldsymbol r,\boldsymbol r',\omega)\right] =\sum_s f_s(\boldsymbol r)f_s^*(\boldsymbol r')\delta(\omega-\varepsilon_s) \quad\text{（离散谱）} \tag{3.608} $$
现在需要进一步推导 $G$ 所遵循的运动方程。体系的多体哈密顿量在粒子数空间中表示为
$$ \begin{aligned} \hat H={}&\int\mathrm d\boldsymbol r\, \psi^\dagger(\boldsymbol r,t)\hat h_0(\boldsymbol r)\psi(\boldsymbol r,t)\\ &+\frac12\iint\mathrm d\boldsymbol r\,\mathrm d\boldsymbol r'\, \psi^\dagger(\boldsymbol r,t)\psi^\dagger(\boldsymbol r',t) v(\boldsymbol r,\boldsymbol r')\psi(\boldsymbol r',t)\psi(\boldsymbol r,t) \end{aligned} \tag{3.609} $$
式中：$\hat{H}_0$ 为单体算符，${\hat{H}_0=-\boldsymbol{\nabla}^{2}/2+V_{\mathrm{ext}}}$。场算符依据海森堡方程进行演化，即
$$ \mathrm{i}\frac{\partial\psi(\boldsymbol{r},t)}{\partial t}=[\psi(\boldsymbol{r},t),\hat{H}] \tag{3.610} $$
将式（3.609）代入式（3.610），可得
$$ \begin{aligned} &\bigl(\mathrm i\partial_t-\hat h_0(\boldsymbol r)\bigr) G(\boldsymbol r t,\boldsymbol r't')\\ &\quad+\mathrm i\int\!\mathrm d\boldsymbol r''\; v(\boldsymbol r,\boldsymbol r'') \langle N|T[\psi^\dagger(\boldsymbol r''t)\psi(\boldsymbol r''t) \psi(\boldsymbol r t)\psi^\dagger(\boldsymbol r't')]|N\rangle =\delta(\boldsymbol r-\boldsymbol r')\delta(t-t'). \end{aligned} \tag{3.611} $$
$\mathrm{i}^{2}\langle N\,|\,T[\psi^{\dagger}(\boldsymbol{r}''t)\psi(\boldsymbol{r}''t)\psi(\boldsymbol{r}t)\psi^{\dagger}(\boldsymbol{r}'t')]\,|\,N\rangle$ 正是二体格林函数 $G_2(1,3,2,3^{\dagger})$ 的定义。由方程（3.611）可以看到，想求得 $G$，需要先知道 $G_2$，而求解 $G_2$ 又必须知道 $G_3$，依此类推。这就是所谓的 BBGKY 级列（BBGKY hierarchy）。为了避免这个困难，用自能算符 $\varSigma$ 代替二体格林函数，则有\cite{aryasetiawan1998method}
$$ \begin{aligned} &\left(\mathrm i\partial_t-\hat h_0(\boldsymbol r)-V_{\mathrm H}(\boldsymbol r)\right) G(\boldsymbol r t,\boldsymbol r't')\\ &\quad-\int\!\mathrm dt''\,\mathrm d\boldsymbol r''\; \Sigma(\boldsymbol r t,\boldsymbol r''t'') G(\boldsymbol r''t'',\boldsymbol r't') =\delta(\boldsymbol r-\boldsymbol r')\delta(t-t'). \end{aligned} \tag{3.612} $$
这就是所求的 Dyson 方程。需要注意，与式（3.611）左端第二项相比，方程（3.612）中的自能项已经抛除了 Hartree 项的贡献。对式（3.612）做傅里叶变换，即可得到 Dyson 方程在频率空间下的形式
$$ (\omega-\hat{H}_0-V_{\mathrm{H}})G(\boldsymbol{r},\boldsymbol{r}',\omega)-\int\mathrm{d}\boldsymbol{r}''\,\varSigma(\boldsymbol{r},\boldsymbol{r}'',\omega)G(\boldsymbol{r}'',\boldsymbol{r}',\omega)=\delta(\boldsymbol{r}-\boldsymbol{r}') \tag{3.613} $$
如果设 $G_0$ 为 $\varSigma=0$ 时对应的格林函数，则可以将式（3.613）重新写为
$$ G(12)=G_0(12)+\int G_0(13)\varSigma(34)G(42)\,\mathrm{d}(34) \tag{3.614} $$
其中，数字“1”代表位置、时间和自旋态 $\{\boldsymbol{r}_1,t_1,\sigma_1\}$。显然，$G_0(12)$ 取决于不同的参考体系。在实践中参考体系往往选取利用标准 Kohn-Sham 方程求得的基态，因此一般取
$$ \hat{H}_0=-\boldsymbol{\nabla}^{2}/2+V_{\mathrm{ext}}+V_{\mathrm{xc}} \tag{3.615} $$
此时 $G_0(12)$ 相对应的自能算符为 $\Delta\varSigma=\varSigma-V_{\mathrm{xc}}$。
GW 方法
Hedin 方程
目前自能项的计算都基于 Hedin 于 1965 年提出的方程组\cite{hedin1965method}：
$$ \begin{cases} \Sigma(1,2)=\mathrm i\!\displaystyle\int\!G(1,3)W(1^+,4)\Gamma(3,2,4)\,\mathrm d(34),\\ G(1,2)=G_{\mathrm H}(1,2)+\displaystyle\int\!G_{\mathrm H}(1,3)\Sigma(3,4)G(4,2)\,\mathrm d(34),\\ P(1,2)=-\mathrm i\!\displaystyle\int\!G(1,3)G(4,1)\Gamma(3,4,2)\,\mathrm d(34),\\ W(1,2)=v(1,2)+\displaystyle\int\!v(1,3)P(3,4)W(4,2)\,\mathrm d(34),\\ \Gamma(1,2,3)=\delta(1,2)\delta(1,3) +\displaystyle\int\!\frac{\delta\Sigma(1,2)}{\delta G(4,5)} G(4,6)G(7,5)\Gamma(6,7,3)\,\mathrm d(4567). \end{cases} \tag{3.616} $$
式（3.616）中的 $G_{\mathrm H}$ 取仅含外势和 Hartree 势的参考传播子；若改用 Kohn–Sham 传播子，应在 Dyson 方程中将 $\Sigma$ 换成 $\Sigma-V_{\mathrm{xc}}$。该方程组的计算就是常常遇到的所谓 Hedin 五边形（见图 3.19）问题，其中 $\varGamma$ 称为顶点函数（vertex function）。方程组（3.616）中的五个方程原则上可以迭代求解，但是实际操作中无法求得 $G$ 的精确值，而且计算量很大，步骤也比较繁难，所以目前该方程组求解的实现都是基于无规相近似（RPA）所得到的简化模型。
图 3.19 Hedin 五边形
GW 近似
利用 RPA 方法可以极大地简化 Hedin 五边形。RPA 相当于忽略顶点修正，即
$$ \varGamma(123)\approx\delta(12)\delta(13) \tag{3.617} $$
于是有
$$ \begin{cases} \Sigma(1,2)=\mathrm iG(1,2)W(1^+,2),\\ G=G_{\mathrm H}+G_{\mathrm H}\Sigma G,\\ P(1,2)=-\mathrm iG(1,2)G(2,1^+),\\ W=v+vPW. \end{cases} \tag{3.618} $$
即在 RPA 下，可以将自能项近似地表示为准粒子格林函数 $G$ 及屏蔽库仑势 $W$ 的卷积。而且忽略顶点修正之后，剩下的四个方程组成了一个封闭的方程组，可大大简化 $G$ 的求解过程。方程组（3.618）最后一个关于 $W$ 的计算式比较复杂，这里给出另一种表示方法：
$$ \begin{cases} W(12)=\displaystyle\int\mathrm{d}(3)\,\epsilon^{-1}(13)v(32)\\ \epsilon(12)=\delta(12)-\displaystyle\int\mathrm{d}(3)\,v(1,3)P(3,2) \end{cases} \tag{3.619} $$
设自能算符 $\varSigma$ 已知，则将格林函数的谱表示代入方程（3.613），可得
$$ (H_0+V_{\mathrm{H}})\psi_s(\boldsymbol{r})+\int\varSigma(\boldsymbol{r},\boldsymbol{r}',\omega)\psi_s(\boldsymbol{r}')\,\mathrm{d}\boldsymbol{r}'=\varepsilon_s\psi_s \tag{3.620} $$
可以看到，该方程与常用的 Kohn-Sham 方程极为类似，二者唯一的区别在于这里用自能算符 $\varSigma$ 代替了 Kohn-Sham 方程中的交换关联势 $V_{\mathrm{xc}}$。因为 $V_{\mathrm{xc}}$ 并没有包含动态介电函数 $\epsilon$，所以是一个静态的作用，意味着第 $i$ 个态上电子的激发不会引起其他电子状态的改变，即没有计入体系对外势场的响应。而 GW 近似中的自能项 $\varSigma$ 则通过 $\epsilon$ 考虑了这种响应。因此，准粒子近似（QPA）可以改进 LDA 或 GGA 泛函的事实也可以理解为，QPA 的自洽场方程中所用的自能项 $\varSigma$ 可以比 $V_{\mathrm{xc}}$ 更好地描述单粒子激发的能量变化。虽然已经做了上述简化，但是到目前为止仍然不知道初始的格林函数 $G$ 应该如何选取。事实上，上述方程的解对初始值比较敏感，所以必须选取一套比较合理的迭代初始值。Hybertsen 和 Louie 开创性地提出，将用标准 Kohn-Sham 方法得到的一套本征波函数 $\{\psi_s^{\mathrm{KS}}\}$ 作为 $G$ 谱表示方程（3.606）中的 $f_s(\boldsymbol{r})$\cite{hybertsen1986electron}。
对方程组（3.618）中的第一个方程做傅里叶变换，可得
$$ \varSigma(\boldsymbol r,\boldsymbol r',\omega)=\frac{\mathrm i}{2\pi} \int_{-\infty}^{\infty}\mathrm d\omega'\, e^{\mathrm i\omega'0^+}G(\boldsymbol r,\boldsymbol r',\omega+\omega') W(\boldsymbol r,\boldsymbol r',\omega') \tag{3.621} $$
式中
$$ G(\boldsymbol{r},\boldsymbol{r}',\omega+\omega')=\sum_{s}\frac{\psi_s^{\mathrm{KS}}(\boldsymbol{r})\psi_s^{*\,\mathrm{KS}}(\boldsymbol{r}')}{\omega+\omega'-\varepsilon_s+\mathrm{i}\eta\,\mathrm{sgn}(\varepsilon_s-\mu)} \tag{3.622} $$
$$ W(\boldsymbol r,\boldsymbol r',\omega') =\int\mathrm d\boldsymbol r''\,\epsilon^{-1}(\boldsymbol r,\boldsymbol r'',\omega')v(\boldsymbol r'',\boldsymbol r') \tag{3.623} $$
需要注意的是，严格来说，由于参考态中已经包含了多体项 $V_{\mathrm{xc}}$ 的贡献，所以此时的自能项应为 3.7.2 节最后给出的 $\Delta\varSigma$。在不引起误解的前提下，为了使形式简洁，在此后的讨论中，我们仍然用 $\varSigma$ 表示自能项。式（3.622）给出的仅仅是 $G(\boldsymbol{r},\boldsymbol{r}',\omega)$ 的一个初始猜测。严格的 GW 近似要求从式（3.622）出发迭代地求解 $G$ 和 $W$。为了简化计算，减少计算时间，也可以在整个计算过程中不更新这两个量，这样的简化称为 $G_0W_0$ 近似。将式（3.615）和式（3.621）代入方程（3.620），可得 $G_0W_0$ 近似下的准粒子能量 $\varepsilon_s^{\mathrm{QP}}$，且
$$ \varepsilon_s^{\mathrm{QP}}=\varepsilon_s^{\mathrm{KS}}+\langle\psi_s^{\mathrm{KS}}\,|\,\varSigma(\varepsilon_s^{\mathrm{QP}})-V_{\mathrm{xc}}\,|\,\psi_s^{\mathrm{KS}}\rangle \tag{3.624} $$
因为 $\varSigma$ 是待求本征能量 $\varepsilon_s^{\mathrm{QP}}$ 的函数，所以必须用迭代办法求解上述方程。一般情况下采取线性化手段避免迭代过程。这里给出最后的结果，即
$$ \varepsilon_s^{\mathrm{QP}}\simeq\varepsilon_s^{\mathrm{KS}}+Z_s \langle\psi_s^{\mathrm{KS}}|\varSigma(\varepsilon_s^{\mathrm{KS}})-V_{\mathrm{xc}}|\psi_s^{\mathrm{KS}}\rangle \tag{3.625} $$
式中
$$ Z_s=\left[1-\mathrm{Re}\langle\psi_s^{\mathrm{KS}}\,|\,\frac{\partial\varSigma(\omega)}{\partial\omega}\,|\,\varepsilon_s^{\mathrm{KS}}\,|\,\psi_s^{\mathrm{KS}}\rangle\right]^{-1} \tag{3.626} $$
平面波基框架下 GW 计算的实现
在前面的讨论中已经知道了 GW 近似（或至少是 $G_0W_0$ 近似）的基本公式。因为一般选择 Kohn-Sham 方程所得到的本征谱作为参考态，所以在 GW 近似具体的算法实现中，仍然需要选择最方便的基函数。使用平面波基可以得到比较简单、直接的计算公式，同时辅以成熟的快速傅里叶变换技术，在程序实现上有着极大的便利。首先给出几个关键物理量的表达式\cite{shishkin2006implementation}：
$$ W_{\boldsymbol{q}}(\boldsymbol{G},\boldsymbol{G}',\omega)=4\pi e^{2}\frac{1}{|\,\boldsymbol{q}+\boldsymbol{G}\,|}\epsilon_{\boldsymbol{q}}^{-1}(\boldsymbol{G},\boldsymbol{G}',\omega)\frac{1}{|\,\boldsymbol{q}+\boldsymbol{G}'\,|} \tag{3.627} $$
$$ \epsilon_{\boldsymbol{q}}(\boldsymbol{G},\boldsymbol{G}',\omega)=\delta_{\boldsymbol{G}\boldsymbol{G}'}-\frac{4\pi e^{2}}{|\,\boldsymbol{q}+\boldsymbol{G}\,|\,|\,\boldsymbol{q}+\boldsymbol{G}'\,|}\chi_{\boldsymbol{q}}^{0}(\boldsymbol{G},\boldsymbol{G}',\omega) \tag{3.628} $$
$$ \chi_{\boldsymbol{q}}^{0}(\boldsymbol{G},\boldsymbol{G}',\omega)=\frac{1}{\varOmega_{\mathrm{cell}}}\sum_{nn'\boldsymbol{k}}2w_{\boldsymbol{k}}(f_{n'\boldsymbol{k}-\boldsymbol{q}}-f_{n\boldsymbol{k}})\times\frac{\langle\psi_{n'\boldsymbol{k}-\boldsymbol{q}}\,|\,\mathrm{e}^{-\mathrm{i}(\boldsymbol{q}+\boldsymbol{G})\boldsymbol{r}}\,|\,\psi_{n\boldsymbol{k}}\rangle\langle\psi_{n\boldsymbol{k}}\,|\,\mathrm{e}^{\mathrm{i}(\boldsymbol{q}+\boldsymbol{G}')\boldsymbol{r}'}\,|\,\psi_{n'\boldsymbol{k}-\boldsymbol{q}}\rangle}{\omega+\varepsilon_{n'\boldsymbol{k}-\boldsymbol{q}}-\varepsilon_{n\boldsymbol{k}}+\mathrm{i}\eta\,\mathrm{sgn}(\varepsilon_{n'\boldsymbol{k}-\boldsymbol{q}}-\varepsilon_{n\boldsymbol{k}})} \tag{3.629} $$
式中：$w_{\boldsymbol{k}}$ 是第一布里渊区内 $\boldsymbol{k}$ 点的权重；$\langle\psi_{n'\boldsymbol{k}-\boldsymbol{q}}\,|\,\mathrm{e}^{-\mathrm{i}(\boldsymbol{q}+\boldsymbol{G})\boldsymbol{r}}\,|\,\psi_{n\boldsymbol{k}}\rangle$ 称为交换电荷密度。
在具体实现中，为了使函数有良好的行为，一般对式（3.627）扣除库仑势 $V_{\boldsymbol{q}}$，有
$$ \overline{W}_{\boldsymbol{q}}(\boldsymbol{G},\boldsymbol{G}',\omega)=4\pi e^{2}\frac{1}{|\,\boldsymbol{q}+\boldsymbol{G}\,|}[\epsilon_{\boldsymbol{q}}^{-1}(\boldsymbol{G},\boldsymbol{G}',\omega)-\delta_{\boldsymbol{G}\boldsymbol{G}'}]\frac{1}{|\,\boldsymbol{q}+\boldsymbol{G}'\,|} \tag{3.630} $$
而相应的自能项 $\bar{\varSigma}$ 为\cite{shishkin2006implementation}
$$ \begin{aligned} \bar{\varSigma}(\omega)_{n\boldsymbol{k},n\boldsymbol{k}}={}&\frac{1}{\varOmega_{\mathrm{cell}}}\sum_{\boldsymbol{q}\boldsymbol{G},\boldsymbol{G}'}\sum_{n'}\frac{\mathrm{i}}{2\pi}\int_{0}^{\infty}\mathrm{d}\omega'\,\overline{W}_{\boldsymbol{q}}(\boldsymbol{G},\boldsymbol{G}',\omega)\times\langle\psi_{n'\boldsymbol{k}-\boldsymbol{q}}\,|\,\mathrm{e}^{\mathrm{i}(\boldsymbol{q}+\boldsymbol{G})\boldsymbol{r}}\,|\,\psi_{n\boldsymbol{k}}\rangle\langle\psi_{n\boldsymbol{k}}\,|\,\mathrm{e}^{-\mathrm{i}(\boldsymbol{q}+\boldsymbol{G}')\boldsymbol{r}'}\,|\,\psi_{n'\boldsymbol{k}-\boldsymbol{q}}\rangle\\ &\times\left(\frac{1}{\omega+\omega'-\varepsilon_{n'\boldsymbol{k}-\boldsymbol{q}}+\mathrm{i}\eta\,\mathrm{sgn}(\varepsilon_{n'\boldsymbol{k}-\boldsymbol{q}}-\mu)}+\frac{1}{\omega-\omega'-\varepsilon_{n'\boldsymbol{k}-\boldsymbol{q}}+\mathrm{i}\eta\,\mathrm{sgn}(\varepsilon_{n'\boldsymbol{k}-\boldsymbol{q}}-\mu)}\right) \end{aligned} \tag{3.631} $$
最后，在计算能量的时候需要加入交换项的贡献，即
$$ \varepsilon_s^{\mathrm{QP}}=\varepsilon_s^{\mathrm{KS}}+Z_s\langle\psi_s^{\mathrm{KS}}\,|\,\bar{\varSigma}(\varepsilon_s^{\mathrm{KS}})-V_{\mathrm{xc}}+V_x\,|\,\psi_s^{\mathrm{KS}}\rangle \tag{3.632} $$
从原则上讲，通过式（3.628）至式（3.632），已经可以实现平面波基框架下的 GW 计算。但是在具体实现中仍然需要注意下面两个问题。
1. 交换电荷
交换电荷密度矩阵并不是在所有的 GW 实现中都需要特别注意。但是平面波基组总是与赝势方法联系在一起，所以实际上参与构建格林函数 $G$ 以及自能项 $\varSigma$ 的均为相应的赝波函数，而且没有芯区的电子轨道信息。但是方程（3.631）要求参与计算的是精确的波函数，因此只有采用 PAW 方法才可以满足这个条件。对 PAW 方法的详细分析已经超出了本书的范围。这里我们只给出几个重要的公式，更详细的讨论可参阅文献\cite{shishkin2006implementation,gajdos2006linear}。PAW 方法将准确的单电子波函数表示为
$$ |\,\psi_{n\boldsymbol{k}}\rangle=|\,\tilde{\psi}_{n\boldsymbol{k}}\rangle+\sum_{i}(|\,\phi_i\rangle-|\,\tilde{\phi}_i\rangle)\langle\tilde{p}_i\,|\,\tilde{\psi}_{n\boldsymbol{k}}\rangle \tag{3.633} $$
式中：$|\,\tilde{\psi}_{n\boldsymbol{k}}\rangle$ 为赝波函数；$|\,\phi_i\rangle$ 为给定原点 $\boldsymbol{R}_i$、角动量 $l_i$ 以及非自旋极化参考能 $\epsilon_i$ 的径向薛定谔方程的全电子径向解；$|\,\tilde{\phi}_i\rangle$ 为相应的径向解；$|\,\tilde{p}_i\rangle$ 为投影算符，且满足
$$ \langle\tilde{p}_i\,|\,\tilde{\phi}_j\rangle=\delta_{ij} \tag{3.634} $$
赝波函数 $|\,\tilde{\psi}_{n\boldsymbol{k}}\rangle$ 和投影算符 $|\,\tilde{p}_i\rangle$ 可以表示成 Bloch 波的形式，即
$$ |\,\tilde{\psi}_{n\boldsymbol{k}}\rangle=\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}}\,|\,\tilde{u}_{n\boldsymbol{k}}\rangle,\quad|\,\tilde{p}_{n\boldsymbol{k}}\rangle=\mathrm{e}^{-\mathrm{i}\boldsymbol{k}(\boldsymbol{r}-\boldsymbol{R}_i)}\,|\,\tilde{p}_i\rangle $$
因为有 $|\,\psi_{n\boldsymbol{k}}\rangle$、$|\,\tilde{p}_i\rangle$ 等项的存在，PAW 方法中电荷密度矩阵元 $\rho_{nn'}$ 有比较复杂的形式。为简单起见，暂不考虑对 $\boldsymbol{k}$ 的求和，即第一布里渊区内只有 $\varGamma$ 点，则 $\rho_{nn'}$ 可表示为\cite{paier2005perdew}
$$ \begin{aligned} \rho_{nn'}(\boldsymbol{r})&=\tilde{\rho}_{nn'}(\boldsymbol{r})-\tilde{\rho}_{nn'}^{1}(\boldsymbol{r})+\rho_{nn'}^{1}(\boldsymbol{r})\\ &=\langle\tilde{\psi}_n\,|\,\boldsymbol{r}\rangle\langle\boldsymbol{r}\,|\,\tilde{\psi}_{n'}\rangle-\sum_{ij}\langle\tilde{\phi}_i\,|\,\boldsymbol{r}\rangle\langle\boldsymbol{r}\,|\,\tilde{\phi}_j\rangle\langle\tilde{\psi}_n\,|\,\tilde{p}_i\rangle\langle\tilde{p}_j\,|\,\tilde{\psi}_{n'}\rangle\\ &\quad+\sum_{ij}\langle\phi_i\,|\,\boldsymbol{r}\rangle\langle\boldsymbol{r}\,|\,\phi_j\rangle\langle\tilde{\psi}_n\,|\,\tilde{p}_i\rangle\langle\tilde{p}_j\,|\,\tilde{\psi}_{n'}\rangle \end{aligned} \tag{3.635} $$
式中：$\tilde{\rho}_{nn'}(\boldsymbol{r})$ 为赝波函数的贡献；$\tilde{\rho}_{nn'}^{1}(\boldsymbol{r})$ 和 $\rho_{nn'}^{1}(\boldsymbol{r})$ 分别为芯区全电子补偿项以及赝势补偿项（上标“1”代表单中心局域量），这两种补偿类似于超软赝势（USPP）中为修正芯区电子总数不足而做的补偿。显然，由于补偿项的存在，方程（3.629）和方程（3.631）中的交换电荷密度矩阵也有比较复杂的形式\cite{shishkin2006implementation}：
$$ \begin{aligned} \langle\psi_{n'\boldsymbol{k}-\boldsymbol{q}}\,|\,\mathrm{e}^{-\mathrm{i}(\boldsymbol{q}+\boldsymbol{G})\boldsymbol{r}}\,|\,\psi_{n\boldsymbol{k}}\rangle\approx{}&\langle\tilde{u}_{n'\boldsymbol{k}-\boldsymbol{q}}\,|\,\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\boldsymbol{r}}\,|\,\tilde{u}_{n\boldsymbol{k}}\rangle+\sum_{ij,\mathrm{LM}}\langle\tilde{u}_{n'\boldsymbol{k}-\boldsymbol{q}}\,|\,\tilde{p}_{i\boldsymbol{k}-\boldsymbol{q}}\rangle\langle\tilde{p}_{j\boldsymbol{k}}\,|\,\tilde{u}_{n\boldsymbol{k}}\rangle\\ &\times\int\mathrm{e}^{-\mathrm{i}\boldsymbol{q}(\boldsymbol{r}-\boldsymbol{R}_i)}\hat{Q}_{ij}^{\mathrm{LM}}(\boldsymbol{r}-\boldsymbol{R}_i)\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r} \end{aligned} \tag{3.636} $$
式中：$\hat{Q}_{ij}^{\mathrm{LM}}$ 为单中心补偿电荷的多级展开。实际上，我们不仅需要计算交换电荷密度矩阵，为了通过方程（3.632）计算 GW 近似下的本征能级，还需要计算全电子波函数的交换能 $E_{\mathrm{xx}}$。在 PAW 框架下 $E_{\mathrm{xx}}$ 的计算式也比较复杂。关于 $\hat{Q}_{ij}^{\mathrm{LM}}$ 和 $E_{\mathrm{xx}}$ 的具体讨论已超出本书的范围，请参考文献\cite{kresse1999from,gajdos2006linear,paier2005perdew}。
2. 等离激元-极点近似
求解自能项 $\varSigma$ 最困难之处在于确定屏蔽库仑势 $W$。这是因为 $W$ 的计算要求我们知道微观介电函数矩阵的所有元素。当然，可以用方程（3.628）和方程（3.629）在每个 $\omega$ 下逐个求解矩阵元 $\epsilon_{\boldsymbol{q}}^{-1}(\boldsymbol{G},\boldsymbol{G}',\omega)$，但是更为常用的方法是利用所谓等离激元-极点近似（plasmon-pole approximation），将 $\epsilon_{\boldsymbol{q}}^{-1}(\boldsymbol{G},\boldsymbol{G}',\omega)$ 写为含参的解析表达式，再通过若干已知条件求解参数，从而得到整个微观介电函数的逆矩阵。该方法最初由 Hybertsen 和 Louie 实现\cite{hybertsen1986electron}。他们通过静态介电函数以及 $f$ 求和法则来确定参数。将 $\epsilon_{\boldsymbol{q}}^{-1}(\boldsymbol{G},\boldsymbol{G}',\omega)$ 的实部和虚部分别写为
$$ \mathrm{Re}\,\epsilon_{\boldsymbol{q}}^{-1}(\boldsymbol{G},\boldsymbol{G}',\omega)=\delta_{\boldsymbol{G}\boldsymbol{G}'}+\frac{\varOmega_{\boldsymbol{q}}^{2}(\boldsymbol{G},\boldsymbol{G}')}{\omega^{2}-\tilde{\omega}_{\boldsymbol{q}}^{2}(\boldsymbol{G},\boldsymbol{G}')} \tag{3.637} $$
$$ \operatorname{Im}\epsilon_{\boldsymbol q}^{-1}(\boldsymbol G,\boldsymbol G',\omega) =A_{\boldsymbol q}(\boldsymbol G,\boldsymbol G') \left[\delta(\omega-\tilde\omega_{\boldsymbol q}) -\delta(\omega+\tilde\omega_{\boldsymbol q})\right] \tag{3.638} $$
式中：$\varOmega_{\boldsymbol{q}}(\boldsymbol{G},\boldsymbol{G}')$ 称为等效裸等离激元频率（effective bare plasmon frequency），有
$$ \varOmega_{\boldsymbol{q}}^{2}(\boldsymbol{G},\boldsymbol{G}')=\omega_p^{2}\frac{(\boldsymbol{q}+\boldsymbol{G})\cdot(\boldsymbol{q}+\boldsymbol{G}')}{|\,\boldsymbol{q}+\boldsymbol{G}'\,|^{2}}\frac{\rho(\boldsymbol{G}-\boldsymbol{G}')}{\rho(0)} \tag{3.639} $$
其中
$$ {\omega_p^2=4\pi\rho(0)} $$
由 $\varOmega_{\boldsymbol{q}}(\boldsymbol{G},\boldsymbol{G}')$ 以及静态（$\omega=0$）介电函数矩阵，对于每一组 $\{\boldsymbol{q},\boldsymbol{G},\boldsymbol{G}'\}$，可得待定参数 $\tilde{\omega}$ 及 $A$：
$$ \tilde\omega_{\boldsymbol q}^2(\boldsymbol G,\boldsymbol G') =\frac{\Omega_{\boldsymbol q}^2(\boldsymbol G,\boldsymbol G')} {\delta_{\boldsymbol G\boldsymbol G'}-epsilon_{\boldsymbol q}^{-1}(\boldsymbol G,\boldsymbol G',0)} \tag{3.640} $$
$$ A_{\boldsymbol q}(\boldsymbol G,\boldsymbol G') =-\frac{\pi\Omega_{\boldsymbol q}^2(\boldsymbol G,\boldsymbol G')} {2\tilde\omega_{\boldsymbol q}(\boldsymbol G,\boldsymbol G')} \tag{3.641} $$
当前比较流行的软件包，如 Abinit 以及 GPAW 等则用了两个频率下 $\epsilon^{-1}$ 的结果求解 $\tilde{\omega}$ 和 $A$。这里给出最后结果：
$$ \tilde\omega_{\boldsymbol q}^2(\boldsymbol G,\boldsymbol G') =E_0^2\frac{\delta_{\boldsymbol G\boldsymbol G'}-epsilon_{\boldsymbol q}^{-1}(\boldsymbol G,\boldsymbol G',\mathrm iE_0)} {\epsilon_{\boldsymbol q}^{-1}(\boldsymbol G,\boldsymbol G',\mathrm iE_0) -\epsilon_{\boldsymbol q}^{-1}(\boldsymbol G,\boldsymbol G',0)} \tag{3.642} $$
$$ A_{\boldsymbol q}(\boldsymbol G,\boldsymbol G') =-\frac{\pi\tilde\omega_{\boldsymbol q}(\boldsymbol G,\boldsymbol G')}{2} \left[\delta_{\boldsymbol G\boldsymbol G'}-epsilon_{\boldsymbol q}^{-1}(\boldsymbol G,\boldsymbol G',0)\right] \tag{3.643} $$
其中所选频率分别为 0 和 $\mathrm{i}E_0$。$E_0$ 的取值需要仔细选择，一般取 1 Hartree。
Bethe-Salpeter 方程
Bethe-Salpeter 方程（Bethe-Salpeter equation，BSE）是描述两体相互作用系统的一种方法。它源自量子场论，是量子电动力学（QED）中的一个重要工具，用于描述粒子间的相互作用。在凝聚态物理和材料科学中，BSE 被用来求体系的激发态性质，例如光吸收谱，包括直接和间接带隙材料的激发态。
BSE 通过考虑电子-空穴相互作用，去除了 DFT 和 GW 近似中的一些限制。在这些方法中，电子被视为独立粒子，在外加势场下运动。然而，电子并不是独立的，而是相互作用的。这种相互作用在许多情况下对物理性质，例如激发态性质、超导性和磁性等有着决定性的影响。因此，为了得到更准确的结果，我们需要一种能够考虑电子相互作用的方法。BSE 就是这样一种方法。
BSE 建立在 GW 近似的基础之上。在 GW 近似中，电子的自能被写成了一个有效势的形式，这个有效势是所有电子的库仑相互作用的平均。这意味着 GW 近似中的电子被看成是在其他电子产生的平均场中运动的。这是一个合理的近似， 其频率依赖的屏蔽相互作用 $W(\omega)$ 已描述动态屏蔽。BSE 在准粒子描述上加入电子–空穴相互作用，可得到更准确的激发态。BSE 可以看作描述激子的一个方程，其中激子是电子-空穴对。电子从价带被激发到导带，留下了一个空穴，电子和空穴之间由于库仑相互作用而形成了一个激子。 BSE 在选定的准粒子能量与相互作用核近似下，计算激子能量及其波函数，并据此预测光吸收谱；结果精度取决于这些近似和数值收敛。
BSE 的数值解可以通过迭代的方法得到。一般的步骤是，首先使用 DFT 得到基态波函数和能量，然后用 GW 近似修正单粒子能级，最后通过解 BSE 得到激子能级。这个过程通常需要大量的计算资源，但是随着算法和计算资源的发展，BSE 已经在越来越多的系统上得到了应用。
BSE 的具体形式为
$$ L(1234)=G(13)G(24)+\int\mathrm{d}(5678)\,G(15)G(25)K(5678)L(7834) \tag{3.644} $$
式中：核 $K$ 包含准粒子相互作用，分为电子-空穴交换作用 $v$ 与电子-空穴吸引作用 $W$，具体可表示为
$$ K(5678)=v(57)\delta(56)\delta(78)-W(56)\delta(57)\delta(68) \tag{3.645} $$
图 3.20 为 BSE 相应的费曼图。将 $L(1234)$ 变换到频率空间，则有
$$ L(1234,\omega)=\frac{1}{H^{2p}-\omega} \tag{3.646} $$
式中：$H^{2p}$ 是有效两体哈密顿量，且有
$$ H_{n_1n_2}^{n_3n_4,2p}=(\varepsilon_{n_2}-\varepsilon_{n_1})\delta_{n_1n_3}\delta_{n_2n_4}+(f_{n_1}-f_{n_2})(V_{n_1n_2}^{n_3n_4}-W_{n_1n_2}^{n_3n_4}(\omega)) \tag{3.647} $$
其中
$$ V_{n_1n_2}^{n_3n_4}=\iint\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}'\,\psi_{n_1}(\boldsymbol{r})\psi_{n_2}^{*}(\boldsymbol{r})\frac{1}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\psi_{n_3}^{*}(\boldsymbol{r}')\psi_{n_4}(\boldsymbol{r}') \tag{3.648} $$
$$ W_{n_1n_2}^{n_3n_4}(\omega)=\iint\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}'\,\psi_{n_1}(\boldsymbol{r})\psi_{n_2}^{*}(\boldsymbol{r})\frac{\epsilon^{-1}(\boldsymbol{r},\boldsymbol{r}',\omega)}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\psi_{n_3}^{*}(\boldsymbol{r}')\psi_{n_4}(\boldsymbol{r}') \tag{3.649} $$
图 3.20 BSE 相应的费曼图
做简化时，一般将方程（3.649）中的介电函数用静态介电函数 $\epsilon^{-1}(\boldsymbol{r},\boldsymbol{r}',\omega=0)$ 近似。在平面波基组中，可以将 $H^{2p}$ 的各个 $\boldsymbol{q}$ 分量表示如下（忽略上标 $2p$）\cite{onida2002electronic,benedict1998optical,rohlfing1998electron}：
$$ \begin{aligned} H^{\mathrm{BSE}}_{vc\boldsymbol k,v'c'\boldsymbol k'}(\boldsymbol Q) ={}&\bigl(\varepsilon_{c,\boldsymbol k+\boldsymbol Q}^{\mathrm{QP}} -\varepsilon_{v,\boldsymbol k}^{\mathrm{QP}}\bigr) \delta_{vv'}\delta_{cc'}\delta_{\boldsymbol k\boldsymbol k'}\\ &+K^{x}_{vc\boldsymbol k,v'c'\boldsymbol k'}(\boldsymbol Q) +K^{d}_{vc\boldsymbol k,v'c'\boldsymbol k'}(\boldsymbol Q). \end{aligned} \tag{3.650} $$
$$ \begin{aligned} K^x_{vc\boldsymbol k,v'c'\boldsymbol k'}(\boldsymbol Q) =\iint\!\mathrm d\boldsymbol r\,\mathrm d\boldsymbol r'\; &\psi^*_{c,\boldsymbol k+\boldsymbol Q}(\boldsymbol r) \psi_{v,\boldsymbol k}(\boldsymbol r) \frac{1}{|\boldsymbol r-\boldsymbol r'|}\\[-1mm] &\times\psi^*_{v',\boldsymbol k'}(\boldsymbol r') \psi_{c',\boldsymbol k'+\boldsymbol Q}(\boldsymbol r'). \end{aligned} \tag{3.651} $$
$$ \begin{aligned} K^d_{vc\boldsymbol k,v'c'\boldsymbol k'}(\boldsymbol Q) =-\iint\!\mathrm d\boldsymbol r\,\mathrm d\boldsymbol r'\; &\psi^*_{c,\boldsymbol k+\boldsymbol Q}(\boldsymbol r) \psi_{c',\boldsymbol k'+\boldsymbol Q}(\boldsymbol r)\\[-1mm] &\times W(\boldsymbol r,\boldsymbol r';0) \psi^*_{v',\boldsymbol k'}(\boldsymbol r') \psi_{v,\boldsymbol k}(\boldsymbol r'). \end{aligned} \tag{3.652} $$
为了得到 BSE，必须利用 GW 近似所得到的本征能级以及由 Kohn-Sham 方程所得到的本征波函数。因此，Rubio 等人将 Kohn-Sham 方法、GW 方法及 BSE 方法三种方法形象地称为计算体系激发谱的三级助推火箭。与 3.7.3 节的讨论类似，可以通过谱表示构造 $L(1234,\omega)$，其中谱函数 $A$ 由下面的本征值方程给出：
$$ \sum_{v'c'\boldsymbol k'}H^{\mathrm{BSE}}_{vc\boldsymbol k,v'c'\boldsymbol k'}(\boldsymbol Q) A^S_{v'c'\boldsymbol k'}(\boldsymbol Q) =\Omega_S(\boldsymbol Q)A^S_{vc\boldsymbol k}(\boldsymbol Q) \tag{3.653} $$
需要注意，完整 BSE 的共振与反共振耦合矩阵通常不是普通厄米矩阵，但这并不意味着其激发能必为复数。对于稳定体系的静态 BSE，物理激发能可以是实数；激子寿命须由另外的衰减机制或相应的复自能计算。为了简化计算，只考虑 $H^{2p}$ 的共振部分，即由被占据的价带 $v$ 向非占据的导带 $c$ 的跃迁。由此可以得到 BSE 近似下的宏观介电函数：
$$ \begin{aligned} \epsilon_M(\omega)&=\lim_{\boldsymbol q\to0} \frac{1}{\epsilon^{-1}_{\boldsymbol q}(\boldsymbol G=0,\boldsymbol G'=0;\omega)},\\ \operatorname{Im}\epsilon_M(\omega)&=\frac{8\pi^2e^2}{\omega^2} \sum_S\left|\sum_{vc\boldsymbol k}A^S_{vc\boldsymbol k} \hat{\boldsymbol e}\!\cdot\!\langle v\boldsymbol k|\boldsymbol v|c\boldsymbol k\rangle\right|^2 \delta(\omega-\Omega_S),\quad\boldsymbol Q\simeq0. \end{aligned} \tag{3.654} $$
$\epsilon_M(\omega)$ 的虚部给出了体系的光吸收谱。更为详细的讨论可参阅文献\cite{onida2002electronic,gajdos2006linear}。