以金刚石结构 Si 为例,演示了用 GPAW 计算材料电子能带结构的过程,对结果进行了简单分析.

背景

一年前我在 GPAW笔记(一)——安装及测试 一文中介绍了 GPAW 的安装和简单测试。当时主要是为了比较 GPAW, abinit 和 VASP 的 GW 效率,于是稍微研究了一下 GPAW 中参数含义,运行了一些官方网站上的脚本。在相近的参数下作了三个程序的 GW 计算交差后,就没有再管,当时也没有留下清楚的学习笔记.

由于最近实际计算的需要,我又重新学习 GPAW, 着重于 DFT 能带计算,GW 和 BSE. 这一篇笔记主要记录了基础 DFT 计算部分的学习,包括 SCF 和能带,简单分析标准输出。这里使用的结构是金刚石型的硅,晶格常数 $a=5.43$ (A).

SCF 计算

从官网例子中改编得的一个 SCF 计算脚本如下.

from ase.build import bulk
from gpaw import GPAW, PW, FermiDirac

si = bulk ('Si', 'diamond', a=5.43)
calc = GPAW (mode=PW (200), xc='PBE', kpts={'size': (8,8,8), 'gamma': True},
random=True, occupations=FermiDirac (0.01), txt='Si_gs.txt')
si.calc = calc
si.get_potential_energy ()
calc.write ('Si_gs.gpw')

对每行的说明:

  • 行 4: 用 bulk 方法构造 Si 模型
  • 行 5-6: 初始化 GPAW. 使用平面波基组 (截断为 200 eV). 泛函为 PBE. 布里渊区采样为 $\Gamma$ 中心 $8\times8\times8$ 的均匀格点。初始波函数随机生成。占据数计算采用 Fermi-Dirac 分布,. txt 选项指定输出文件.
  • 行 7: 将 Si 模型的 calculator 指定为刚刚初始化好的 GPAW.
  • 行 8: 开始 SCF 过程。具体来说,Atoms.get_potential_energy 方法用于获取 Atoms._calcenergy 属性,操作上是从 Calculator.results 字典中提取关键字 energy 值。对于未收敛的 GPAW calculator, energy 关键字不存在,此时会执行 GPAW.calculate 进行 SCF 循环,直到收敛,添加 energy 并返回之.
  • 行 9: 将计算过程的所有信息写入 Si_gs.gpw 文件中.

标准输出分析

标准输出保存在 Si_gs.txt 中。主要看三部分。在 200 行附近给出 SCF 迭代的流程

                     log10-error:    total        iterations:
time wfs density energy fermi poisson
iter: 1 14:49:29 +1.03 +inf -10.866303 7
iter: 2 14:49:30 -1.97 -0.80 -10.866989 2
iter: 3 14:49:30 -1.01 -0.81 -11.320728 22
iter: 4 14:49:31 -1.94 -0.56 -10.847586 3
iter: 5 14:49:31 -2.71 -1.11 -10.806417 8
iter: 6 14:49:31 -2.69 -1.35 -10.780669 9
iter: 7 14:49:32 -3.97 -2.59 -10.780801 0
iter: 8 14:49:32 -5.11 -2.71 -10.780798 0
iter: 9 14:49:33 -5.83 -2.71 -10.780763 0
iter: 10 14:49:33 -6.76 -3.19 -10.780765 0
iter: 11 14:49:34 -8.02 -3.59 -10.780765 0
iter: 12 14:49:34 -6.82 -3.60 -10.780765 0
iter: 13 14:49:34 -8.06 -3.60 -10.780764 0
iter: 14 14:49:35 -10.08 -3.78 -10.780764 0
iter: 15 14:49:35 -8.45 -3.78 -10.780764 0
iter: 16 14:49:36 -8.00 -3.68 -10.780764 0
iter: 17 14:49:36 -8.17 -3.80 -10.780764 0
iter: 18 14:49:37 -9.30 -3.92 -10.780764 0
iter: 19 14:49:37 -8.99 -3.95 -10.780764 0
iter: 20 14:49:37 -9.83 -4.06 -10.780764 0

Converged after 20 iterations.

往下一点是相对 PAW 原子的总能量成分分析。从 GPAW 分析来看,Si 在形成晶体后,电子动能增加,势能与交换关联能降低。后面两者容易理解,但原子形成固体后动能增加这一点从化学成键相悖。只能说在赝势 PAW 下,直接的能量组成并没有物理意义. VASP 的 OUTCAR 中的能量组成分析没有给出动能的部分.

