专栏名称: 研之成理

夯实基础，让基础成就辉煌；传递思想，让思想改变世界。“研之成理科研平台”立足于科研基础知识与科研思想的传递与交流，旨在创建属于大家的科研乐园！主要内容包括文献赏析，资料分享，科研总结，论文写作，软件使用等。科研路漫漫，我们会一路陪伴你！

理论计算干货：VASP计算二维材料的载流子迁移率

研之成理 · 公众号 · 科研 · 2019-11-01 10:29

正文

1、前言

载流子迁移率 通常指半导体内部电子和空穴整体的运动快慢情况，是衡量半导体器件性能的重要物理量。2004年，石墨烯的成功剥离引起了研究人员对于二维材料性质探索的浓厚兴趣。石墨烯、黑磷等二维材料展现出的高载流子迁移率是其中的一个重要研究课题，科研人员在理论计算方面已经做了大量的工作。由于电子在运动过程中不仅受到外电场力的作用，还会不断的与晶格、杂质、缺陷等发生无规则的碰撞，大大增加了理论计算的难度。

目前 计算载流子迁移率 比较常用的理论是 形变势理论 和 玻尔兹曼输运理论 ，前者没有考虑电子和声子(晶格振动)以及电子与电子之间的相互作用等因素，计算结果存在一定的误差，但笔者的计算结果与实验值在数量级上是吻合的；玻尔兹曼输运理论的一种计算考虑了电子-声子的相互作用，基于第一性原理计算和最大局域化 Wannier 函数插值方法，借助于 Quantum-ESPRESSO 和 EPW 软件可以完成载流子迁移率计算。缺点是计算量太大，一般的课题组很难承受起高昂的计算费用，另外 EPW 软件对于二维材料的计算存在部分问题，在其官方论坛也有讨论，计算过程在后续文章中会提到。

本文以形变势理论方法为基础，详细介绍了二维 InSe 的电子和空穴的有效质量与载流子迁移率的计算方法。

2、理论基础

基于 Bardeen 和 Shockley 提出的形变势理论，二维材料载流子迁移率可以根据下式计算：

其中， m ^∗ 是传输方向上的有效质量， T 是温度， k _B 是玻尔兹曼常数。 E ₁ 表示沿着传输方向上位于价带顶 (VBM) 的空穴或聚于导带底 (CBM) 的电子的形变势常数，由公式 E ₁ =Δ E /(Δ l / l ₀ ) 确定，Δ E 为在压缩或拉伸应变下 CBM 或 VBM 的能量变化， l ₀ 是传输方向上的晶格常数，Δ l 是 l ₀ 的变形量。 m _d 是载流子的平均有效质量，由下面公式定义。

C _2D 是均匀变形晶体的弹性模量，对于 2D 材料，弹性模量可以通过下面公式来计算，

其中 E 是总能量， S ₀ 是优化后的面积。

下面对公式中的单位(量纲)做一个简单换算，具体如下：

换算过程：

3、计算与数据处理工具

VASP.5.4.4软件可以手动控制优化晶格方向
OriginLab软件
Excel
Materials Studio软件
正格矢到倒格矢转化脚本，来源于小木虫(见下面链接)

http://muchong.com/bbs/viewthread.php?tid=7149817&fpage=1

