GPAW笔记(四)——GW计算

本文将简单介绍GPAW中的进行准粒子计算的G0W0类, 并使用GWBands类制作单层二硫化钼的能带结构图.

前言

无论是否自洽, GW计算都需要一组单粒子态作为输入, 这组单粒子态中需要包含大量的非占据态. 目前主流是采用Kohn-Sham DFT产生的Kohn-Sham本征态. 这里主要介绍非自洽计算情形, 此时GW结果受到输入影响.

为了得到包含较多非占据态的本征态一般有两种方法. 第一种是在基态计算中包含大量的非占据态, 这种情况电子步迭代效率低. 第二种是首先在较少非占据态下得到收敛的电子密度(电子密度只与占据态有关), 然后在固定电子密度下对角化Kohn-Sham单电子Hamiltonian. 对于局域和半局域泛函, 当电子密度固定时, 久期行列式就完全确定了, 一步即可得到等于基组数量的非占据态. 对于非局域泛函, 也只需要几步以收敛非局域势算符.

这里用的例子是单层二硫化钼. 首先用ASE构造二硫化钼模型. 晶格常数相比实验上体相稳定结构略有拉伸.

from ase.build import mx2
mos2 = mx2(formula='MoS2', kind='2H', a=3.19, thickness=3.127, size=(1, 1, 1), vacuum=5.0)

基态计算与对角化

关于基态和对角化计算已经在GPAW 笔记 (二)——DFT 自洽场与能带计算和GPAW笔记(三)——求解器对全哈密顿量对角化的影响中提及. 这里用了一个比较小的平面波截断以减小计算量. 而k值设的比较大. 若k太小, 后面用GWBands作能带插值时VBM和CBM位置会出错.

from gpaw import GPAW, PW, FermiDirac

Ecut = 300
kx = 18
calc = GPAW(mode=PW(Ecut), xc='PBE',
            convergence={"density": 1e-8},
            kpts={'size': (kx, kx, 1), 'gamma': True},
            occupations=FermiDirac(0.001), parallel={'domain': 1},
            txt='gs.txt')
mos2.set_calculator(calc)
# 基态SCF
mos2.get_potential_energy()
calc.write('gs.gpw')
# 对角化哈密顿量, 将波函数和能量写入文件, 以供GW读取
calc.diagonalize_full_hamiltonian()
calc.write('fulldiag.gpw', 'all')

GW计算

这里用包含4s4p的Mo PAW setup, 因此体系共有14+6+6=26个价电子, 在自旋非极化下有13个占据态, 因此价带顶的指标是12. 现在计算其前5个价带和导带.

from gpaw.response.g0w0 import G0W0

diagfile = 'fulldiag.gpw'
ecut = 50

gw = G0W0(calc=diagfile, bands=(8, 18), # VB at index 12
          method="G0W0", ecut=50, nblocksmax=True,
          truncation='2D', q0_correction=True,
          domega0=0.03, omega2=10,
          filename='g0w0', savepckl=True)
gw.calculate()

其中nblocksmax设为True时GPAW将最大化响应函数chi0的分块, 减小每个进程的内存消耗. 8进程测试的内存占用和GW总用时结果如下表, 开启nblocksmax后预测的单进程内存消耗减小, 但是计算消耗时间更长.

nlocksmax Estimate Mem. per proc. (MB) wall time (s)
True 30 1923
False 233 1429

truncation和q0_correction是针对低维体系的参数, 前者加快对于真空层厚度的收敛, 后者则是加快对面内布里渊区采样格点的收敛. 具体可以参考文献1和3.

参数domega0和omega2与频率积分有关. domega0决定第一个频率点的位置, omega2决定在何处倍增格点间距. 一般来说带隙越小, domega0就要取得越小以对低频响应采样充分. 显然的, domega0越小, omega2越大, 频率格点数越多, 计算量越大. 目前这个设置可以使K点带隙收敛到2 meV以内.

作能带图

import pickle
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
from gpaw.response.gw_bands import GWBands

calcfile = 'fulldiag.gpw'
gwpckl = 'g0w0_results.pckl'

