vasp 系列 002. 通过 vaspkit 生成的 BNAD.dat 数据文本计算有效质量

1. 背景

1. 现在使用 vasp 软件计算，基本都是用 vaspkit 数据后处理工具。vaspkit 工具有计算有效质量的功能，但是需要重新计算能带。
2. 在计算每一个材料的时候，绝大多数情况都要计算能带，而能带的数据可以使用 vaspkit 来提取，得到一个 BAND.dat 的能带数据文件。

BAND.dat 长这样：

#K-Path      Energy-Level
# NKPTS & NBANDS: 240  96
# Band-Index    1
   0.00000    -12.501256
   0.01166    -12.500967
   .
   .
   .
   2.51673    -12.501256
 
 
# Band-Index    2
   2.51673    -12.500515
   2.50327    -12.500128
   .
   .
   .
   2.51673    -12.501256
....

2. 环境

安装 Anaconda （一个集成了 Python 绝大多数库的工具）

---

---

---

3. 有效质量相关理论

3.1 要明白：

只有 有带隙 的材料才能用这个脚本计算有效质量，最好是 绝缘体 和 半导体 ，金属的价带和导带是有重叠的是没法用下面的脚本计算有效质量。

3.2 倒格矢

3.2.1. 正格矢

$${\mathbf{R}}_{\mathbf{n}}=n_1{\mathbf{a}}_1+n2{\mathbf{a}}_2+n3{\mathbf{a}}_3$$

下面例子的单位是 埃（1E-10 米）

---

网上最流行的二维 InSe 教程的 POSCAR

SYSTEM = InSe                           
   1.00000000000000     
     4.0927800000000003    0.0000000000000000    0.0000000000000000
    -2.0463909999999998    3.5444536547999999    0.0000000000000000
     0.0000000000000000    0.0000000000000000   20.0000000000000000
   Se   In
     2     2
Direct
 -0.0000000000000000  0.0000000000000000  0.6346282943308172
  0.0000000000000000 -0.0000000000000000  0.3653717056691828
  0.6666666870000029  0.3333333429999996  0.5704161898487420
  0.6666666870000029  0.3333333429999996  0.4295838101512581

$${\mathbf{a}}_1=(4.09278,0,0)$$
$${\mathbf{a}}_2=(-2.0463,3.5444536,0)$$
$${\mathbf{a}}_3=(0,0,20)$$
$$|{\mathbf{a}}_2|=\sqrt{(-2.0463{)}^2+(3.5444536{)}^2}=4.0927824$$

单位是 埃，编辑器打不出来。

---

3.2.2. 倒格矢（单位：1/正格矢单位，比如下面的例子：1/angstrom）

三维情况下：
$${\mathbf{G}}_{\mathbf{h}}=h_1{\mathbf{b}}_1+h_2{\mathbf{b}}_2+h_3{\mathbf{b}}_3$$

---

其中：
$${\mathbf{b}}_1=2\pi \frac{({\mathbf{a}}_2\times {\mathbf{a}}_3)}{{\mathbf{a}}_1\cdot ({\mathbf{a}}_2\times {\mathbf{a}}_3)},$$
$${\mathbf{b}}_2=2\pi \frac{({\mathbf{a}}_3\times {\mathbf{a}}_1)}{{\mathbf{a}}_1\cdot ({\mathbf{a}}_2\times {\mathbf{a}}_3)},$$
$${\mathbf{b}}_3=2\pi \frac{({\mathbf{a}}_1\times {\mathbf{a}}_2)}{{\mathbf{a}}_1\cdot ({\mathbf{a}}_2\times {\mathbf{a}}_3)},$$

---

或者（有的书是没有乘以 2 π 的，比如 原版 的黄昆的书）：
$${\mathbf{b}}_1=\frac{({\mathbf{a}}_2\times {\mathbf{a}}_3)}{{\mathbf{a}}_1\cdot ({\mathbf{a}}_2\times {\mathbf{a}}_3)},$$
$${\mathbf{b}}_2=\frac{({\mathbf{a}}_3\times {\mathbf{a}}_1)}{{\mathbf{a}}_1\cdot ({\mathbf{a}}_2\times {\mathbf{a}}_3)},$$
$${\mathbf{b}}_3=\frac{({\mathbf{a}}_1\times {\mathbf{a}}_2)}{{\mathbf{a}}_1\cdot ({\mathbf{a}}_2\times {\mathbf{a}}_3)},$$

