GW近似 - PhyRevue

普通 Kohn–Sham DFT ：

\left[ -\frac{1}{2}\nabla^2 +V_{\mathrm{ext}} +V_H +V_{\mathrm{xc}} \right] \psi_{n\mathbf{k}} = \varepsilon^{\mathrm{KS}}_{n\mathbf{k}} \psi_{n\mathbf{k}}

其中：

$V_{\mathrm{ext}}$ ：原子核对电子的外势；
$V_H$ ：电子密度产生的经典静电库仑势；
$V_{\mathrm{xc}}$ ：交换-关联势，负责处理 Hartree 势中没有包含的量子多体效应。

KS 方程主要帮助我们得到准确的基态电子密度，但其本征值 $\varepsilon^{\mathrm{KS}}_{n\mathbf{k}}$ 并不能保证是真实的电子激发能。

例如，PBE 使用的交换-关联势 $V_{\mathrm{xc}}$ 主要根据某个位置附近的电子密度及其梯度来近似描述复杂的电子相互作用。但是，当“往材料中加入一个电子”后，体系中其他电子会发生重新分布，其产生的多体响应具有明显的非局域性和动态性，而 $V_{\mathrm{xc}}$ 对这些效应的描述有限。

GW 近似：

假设往材料中加入一个电子，它会排斥周围的电子，使周围电子重新分布，并形成一个极化云。

因此，实际传播的不再是一个完全独立的裸电子，而可以理解为：

\text{电子} + \text{周围电子产生的极化云} \Longrightarrow \text{quasiparticle（准粒子）}

由 Dyson 方程推导得到的准粒子方程：描述准粒子（quasiparticle）的能量和波函数

\left[ -\frac12\nabla^2 +V_{\mathrm{ext}} +V_H \right]\psi^{QP}_n(\mathbf r) + \int \Sigma(\mathbf r,\mathbf r',E_n^{QP}) \psi^{QP}_n(\mathbf r') d\mathbf r' = E_n^{QP}\psi^{QP}_n(\mathbf r)

这里的 $\Sigma=\text{self-energy}$ ，即自能。

对比准粒子方程和 Kohn–Sham DFT 方程，可以将自能 $\Sigma$ 理解为对交换-关联势 $V_{\mathrm{xc}}$ 的升级：

V_{\mathrm{xc}} \longrightarrow \Sigma(E)

其中 $V_{\mathrm{xc}}$ 通常是局域或半局域、静态的势，而 $\Sigma$ 一般是非局域、能量依赖的算符。

GW 近似最核心的公式为：

\Sigma(\mathbf r,\mathbf r',\omega) = \frac{i}{2\pi} \int d\omega' G(\mathbf r,\mathbf r',\omega+\omega') W(\mathbf r,\mathbf r',\omega')

即：

\Sigma=iGW

其中：

$G$ ：单粒子格林函数；
$W$ ：屏蔽库仑相互作用。

在 $G_0W_0$ 中，通常从 DFT 得到的波函数和本征值构造非相互作用格林函数 $G_0$ ：

G_0(\mathbf r,\mathbf r',\omega) = \sum_{n\mathbf k} \frac{ \psi_{n\mathbf k}(\mathbf r) \psi^*_{n\mathbf k}(\mathbf r') }{ \omega-\varepsilon_{n\mathbf k} \pm i\eta }

屏蔽库仑相互作用 $W$ 为：