#! /usr/bin/python# This program reads in base vectors from a given file, calculates reciprocal 
vectors# then writes to outfile in different units# LinuxUsage: crecip.py infile outfile# Note:  the infile must be in the form below:#   inunit  ang/bohr#   _begin_vectors#     46.300000000   0.000000000   0.000000000#      0.000000000  40.500000000   0.000000000#      0.000000000   0.000000000   10.000000000#   _end_vectors# # Note: LATTICE VECTORS ARE SPECIFIED IN ROWS !def GetInUnit( incontent ):  inunit = ""  for line in incontent:      if line.find("inunit") == 0:          inunit = line.split()[1]          break




    
  return inunitdef GetVectors( incontent ):  indstart = 0  indend = 0  for s in incontent:      if s.find("_begin_vectors") != -1:          indstart = incontent.index(s)      else:          if s.find("_end_vectors") != -1:              indend = incontent.index(s)  result = []  for i in range( indstart + 1, indend ):      line = incontent[i].split()      result.append( [ float(line[0]), float(line[1]), float(line[2]) ] )  return resultdef Ang2Bohr( LattVecAng ):  LattVecBohr = LattVecAng  for i in range(0,3):      for j in range(0,3):          LattVecBohr[i][j] = LattVecAng[i][j] *  1.8897261246  return LattVecBohrdef DotProduct( v1, v2 ):  dotproduct = 0.0  for i in range(0,3):      dotproduct = dotproduct + v1[i] * v2[i]  return dotproductdef CrossProduct( v1, v2 ):  # v3 = v1 WILL LEAD TO WRONG RESULT  v3 = []  v3.append( v1[1] * v2[2] - v1[2] * v2[1] )  v3.append( v1[2] * v2[0] - v1[0] * v2[2] )  v3.append( v1[0] * v2[1] - v1[1] * v2[0] )  return v3def CalcRecVectors( lattvec ):  pi = 3.141592653589793  a1 = lattvec[0]  a2 = lattvec[1]  a3 = lattvec[2]  b1 = CrossProduct( a2, a3 )  b2 = CrossProduct( a3, a1 )  b3 = CrossProduct( a1, a2 )  volume = DotProduct( a1, CrossProduct( a2, a3 ) )  RecVec = [ b1, b2, b3 ]  # it follows the definition for b_j: a_i * b_j = 2pi * delta(i,j)  for i in range(0,3):      for j in range(0,3):          RecVec[i][j] = RecVec[i][j] * 2 * pi / volume  return RecVec    def main(argv = None):    argv = sys.argv  infilename  = argv[1]  outfilename = argv[2]      pi = 3.141592653589793  bohr2ang = 0.5291772109253  ang2bohr = 1.889726124546      infile  = open(infilename,"r")  incontent = infile.readlines()  infile.close()      inunit = GetInUnit( incontent )  LattVectors = GetVectors( incontent )  # convert units from ang to bohr  if inunit == "ang":      LattVectors = Ang2Bohr( LattVectors )    




    
  # calculate reciprocal vectors in 1/bohr  RecVectors = CalcRecVectors( LattVectors )      # open outfile for output  ofile = open(outfilename,"w")      # output lattice vectors in bohr  ofile.write("lattice vectors in bohr:\n")  for vi in LattVectors:      ofile.write("%14.9f%14.9f%14.9f\n" % (vi[0], vi[1], vi[2]))  ofile.write("\n")      # output lattice vectors in ang  convfac = bohr2ang  ofile.write("lattice vectors in ang:\n")  for vi in LattVectors:      ofile.write("%14.9f%14.9f%14.9f\n" % (vi[0]*convfac, vi[1]*convfac,
 




vi[2]*







convfac))  ofile.write("\n")      # output reciprocal vectors in 1/bohr  ofile.write("reciprocal vectors in 1/bohr:\n")  for vi in RecVectors:      ofile.write("%14.9f%14.9f%14.9f\n" % (vi[0], vi[1], vi[2]))  ofile.write("\n")      # output reciprocal vectors in 1/ang  convfac = ang2bohr  ofile.write("reciprocal vectors in 1/ang:\n")  for vi in RecVectors:      ofile.write("%14.9f%14.9f%14.9f\n" % (vi[0]*convfac, vi[1]*convfac, vi[2]*
convfac))  ofile.write("\n")      # output reciprocal vectors in 2pi/bohr  convfac = 1.0/(2.0*pi)  ofile.write("reciprocal vectors in 2pi/bohr:\n")  for vi in RecVectors:      ofile.write("%14.9f%14.9f%14.9f\n" % (vi[0]*convfac, vi[1]*convfac, vi[2]*convfac))  ofile.write("\n")  # output reciprocal vectors in 2pi/ang  convfac = ang2bohr/(2.0*pi)  ofile.write("reciprocal vectors in 2pi/ang:\n")  for vi in RecVectors:      ofile.write("%14.9f%14.9f%14.9f\n" % (vi[0]*convfac, vi[1]*convfac, vi[2]*convfac))      # close  ofile.close()      return 0if __name__ == "__main__": import sys sys.exit(main())

