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 |
对每行的说明:
- 行 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: |
往下一点是相对 PAW 原子的总能量成分分析。从 GPAW 分析来看,Si 在形成晶体后,电子动能增加,势能与交换关联能降低。后面两者容易理解,但原子形成固体后动能增加这一点从化学成键相悖。只能说在赝势 PAW 下,直接的能量组成并没有物理意义. VASP 的 OUTCAR 中的能量组成分析没有给出动能的部分.
Energy contributions relative to reference atoms: (reference = -15772.688500) |
再往下是关于能带结构的简单信息。包括费米能级和前两个 k 点上的价带导带本征值和占据数。需要注意的是这里的 Occupancy 是该自旋轨道上的电子数乘以 k 点分数权重的值。例如 $\Gamma$ 点权重是固定的 1/512, 价带全占满有 2 个电子,2 / 512 = 0.00391。如果是自旋极化计算,这部分会给出两个自旋通道的占据数,占据数会变成非极化的一半. $\Gamma$ 点带隙为 2.56 eV.
Fermi level: 5.73099 |
能带计算
PBE 能带计算需要使用 SCF 得到的电子密度。官方网站给出的一个能带计算例子如下
calc = GPAW ('Si_gs.gpw', nbands=16, fixdensity=True, symmetry='off', |
这里主要需要理解 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 能带计算结果 |
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 ([[ |
获取各 k 点的直接带隙
vb = int (calc.get_number_of_electrons ()/2.0) - 1 |
手动制作能带图的方法可以参考 Band Structures With GPAW - Mantid Project.