摘要

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

背景

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

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

SCF计算

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

1
2
3
4
5
6
7
8
9
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迭代的流程

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
                     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中的能量组成分析没有给出动能的部分.

1
2
3
4
5
6
7
8
9
10
11
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$点带隙为2.56 eV.

1
2
3
4
5
6
7
8
9
10
11
12
13
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得到的电子密度. 官方网站给出的一个能带计算例子如下.

1
2
3
4
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类. 使用方式如下

1
2
3
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 获取本征值 (ik=0, ispin=0) 1D array
get_ibz_k_points gpaw.paw.PAW 获取不可约k点 2D array
get_bz_k_points gpaw.paw.PAW 获取所有k点 2D array

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

1
2
3
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点的直接带隙

1
2
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