---

---

---

下面以上面的 InSe 的 POSCAR 为例子，按照有 × 2π 的情况计算倒格矢

1. 计算分母，也就是体积：
   $$V={\mathbf{a}}_1\cdot ({\mathbf{a}}_2\times {\mathbf{a}}_3)$$

In [10]: import numpy as np
In [11]: a1 = np.array([4.09278, 0, 0])
    ...: a2 = np.array([-2.04639, 3.54445365, 0])
    ...: a3 = np.array([0, 0, 20])
    ...: 
    ...: # 求分母，也就是体积
    ...: V = np.dot(a1, np.cross(a2, a3))
    ...: print(V)
290.13338019294

2. 计算各个倒格矢

In [16]: a1 = np.array([4.09278, 0, 0])
    ...: a2 = np.array([-2.04639, 3.54445365, 0])
    ...: a3 = np.array([0, 0, 20])
    ...: 
    ...: # 求分母，也就是体积
    ...: V = np.dot(a1, np.cross(a2, a3))
    ...: 
    ...: b1 = 2 * np.pi * np.cross(a2, a3) / V
    ...: b2 = 2 * np.pi * np.cross(a3, a1) / V
    ...: # 二维情况不考虑 b3
    ...: # b3 = 2 * np.pi * np.cross(a1, a2) / V
    ...: 
    ...: # 打印倒格矢向量形式
    ...: print(b1)
    ...: print(b2)
    ...: # print(b3)
    ...: 
    ...: # 打印倒格矢向量的模（长度）
    ...: print(np.linalg.norm(b1))
    ...: print(np.linalg.norm(b2))
    ...: # print(np.linalg.norm(b3)
[ 1.53518765  0.88634045 -0.        ]
[0.         1.77268091 0.        ]
1.772681730317445
1.772680905893518

可以看到这个例子的倒格矢$${\mathbf{b}}_1,{\mathbf{b}}_2$$ 的长度都是 1.77268 (1/angstrom)。

3.3 有效质量的定义（单位下面具体讲）

$$\frac{1}{m^{\ast }}=\frac{1}{{\hslash }^2}\frac{{\partial }^2{E}_k}{\partial {\mathbf{k}}^2}$$

解释：

在能带的极值处对 倒格矢坐标 求二阶导，然后代进公式计算，便可求得有效质量

其中：

如果计算的是电子有效质量，则取得能带是最低的导带，在极小值处附近取点做二阶拟合，或者求二阶导。

如果计算的是空穴有效质量，则取得能带是最高的价带，在极大值处附近取点做二阶拟合，或者求二阶导。

比如在 InSe 的例子中：（红圈的是求电子有效质量，最低的导带极小值附近；绿框的是求空穴有效质量，最高的价带极大值附近），那么计算的时候就是通过找到这几处地方求其二阶拟合（或二阶导数）。

单位（根据能带图）：
$$\hslash :1.0545\times 10^{-34}J\cdot s$$
$$E_k:eV$$
$$\mathbf{k}:1/angstrom$$

---

---

3.4 有效质量的单位：

$$\frac{1}{m^{\ast }}:\frac{eVangstro{m}^2}{(J^2\cdot s^2)}$$

现在问题是为什么 E_k 对 k 的二阶导的单位不是分子分母都是平方，而是只有分母平方，而分子没有平方。
在这里举个加速度的例子：
$$a=\frac{d^{2}x}{d{t}^2}$$
和加速度的单位一样的道理：
$$a的加速度不为m^2{s}^{-2}，而是m{s}^{-2}$$

所以， 有效质量的单位是：
$$m^{\ast }:J^2\cdot s^2(eVangstro{m}^2{)}^{-1}$$
这个单位比较难看，并且化简不方便。但是也可以化简：
$$J:kg\cdot m^2{s}^{-2}$$
$$eV:1.6\times 10^{-19}J$$
$$m^{\ast }:J^2\cdot s^2(eVangstro{m}^2{)}^{-1}=\frac{J\cdot (kg\cdot m^2{s}^{-2})s^2}{1.6\times 10^{-19}J\cdot (10^{-10}m{)}^2}=(\ast \ast )kg$$
最后得到的是质量单位 kg，但数值很复杂。

---