Energy contributions relative to reference atoms: (reference = -15772.688500)

Kinetic: +15.782785
Potential: -13.855376
External: +0.000000
XC: -12.678254
Entropy (-ST): -0.000000
Local: -0.029919
--------------------------
Free energy: -10.780764
Extrapolated: -10.780764

再往下是关于能带结构的简单信息。包括费米能级和前两个 k 点上的价带导带本征值和占据数。需要注意的是这里的 Occupancy 是该自旋轨道上的电子数乘以 k 点分数权重的值。例如 $\Gamma$ 点权重是固定的 1/512, 价带全占满有 2 个电子,2 / 512 = 0.00391。如果是自旋极化计算,这部分会给出两个自旋通道的占据数,占据数会变成非极化的一半. $\Gamma$ 点带隙为 2.56 eV.

Fermi level: 5.73099

Showing only first 2 kpts
Kpt Band Eigenvalues Occupancy
0 2 5.31933 0.00391
0 3 5.31934 0.00391
0 4 7.87829 0.00000
0 5 7.87829 0.00000

1 2 4.51004 0.02344
1 3 4.51006 0.02344
1 4 7.32815 0.00000
1 5 9.02227 0.00000

能带计算

PBE 能带计算需要使用 SCF 得到的电子密度。官方网站给出的一个能带计算例子如下

calc = GPAW ('Si_gs.gpw', nbands=16, fixdensity=True, symmetry='off',
kpts={'path': 'GXWKL', 'npoints': 60}, convergence={'bands': 8})
calc.get_potential_energy ()
calc.write ('Si_bs.gpw')

这里主要需要理解 GPAW 一行

  • 第一个参数 Si_gs.gpw 为读入文件.
  • nbands 为 SCF 迭代所包含的能带数.
  • fixdensity=True, 顾名思义,固定电子密度.
  • kpts 为一个字典,包含 pathnpoints 两个关键字. path 包含 BZ 特殊点记号。这里的路径为 G-X-W-K-L, 根据总 k 点数 npoints 自动设置所需要计算的 k 点,不会重复计算特殊 k 点。特殊点之间的 k 点数不是均匀的.
  • convergence 为一个字典,包含关键字 band, 值为 8. 表示使用最低的 8 个能级的波函数的 Kohn-Sham 方程余矢量模方作为收敛判据。可以使用 ‘all’ 和 -8 之类的负值,后者表示收敛除最后 8 个外的所有能级。具体参考 Accuracy of the self-consistency cycle 条目.

get_potential_energy 开始 SCF 循环,结束后将 calculator 写入 Si_bs.gpw.

GPAW 提供了能带分析的帮助方法 band_structure. 该方法继承自 ASE 的 Calculator 类。使用方式如下

calc = GPAW ('Si_bs.gpw')    # 读取 Si 能带计算结果
bs = calc.band_structure () # 调用 band_structure 方法
bs.plot (show=True, emax=10.0, filename="Si_bs.png") # 作图

band_structure 方法返回一个 ase.dft.band_structure.BandStructure 类。得到能带图如下,和 官网例子 是相同的.

与电子结构相关的 GPAW 方法

BandStructure 作图非常方便,但如果是自己作图或者需要能级数据做进一步处理,就需要从 GPAW 对象直接获取 k 点和 Kohn-Sham 本征值。一些可能用到的方法如下

方法 继承 作用 参数 返回
get_fermi_level gpaw.paw.PAW 获取 Fermi 能级 float
get_number_of_spins gpaw.paw.PAW 获取自旋通道数量 int
get_number_of_electrons gpaw.paw.PAW 获取价电子总数 float
get_eigenvalues gpaw.paw.PAW 获取本征值 (eV) (kpt=0, spin=0) 1D array
get_occupation_numbers gpaw.paw.PAW 获取轨道占据数 (kpt=0, spin=0) 1D array
get_ibz_k_points gpaw.paw.PAW 获取不可约 k 点 2D array
get_bz_k_points gpaw.paw.PAW 获取所有 k 点 2D array

例如,获取所有能级到 eigens 数组

eigens = np.array ([[
calc.get_eigenvalues (ik, isp) for ik,_ in enumerate (calc.get_ibz_k_points ())
] for isp in range (calc.get_number_of_spins ())])

获取各 k 点的直接带隙

vb = int (calc.get_number_of_electrons ()/2.0) - 1
direct_gaps = eigens [0, :, vb+1] - eigens [0, :, vb]

手动制作能带图的方法可以参考 Band Structures With GPAW - Mantid Project.

Comments