4、二维 InSe 有效质量计算过程

4.1 建模

由于计算过程中需要对二维 InSe 施加应变，但二维 InSe 原胞是六角结构，不容易施加应变。但是侯柱峰老师讲了对石墨烯原胞施加应变的方法，笔者认为虽然可行，但过于繁琐，故不采用此法。我们可以利用根号建模的方法讲六角结构 InSe原胞变为方形结构的 InSe 超胞，然后施加应变可大大提高操作效率，但计算量的增加再可接受范围之内。下面给出关键的建模步骤，更多的根号建模部分可参考我的往期博客文章。

切面并构建二维InSe原胞，同时调整晶格基矢，使其变为方形结构：

4.2 结构优化

INCAR

SYSTEM = InSe             ISTART = 0       NWRITE = 2      PREC   = AccurateENCUT  = 500GGA    = PE      NSW    = 200ISIF   = 3         ISYM   = 2        IBRION = 2      NELM   = 80        EDIFF  = 1E-05         EDIFFG = -0.01 ALGO   = Normal            LDIAG  = .TRUE.  LREAL  = .FALSE.          ISMEAR = 0         SIGMA  = 0.05 ICHARG = 2LWAVE  =  .FALSE.          LCHARG = .FALSE.NPAR   = 4

KPOINTS

Monkhorst Pack0Gamma11 7 1.0 .0 .0

POSCAR

Se In1.000 4.083622259999999 -0.000000000000001 0.000000000000000 0.000000000000000 7.073041233239241 0.000000000000000 0.000000000000000 0.000000000000000 25.377516849029359Se In4 4Direct0.5000005000000000 0.1666665000000000 0.5271404971815050 !Se0.0000004999999997 0.6666665000000004 0.5271404971815050 !Se0.5000005000000000 0.1666665000000000 0.3152396685456632 !Se0.0000004999999997 0.6666665000000004 0.3152396685456632 !Se0.4999995000000003 0.8333335000000002 0.4767849697227853 !In-0.0000005000000000 0.3333335000000000 0.4767849697227853 !In0.4999995000000003 0.8333335000000002 0.3655951960043828 !In-0.0000005000000000 0.3333335000000000 0.3655951960043828 !In

OPTCELL

100010000

POTCAR

cat Se/POTCAR In_d/POTCAR > POTCAR

4.3 静态自洽

INCAR

SYSTEM = InSe             ISTART = 0       NWRITE = 2      PREC   = AccurateENCUT  = 500GGA    = PE      NSW    = 0ISIF   = 2         ISYM   = 2       IBRION = -1      NELM   = 80        EDIFF  = 1E-05         EDIFFG = -0.01 ALGO   = Normal            LDIAG  = .TRUE.  LREAL  = .FALSE.          ISMEAR = 0         SIGMA  = 0.05 ICHARG = 2NPAR   = 4

KPOINTS

Monkhorst Pack0Gamma21 13 1.0 .0 .0

POSCAR

cp CONTCAR scf/POSCAR

4.4 能带计算

INCAR

SYSTEM = InSe             ISTART = 1       NWRITE = 2      PREC   = AccurateENCUT  = 500GGA    = PE      NSW    = 0ISIF   = 2         ISYM   = 2       IBRION = -1      NELM   = 80        EDIFF  = 1E-05         EDIFFG = -0.01 ALGO   = Normal            LDIAG  = .TRUE.  LREAL  = .FALSE.          ISMEAR = 0         SIGMA  = 0.05 ICHARG = 2LORBIT = 11LWAVE  =  .FALSE.          LCHARG = .FALSE.




    
NPAR   = 4

KPOINTS