W_{GG'}(\mathbf q,\omega) = \epsilon^{-1}v = \epsilon^{-1}_{GG'}(\mathbf q,\omega) \frac{4\pi}{|\mathbf q+\mathbf G'|^2}

其中：

v(\mathbf r,\mathbf r') = \frac{1}{|\mathbf r-\mathbf r'|}

是裸 Coulomb 相互作用；

\epsilon = \epsilon_{\mathbf G\mathbf G'}(\mathbf q,\omega)

是介电函数；

$\mathbf q$ 是转移动量，对应：

\mathbf k \rightarrow \mathbf k+\mathbf q

而：

\mathbf G,\mathbf G' = \text{倒格矢}

。

$\epsilon^{-1}$ 的来源：

体系受到总电势 $V_{\mathrm{tot}}$ 后，电子密度会发生变化。在线性响应范围内：

\delta n = P V_{\mathrm{tot}}

其中 $P$ 是极化函数（polarizability），描述电子密度对外界电势扰动的响应。

电子密度发生变化后，又会通过 Coulomb 相互作用产生感应电势：

V_{\mathrm{ind}} = v\,\delta n = vP V_{\mathrm{tot}}

其中 $v$ 是裸 Coulomb 相互作用。

总电势由外部电势和感应电势组成：

V_{\mathrm{tot}} = V_{\mathrm{ext}} + V_{\mathrm{ind}} = V_{\mathrm{ext}} + vP V_{\mathrm{tot}}

整理得：

(1-vP)V_{\mathrm{tot}} = V_{\mathrm{ext}}

定义介电函数：

\epsilon = 1-vP

于是：

\epsilon V_{\mathrm{tot}} = V_{\mathrm{ext}}

因此：

V_{\mathrm{tot}} = \epsilon^{-1}V_{\mathrm{ext}}

所以，$\epsilon^{-1}$ 描述了材料中的电子屏蔽如何把原始外部电势 $V_{\mathrm{ext}}$ 转化为体系内部实际感受到的总电势 $V_{\mathrm{tot}}$ 。

也可以从电子的多次响应来理解 $W$ ：

W = v+vPv+vPvPv+\cdots

1.第一次电子响应：

V_{\mathrm{ext}} = v\rho_{\mathrm{ext}}

\rho_{\mathrm{ext}} \xrightarrow{v} V_{\mathrm{ext}} \xrightarrow{P} \delta n^{(1)}

2.新产生的电荷密度产生新的库仑势：

V = (v+vPv)\rho_{\mathrm{ext}}

\delta n^{(1)} \xrightarrow{v} V_{\mathrm{ind}}^{(1)}

3.新的库仑势再次作用于电子，引起第二次响应：

V = (v+vPv+vPvPv)\rho_{\mathrm{ext}}

V_{\mathrm{ind}}^{(1)} \xrightarrow{P} \delta n^{(2)} \xrightarrow{v} V_{\mathrm{ind}}^{(2)}

............

一直继续下去：

V_{\mathrm{tot}} = (v+vPv+vPvPv+\cdots)\rho_{\mathrm{ext}}

定义 $W$ ：

V_{\mathrm{tot}} = W\rho_{\mathrm{ext}}

因此得到：

\begin{aligned} W &= v+vPv+vPvPv+\cdots \\ &= (1-vP)^{-1}v \\ &= \epsilon^{-1}v \end{aligned}

极化函数P

在 $G_0W_0$ 中：

P_0 = -iG_0G_0

而：

G_0(\mathbf r,\mathbf r',\omega) = \sum_{n\mathbf k} \frac{ \psi_{n\mathbf k}(\mathbf r) \psi^*_{n\mathbf k}(\mathbf r') }{ \omega-\varepsilon_{n\mathbf k} \pm i\eta }

将两个 $G_0$ 代入，可以得到独立粒子极化函数：

\Longrightarrow P_0(\mathbf r,\mathbf r',\omega) = 2 \sum_{n\mathbf k} \sum_{m\mathbf k'} \frac{ f_{n\mathbf k} - f_{m\mathbf k'} }{ \omega+ \varepsilon_{n\mathbf k} - \varepsilon_{m\mathbf k'} +i\eta } \, \psi^*_{n\mathbf k}(\mathbf r) \psi_{m\mathbf k'}(\mathbf r) \psi^*_{m\mathbf k'}(\mathbf r') \psi_{n\mathbf k}(\mathbf r')

做 Fourier transform，将其投影到平面波基组：

P^0_{GG'}(\mathbf q,\omega) = \frac{1}{\Omega} \int d\mathbf r \int d\mathbf r' \, e^{-i(\mathbf q+\mathbf G)\cdot\mathbf r} P_0(\mathbf r,\mathbf r',\omega) e^{+i(\mathbf q+\mathbf G')\cdot\mathbf r'}

得到：

\begin{aligned} P^0_{GG'}(\mathbf q,\omega) = \frac{2}{\Omega} \sum_{n\mathbf k} \sum_{m\mathbf k'} & \frac{ f_{n\mathbf k}-f_{m\mathbf k'} }{ \omega+\varepsilon_{n\mathbf k} -\varepsilon_{m\mathbf k'} +i\eta } \\ &\times \int d\mathbf r\, e^{-i(\mathbf q+\mathbf G)\cdot\mathbf r} \psi^*_{n\mathbf k}(\mathbf r) \psi_{m\mathbf k'}(\mathbf r) \\ &\times \int d\mathbf r'\, e^{+i(\mathbf q+\mathbf G')\cdot\mathbf r'} \psi^*_{m\mathbf k'}(\mathbf r') \psi_{n\mathbf k}(\mathbf r'). \end{aligned}

由于晶体具有平移对称性，动量守恒要求：

\mathbf k' = \mathbf k+\mathbf q

因此定义跃迁矩阵元：

M^G_{nm}(\mathbf k,\mathbf q) = \int d\mathbf r\, \psi^*_{n\mathbf k}(\mathbf r) e^{-i(\mathbf q+\mathbf G)\cdot\mathbf r} \psi_{m,\mathbf k+\mathbf q}(\mathbf r)

即：

M^G_{nm}(\mathbf k,\mathbf q) = \left\langle n\mathbf k \left| e^{-i(\mathbf q+\mathbf G)\cdot\mathbf r} \right| m,\mathbf k+\mathbf q \right\rangle

最后得到：

\begin{aligned} P^0_{GG'}(\mathbf q,\omega) &= \frac{2}{\Omega} \sum_{\mathbf k} \sum_{n,m} \frac{ f_{n\mathbf k} - f_{m,\mathbf k+\mathbf q} }{ \omega + \varepsilon_{n\mathbf k} - \varepsilon_{m,\mathbf k+\mathbf q} +i\eta } M^G_{nm}(\mathbf k,\mathbf q) M^{G'*}_{nm}(\mathbf k,\mathbf q) \\ &= \frac{2}{\Omega} \sum_{\mathbf k} \sum_v^{\rm occ} \sum_c^{\rm empty} M^G_{vc} (M^{G'}_{vc})^* \left[ \frac{1} {\omega-\Delta\varepsilon_{vc}+i\eta} - \frac{1} {\omega+\Delta\varepsilon_{vc}-i\eta} \right] \end{aligned}

其中：

\Delta\varepsilon_{vc} = \varepsilon_{c,\mathbf k+\mathbf q} - \varepsilon_{v,\mathbf k}

$\Omega$ 是晶体体积，$f_{n\mathbf k}$ 是占据数。

这里可以直观地理解为：

\boxed{ P_0 = \text{体系中所有可能的占据态} \rightarrow \text{空态跃迁的总响应} }

接下来：

P_0 \rightarrow \epsilon \rightarrow \epsilon^{-1} \rightarrow W

1.从 $P_0$ 构造介电矩阵 $\epsilon$ 。

公式：

\epsilon = 1-vP_0

在平面波基组中：

\epsilon_{GG'}(\mathbf q,\omega) = \delta_{GG'} - \sum_{G''} v_{GG''}(\mathbf q) P^0_{G''G'}(\mathbf q,\omega)

由于裸 Coulomb 相互作用在平面波基组中是对角的：

v_{GG'}(\mathbf q) = \delta_{GG'} \frac{4\pi}{|\mathbf q+\mathbf G|^2}

因此：

\epsilon_{GG'}(\mathbf q,\omega) = \delta_{GG'} - \frac{4\pi} {|\mathbf q+\mathbf G|^2} P^0_{GG'}(\mathbf q,\omega)

对介电矩阵求逆得到：

\epsilon^{-1}_{GG'}(\mathbf q,\omega)

然后构造屏蔽库仑相互作用 $W$ ：

\begin{aligned} W_{GG'}(\mathbf q,\omega) &= \sum_{G''} \epsilon^{-1}_{GG''}(\mathbf q,\omega) v_{G''G'}(\mathbf q)\\ &= \epsilon^{-1}_{GG'}(\mathbf q,\omega) v_{G'}(\mathbf q)\\ &= \epsilon^{-1}_{GG'}(\mathbf q,\omega) \frac{4\pi} {|\mathbf q+\mathbf G'|^2} \end{aligned}

因此整个过程为：

G_0 \Longrightarrow P_0=-iG_0G_0 \Longrightarrow \epsilon=1-vP_0 \Longrightarrow W_0=\epsilon^{-1}v

再带入：

\Sigma=iG_0W_0

得到自能 $\Sigma$ 。

最后求准粒子能量：

E_{n\mathbf k}^{QP} = \varepsilon_{n\mathbf k}^{KS} + Z_{n\mathbf k} \left\langle \psi_{n\mathbf k} \left| \Sigma(E_{n\mathbf k}^{QP}) - V_{xc} \right| \psi_{n\mathbf k} \right\rangle

其中 $Z_{n\mathbf k}$ 是准粒子重整化因子，用来考虑自能 $\Sigma(E)$ 对能量的依赖：

Z_{n\mathbf k} = \left[ 1- \left. \frac{\partial \operatorname{Re}\Sigma_{n\mathbf k}(E)} {\partial E} \right|_{E\approx\varepsilon^{KS}_{n\mathbf k}} \right]^{-1}

因此可以把整个 $G_0W_0$ 的计算流程概括为：

\boxed{ \text{DFT} \rightarrow \{\psi^{KS},\varepsilon^{KS}\} \rightarrow G_0 \rightarrow P_0 \rightarrow \epsilon \rightarrow W_0 \rightarrow \Sigma=iG_0W_0 \rightarrow E^{QP} }