# 设置K点路径
K = np.array([1 / 3, 1 / 3, 0])
M = np.array([0.0, 0.5, 0])
G = np.array([0.0, 0.0, 0.0])
kpoints = np.array([G, K/2, K, M, G])

# 初始化GWBands对象. bandrange要和G0W0设置一致, 不过这里第二个指标是要包含进去的
GW = GWBands(calcfile=calcfile, gwpckl=gwpckl, kpoints=kpoints, bandrange=(8,17))
# 设置
gwbopts={"nk_Int": 100, "interpolate": True, "vac": False}
# 提取PBE, GW和GW-SOC能带插值数据
pbe = GW.get_gw_bands(SO=False, dft=True, **gwbopts)
gw = GW.get_gw_bands(SO=False, **gwbopts)
gwsoc = GW.get_gw_bands(SO=True, **gwbopts)
# 提取K点路径的一维坐标
x_x = gw['x_k']
X = gw['X']/x_x[-1]
x_x /= x_x[-1]
# 对齐VBM
ePBE_kn, eGW_kn, eGWsoc_kn = [d['e_kn'] - d['vbm'] for d in [pbe, gw, gwsoc]]
# 作不同方法得到的能带
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
styles = [
          {'ls': '-.', 'marker': '', 'color': 'k', "lw": 3},
          {'ls': '--', 'marker': '', 'color': '#393b79', "lw": 3},
          {'ls': '-', 'marker': '', 'color': '#d62728', "lw": 3},
         ]
for i, bands in enumerate([ePBE_kn, eGW_kn, eGWsoc_kn]):
    ax.plot(x_x, bands, **styles[i])
# 能量零点
ax.axhline(0.0, color='k', linestyle=':', lw=2)
# 图例
leg_handles = [mpl.lines.Line2D([], [], **style) for style in styles]
leg_labels = [r'PBE', r'G$_0$W$_0$', r'G$_0$W$_0$-SOC']
ax.legend(leg_handles, leg_labels, fontsize=20)

ax.set_xlim(0, x_x[-1])
ax.set_ylim([-2, 4])
ax.set_ylabel('Energy (eV)', fontsize=24)
# 横坐标特殊k点, 并作分割线
for p in X[:-1]:
    plt.axvline(p, color='#AAAAAA', ls='--', lw=2)
labels_K = [r'$\Gamma$', r'$T$', r'$K$', r'$M$', r'$\Gamma$']
plt.xticks(X, labels_K, fontsize=18)

plt.yticks(fontsize=17)
fig.tight_layout()
plt.savefig('MoS2_band_GTKMG.png', dpi=300)
plt.show()

效果如下图. 可以很容易看到, 在不包含SOC情况下, PBE预测单层MoS2具有$\Gamma-K$的间接带隙, 而GW给出的是K上的直接带隙. 包含SOC会导致K点能带裂分, 在GW下进一步增大$\Gamma$点和K点VBM的能差.

PBE, GW方法得到的MoS2能带
PBE, GW方法得到的MoS2能带

另外值得注意的一点是CB在$\Gamma-K$上也有一个能量较低的态(在T=0.5K附近). 在DFT下$T_c$与$K_c$能量差肉眼可见, 而在GW尤其是包含SOC下, 这两个态基本是简并的. 考虑到对未占据态的自能修正为正值, 这说明$K_c$的自能修正要大于$T_c$. 定性分析上, $T_c$主要是Mo的$d_{x^2-y^2}$和$d_{xy}$在面内成键, 而$K_c$是Mo的$d_{z^2}$.(文献2) 后者较为定域, 从GW修正DFT离域误差的角度可定性理解.

参考资料

  1. Ismail-Beigi, S. Phys. Rev. B 73, 233103 (2006)
  2. Zhang, L.; Zunger, A. Nano Lett. 15, 949-957 (2015)
  3. Rasmussen, F. et al. Phys. Rev. B 94, 155406 (2016)
  4. Quasi-particle spectrum of two-dimensional materials - GPAW tutorial

Comments

▲