k-points along high symmetry lines80Line-modeRec 0 0.5 0 !Y 0 0 0 !gamma 0 0 0 !gamma 0.5 0 0 !X 0.5 0 0 !X 0.5 0.5 0 !S 0.5 0.5 0 !S 0 0.5 0 !Y

4.5 有效质量计算

计算说明：对于包含了晶格周期性的有效质量的表达式如下所示：

在带入物理量进行有效质量计算时，涉及单位制问题。一般写输入文件时，长度单位为 Å；而程序输出的能带结构中，能量单位为 eV，计算起来比较繁琐。

原子单位制有两种，一种为 Hartree 原子单位制，另一种为 Rydberg 单位制。这两种单位制的区别在于，Hartree 单位制下基本物理量简单，电子电荷和质量都为1；而 Rydberg 单位制下薛定谔方程简单，系数为1。Hartree 单位制下，一个长度单位等于1L bohr = 0.5292 Å，一个能量单位1E Hartree = 27.21 eV，约化普朗克常数 ℏ=1。这样有效质量表达式中的约化普朗克常数就没了。

根据原胞基矢和正倒格子基矢间对应关系，算出倒格子基矢：

将 VASP 计算 band 时的 k 点坐标(分数坐标)转变为笛卡尔坐标：

根据两点间的距离公式，计算出各 K 点之间的距离：

4.6 求每个 K 点的位置值

根据 VASP 计算能带时各高对称点间均匀撒点，求出每个点的位置值，第一个点设为 0，本例中为 80 个点，在 excel 中进行操作。

因计算时均匀撒点 80 个，故有 79 个小间隔，对于 |YΓ| 来说，每个小间隔为 0.002975210683544，故1-80个点的坐标值都可算出，以此类推，后面的点的坐标在前面点的基础上加上间隔即可(注意：在 80 个点结束处和 81 个点开始处的值是一样的，后面的点类似)。

4.7 画出价带顶和导带底的能带

在 origin 中找出能带数据的价带顶 (VBM) 和导带底 (CBM) 的数据，把上面计算得到的 K 点路径做为 X 轴，VBM 和 CBM 作为 Y 轴，在 origin 中画图如下：

4.8 计算有效质量(以x方向电子的有效质量为例)

首先换算能量单位，由 eV 换算为原子单位制下的能量 CBM/27.21 然后选取 Γ-X 方向上以 Γ 开始的 4-8 个点的数据画能带图，如下:

用 y = a+bx+cx^2 函数拟合，操作如下：

有效质量为 1/2 C ，其中 C =2.69188，算得有效质量为 0.19 ，为电子沿 x 方向的质量。文献计算结果为 0.19，符合一致。

5、二维 InSe 载流子迁移率计算过程

本文载流子迁移率的计算依据的是形变势理论，因此需要对二维 InSe 的 x 方向和 y 方向施加应变，施加应变的范围是 -2%~2%。为了计算的准确性，计算过程中考虑了泊松效应，即对 x 轴施加应变时，固定 x 轴和 z 轴，优化 y 方向的晶格常数；对 y 轴施加应变时，固定 y 轴和 z 轴，优化 x 方向的晶格常数。

5.1 准备输入文件

首先建立 mobility-x 文件夹，然后在该文件夹下建立初始文件夹，命名为 IS，其结构目录如下：

$ tree mobilitymobility├── IS│ ├── 2_scf│ │ ├── band│ │ │ ├── INCAR│ │ │ ├── KPOINTS│ │ │ └── POTCAR│ │ ├── INCAR│ │ ├── KPOINTS│ │ └── POTCAR│ ├── INCAR│ ├── KPOINTS│ ├── OPTCELL│ ├── pbs│ ├── POSCAR│ └── POTCAR└── mobility.sh

说明：

pbs文件

#!/bin/bash#PBS -N mobility#PBS -l nodes=1:ppn=16#PBS -m abe#PBS -j n##PBS -o job.log##PBS -e job.err

理论计算干货：VASP计算二维材料的载流子迁移率

正文

请到「今天看啥」查看全文