GW Approximation - PhyRevue

Ordinary 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}}

where:

$V_{\mathrm{ext}}$ : the external potential from the atomic nuclei acting on the electrons;
$V_H$ : the classical electrostatic Coulomb potential generated by the electron density;
$V_{\mathrm{xc}}$ : the exchange–correlation potential, which accounts for the quantum many-body effects not included in the Hartree potential.

The KS equations mainly help us obtain an accurate ground-state electron density, but their eigenvalues $\varepsilon^{\mathrm{KS}}_{n\mathbf{k}}$ are not guaranteed to be true electron excitation energies.

For example, the exchange–correlation potential $V_{\mathrm{xc}}$ used in PBE approximates the complex electron interactions mainly from the electron density and its gradient near a given position. However, after adding an electron to the material, the other electrons in the system redistribute, and the resulting many-body response is clearly nonlocal and dynamic; $V_{\mathrm{xc}}$ has only a limited ability to describe these effects.

The GW approximation:

Suppose we add an electron to a material. It repels the surrounding electrons, causing them to redistribute and forming a polarization cloud.

Therefore, what actually propagates is no longer a completely independent bare electron; it can be understood as:

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

The quasiparticle equation derived from the Dyson equation: describes the energy and wavefunction of a 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)

Here $\Sigma=\text{self-energy}$ , i.e., the self-energy.

Comparing the quasiparticle equation with the Kohn–Sham DFT equation, the self-energy $\Sigma$ can be understood as an upgrade of the exchange–correlation potential $V_{\mathrm{xc}}$ :

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

Here $V_{\mathrm{xc}}$ is usually a local or semilocal, static potential, whereas $\Sigma$ is generally a nonlocal, energy-dependent operator.

The central formula of the GW approximation is:

\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')

that is:

\Sigma=iGW

where:

$G$ : the single-particle Green's function;
$W$ : the screened Coulomb interaction.

In $G_0W_0$ , the noninteracting Green's function $G_0$ is usually constructed from the wavefunctions and eigenvalues obtained from DFT:

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 }

The screened Coulomb interaction $W$ is:

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

where:

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

is the bare Coulomb interaction;

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

is the dielectric function;

$\mathbf q$ is the transferred momentum, corresponding to:

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

and:

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

Origin of $\epsilon^{-1}$ :

After the system is subjected to the total potential $V_{\mathrm{tot}}$ , the electron density changes. Within the linear response regime:

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

where $P$ is the polarizability, which describes the response of the electron density to an external potential perturbation.

After the electron density changes, it in turn produces an induced potential through the Coulomb interaction:

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

where $v$ is the bare Coulomb interaction.

The total potential consists of the external potential and the induced potential:

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

Rearranging gives:

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

Define the dielectric function:

\epsilon = 1-vP

Thus:

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

Therefore:

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

So $\epsilon^{-1}$ describes how electronic screening in the material converts the original external potential $V_{\mathrm{ext}}$ into the total potential $V_{\mathrm{tot}}$ actually felt inside the system.

$W$ can also be understood through the multiple responses of the electrons:

W = v+vPv+vPvPv+\cdots

First electron response:

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

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

The newly generated charge density produces a new Coulomb potential:

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

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

The new Coulomb potential acts on the electrons again, causing a second response:

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

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

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

Continuing indefinitely:

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

Define $W$ :

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

Therefore:

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

Polarization function P

In $G_0W_0$ :

P_0 = -iG_0G_0

and:

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 }

Substituting the two $G_0$ 's gives the independent-particle polarizability:

\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 transforming it and projecting it onto a plane-wave basis:

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'}

we obtain:

\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}

Because the crystal has translational symmetry, momentum conservation requires:

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

Therefore define the transition matrix element:

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)

that is:

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

Finally we get:

\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}

where:

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

$\Omega$ is the crystal volume, and $f_{n\mathbf k}$ is the occupation number.

Here it can be intuitively understood as:

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

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

Construct the dielectric matrix $\epsilon$ from $P_0$ .

Formula:

\epsilon = 1-vP_0

In a plane-wave basis:

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

Since the bare Coulomb interaction is diagonal in a plane-wave basis:

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

Therefore:

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

Inverting the dielectric matrix gives:

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

Then construct the screened Coulomb interaction $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}

Thus the whole process is:

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

Substituting this back into:

\Sigma=iG_0W_0

gives the self-energy $\Sigma$ .

Finally, the quasiparticle energy is obtained from:

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

where $Z_{n\mathbf k}$ is the quasiparticle renormalization factor, which accounts for the energy dependence of the self-energy $\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}

Therefore the entire $G_0W_0$ computational workflow can be summarized as:

\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} }