本教程旨在帮助初学者快速掌握二维材料热电性质计算方法。欢迎广大学者阅览，如有错误请及时批评指正。教程包括：

结构优化

声子谱和声子态密度计算

从头算分子动力学（AIMD）计算

材料的能带结构（band structure）和态密度（DOS）计算

考虑自旋轨道相互作用（SOC）和范德瓦尔斯力（vdW）下材料的能带结构计算

弛豫时间的计算

热输运性质的计算

电输运参数和ZT值的计算

注：蓝色代表已经更新，请看上一期：

天玑算-科研服务-独家 | （一）二维材料热电性质计算方法

天玑算-科研服务-独家 | (二) 二维材料热电性质计算方法—AIMD，能带，DOS

天玑算科技，为科研助力。今天小天就为大家更新考虑自旋轨道相互作用（SOC）和范德瓦尔斯力（vdW）下材料的能带结构计算和弛豫时间的计算方法。

后续还会继续更新，敬请期待。制作烧脑，码字不易，还请大家多多点赞在看，转发分享哦。

计算教程

五、计算SOC和vdW下材料的能带结构

为了考虑SOC和vdW的影响，我们还要计算SOC和vdW下材料的能带结构。

（一）计算SOC下材料的能带结构

在正常计算能带的基础上，从结构优化开始，在每个INCAR中加入如下几行：

LSORBIT=.T.                                        

SAXIS= 001            

GGA_COMPAT=.F.                

ISPIN=2

注意：Vasp执行文件不能用vasp_stp了，要改用vasp_ncl(Vasp保证是5.4.1版本以上，否则需要重新编译Vasp)

（二）计算vdW下材料的能带结构

在正常计算能带的基础上，从结构优化开始，在每个INCAR中加入一行：

IVDW=12

六、基于DP（Deformation Potential）理论求出弛豫时间

根据DP（deformation potential）理论，可以根据如下公式求出弛豫时间：

其中，图片（普朗克常量）和 图片（玻尔兹曼常数）都是常数，所以也可以写成如下形式：

A是常数，C2D是弹性常数，T是绝对温度，图片是有效质量，E1是DP常数。这些参数的计算方法如下：

（一）计算C2D

根据定义C2D = [图片2E/图片(Δa/a0)2]/S0，我们对优化好的结构在x或者y方向施加应力使它产生-0.03到0.03的应变即用优化好的结构的POSCAR的x或y方向的晶格常数。

如图8，分别×(0.97-1.03)，间隔0.01，这样就产生了14个POSCAR（x方向7个，y方向7个，其中有两个是重复的即不施加应力的两个×1的POSCAR，所以应该有13个POSCAR文件，放在13个文件夹里），

然后分别在每个的INCAR中设置ISIF=2，将得到的结构再分别优化一遍，得到每个稳定结构的能量（OSZICAR中最后一行读取E0），

然后用origin作图，横坐标为-0.03-0.03，间隔0.01，纵坐标为读取的体系的能量，这样在固定方向上就有7行数据然后作图拟合求二次元方程，得到的二次幂乘以2即为所求二阶偏导图片2E/图片(Δa/a0)2，

然后再除以面积（未施加应力下x方向的晶格常数乘以y方向的晶格常数）即得到C2D

图8

（二）计算E1

根据公式E1 = 图片Eedge/图片(Δa/a0)，需要计算每个应变下导带底和价带顶的移动，所以首先要计算这13个结构的能带，然后在得到的本征值文件（EIGENVAL）里面找到能带x和y方向价带顶和导带底对应的点，

方法为先找到其所在能带上的位置，然后找到其所在的路径进而找到对应于EIGENVAL的点，比如如果导带价带都在Γ点，那么对应的点就为(0,0,0)，在EIGENVAL里找到(0,0,0)这个点(如图2所示)，

然后在下面列出的数据中找到有电子和无电子的两个点即为VBM和CBM（当最右边一列从1变成0对应的两行）即所对应的能量（左边的一列，单位是eV），将每个文件的VBM和CBM的能量导入origin中在每个能级下减去真空能级作为纵坐标。横坐标和第一步的一样，为-0.03-0.03，间隔0.01，作出直线，其斜率即为E1。

图9

真空能级的详细求法可以参考https://www.bigbrosci.com/2018/09/03/ex49/

本教程简单整理如下：

在结构优化完毕的基础上，将结构优化中的CONTCAR,POTCAR,KPOINTS拷贝到新文件夹下，CONTCAR重命名为POSCAR，INCAR设置如下：

SYSTEM=AlISe

ENCUT=500

ISTART=0

ICHARG=2

ISMEAR=0

SIGMA=0.05

PREC=Accurate

EDIFF=1.0E-4

IBRION=2

NSW=1

LVHAR=.TRUE.

LWAVE=.FALSE.

LCHARG=.FALSE.

ALGO=Normal

LPLANE=.TRUE.

LREAL=.FALSE.

生成的LOCPOT导入到P4VASP中查看结果：

a)打开 p4vasp 并导入 LOCPOT 文件：file –> loadsystem，然后找到 LOCPOT 所在的目录，点解导入，效果如下

b)开始可视化 Electronic –> Local potential

效果如下：

点击下图中的 Show，会出现三个方向让你选择，一般都是沿着 z 方向的，按照红色箭头点下即可。

效果如下：

c)导出数据

同能带计算，点击：Graph–>Export，得到xx.dat文件

d)分析数据

用 Notepadqq打开xx.dat 文件。

打开之后，我们会发现，数据共有三部分，每两部分之间由空行分开。

这三部分的顺序，如下图，

我们从图中可以看出，基本上在10-18横坐标（第一列）的这个区间范围，三条线的数据差不多是一样的。静电势能的话，可以通过红线(Average那条线)，也就是第二部分来获取。可以查看下14-15时红色曲线的数值，如下图：

可以看出：此材料的真空能级大概为4.84 eV

（三）计算m*

根据公式m*=图片2/(图片2E/图片k2)，我们需要求出未产生应变下能带中x和y方向价带顶和导带底附近曲线的二阶导，

具体操作为：在origin中找到沿x和y方向价带顶和导带底附近的能量数据（一般从CBM或者VBM开始取三四个点就可以了），所得能量除以27.21即为纵坐标，横坐标为倒空间内对应点与Γ点的距离，以Al2I2Se2 monolayer 为例，倒空间内每个点与Γ点的距离求法如下：

Al2I2Se2 monolayer原胞的晶格常数为：

换算成以a.u为单位（1a.u = 0.5291 Å）：

Al2I2Se2 monolayer 能带所对应的K点路径为G(0,0,0) → X(0.5,0,0) → S(0.5,0.5,0) → Y(0,0.5,0)

所以X和G(Γ)点在倒空间内的距离GX=|0.5图片|=0.31875，因为我们设置K-path的时候设置的G → X有20个点，如图10所示，所以G→ X每两个点的距离为图片，同理，GY==|0.5图片|=0.521425，G → Y每两个点的距离为图片=0.0274434。

图10

每个能量对应的距离得到后，就可以用origin作图了，得到半个抛物线，将这半个抛物线做个对称得到一条完整的抛物线，比如对于沿x方向的在Γ点的导带底，我们对沿x方向的导带底的曲线做关于x=0的对称，得到一条完整的抛物线，

然后求这条抛物线的二阶导，将二阶导求倒数即为所求某方向电子或者空穴的有效质量。

迁移率μ满足：

求出弛豫时间后，即可求出迁移率。（可通过把其他论文中的数据带入公式求出常数项A(2421.786)进而求出弛豫时间τ）图片取1.682。