所以，在计算的时候，我们都是用 Hartree atomic units：

https://en.wikipedia.org/wiki/Hartree_atomic_units

$$ReducedPlanckconstant:\hslash =1hartree=1.0545\times 10^{-34}J\cdot s$$
$$Elementarycharge:e=1hartree=1.6\times 10^{-19}C$$
$$Bohrradius:a_0=1hartree=5.29\times 10^{-11}m$$
$$Electronmass:m_{\text{e}}=1hartree=9.1\times 10^{-31}kg$$
$$Energy:E_h=1hartree=27.21eV$$

---

化简后（Hartree atomic units）：
$$\frac{1}{m^{\ast }}=\frac{{\partial }^2{E}_k}{\partial {\mathbf{k}}^2}$$

这时候，E_k 的单位是 Hartree，k 的单位也是 Hartree，h_bar 是 1 省略不写。所以最后得出的有效质
量单位是 Hartree，即单位为元电荷（这也是论文经常看见的）： $$1Hartree=1m_e(m_0)=9.1\times 10^{-31}kg$$

---

---

3.5 有效质量的求法（基于 Hartree atomic units）:

$$\frac{1}{m^{\ast }}=\frac{{\partial }^2{E}_k}{\partial {\mathbf{k}}^2}$$
单位是：
$$1Hartree=1m_e(m_0)=9.1\times 10^{-31}kg$$

1. E_k 对 k 求二阶导，观察极值附近的数据点，可以近似是一条抛物线，可以对其进行二次拟合（为了对应 Origin ，二次系数用 c 表示）：
   $$f(x)=c{x}^2+bx+a\Rightarrow f^{{}^{\prime \prime }}(x)=2c$$
   即：
   $$\frac{1}{m^{\ast }}=\frac{{\partial }^2{E}_k}{\partial {\mathbf{k}}^2}=2c$$
   $$m^{\ast }=\frac{1}{2c}，单位：m_e(m_0)$$
   这也是网上求有效质量的方法：

https://mp.weixin.qq.com/s?__biz=MzI2OTQ4OTExOA==&mid=2247484762&idx=1&sn=fff244bc64af036f52ccecb1bec28eef&chksm=eadec3ebdda94afd846972e5ae363107eacec316c8db7c3b328305451ea66d5870bcc597a516&token=1986072339&lang=zh_CN#rd

或者

https://zhuanlan.zhihu.com/p/62573921 点击转 “有效质量详细教程”

2. 也可以使用插值法直接求解数据的二阶导。

---

---

---

4. BAND.dat 数据文件格式

数据格式（每一条带 240 个 K 点，一共 96 条带，能量单位是 eV，x 轴即 K-path 的单位是 1/A）：

#K-Path      Energy-Level
# NKPTS & NBANDS: 240  96
# Band-Index    1
   0.00000    -12.501256
   0.01166    -12.500967
   .
   .
   .
   2.51673    -12.501256
 
 
# Band-Index    2
   2.51673    -12.500515
   2.50327    -12.500128
   .
   .
   .
   2.51673    -12.501256
....

这里要明白一点：

x 轴（K-PATH）的数值是如何来的！

先看 KPOINTS

K-Path Generated by VASPKIT.
   80
Line-Mode
Reciprocal
   0.0000000000   0.0000000000   0.0000000000     G       
   0.5000000000   0.0000000000   0.0000000000     M              
 
   0.5000000000   0.0000000000   0.0000000000     M              
   0.3333333333   0.3333333333   0.0000000000     K              
 
   0.3333333333   0.3333333333   0.0000000000     K              
   0.0000000000   0.0000000000   0.0000000000     G        

三段高对称点的路径，每一段路径 80 个 K 点，一共 240 个 K 点，所以看 BAND.dat 第 80 个 K 点的 K-PATH 的数值就是高对称点 M (0.5, 0, 0)。

从 3.2.2 的计算得出 b1 = 1.772681730317445 1/A
b1 / 2 = 0.8863408651587225 1/A

显然和 高对称点 M (0.5, 0, 0) 的值一样，所以可以判断出，vaspkit 导出的 K-path 的数据是采用 3.2.2 中有 × 2 π 的形式

5. Python 脚本及其使用见：

使用： Python 脚本通过 vaspkit 生成的能带数据文本 BAND.dat 计算有效质量

脚本链接： 脚本网页下载