GPAW 笔记 (二)——DFT 自洽场与能带计算

以金刚石结构 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._calc 的 energy 属性，操作上是从 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 为一个字典，包含 path 和 npoints 两个关键字. 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

▲