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

结构优化

声子谱和声子态密度计算

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

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

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

弛豫时间的计算

热输运性质的计算

电输运参数和ZT值的计算

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

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

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

天玑算-科研服务-独家 | （三）二维材料热电性质计算方法-能带计算，弛豫时间计算

天玑算科技，为科研助力。今天小天就为大家更新热输运性质，电输运参数和ZT值的计算。

计算教程

七、热输运性质的计算

（一）计算三阶力常数矩阵

准备输入文件：INCAR POSCAR POTCAR KPOINTS，其中POSCAR，POTCAR和KPOINTS与自洽计算一致即可，INCAR设置如下：

SYSTEM=3rd

PREC=High

ISTART=0

ICHARG=2

ISPIN=1

NELM=60;NELMIN=4;NELMDL=-3

EDIFF=1E-7

IALGO=38

ADDGRID=True

LREAL= .F.

NSW=0

IBRION=-1

EDIFFG=-1E-7

ISMEAR=0;SIGMA=0.01

LWAVE=F

LCHARG=F

执行命令： /thirdorder_vasp.py sow x y z -c (x,y,z是三个方向的扩胞倍数，c是考虑多少个临近原子的受力来计算力常数矩阵)

程序运行后的界面如上图所示，其中Automatic cutoff是考虑14个原子间受力对应的截止半径。文件夹中的输出的内容如下图：

这就以意味着需要进行828个POSCAR对应的自洽计算。通过thirdoder官网提供的脚本，生成计算所需的文件夹，脚本如下：

for i in 3RD.POSCAR*;do

  s=$(echo $i|cut -d"." -f3) &&

  d=job-$s &&

  mkdir $d &&

  cp $i $d/POSCAR &&

  cp INCAR POTCAR KPOINTS $d &&

  (cd $d && qsub runvasp.sh)

  done

这脚本可以直接复制到命令栏中，运行完成后的文件夹输出内容为：

批量提交作业，批量提交作业脚本（天河超算）：

然后执行find job* -name vasprun.xml|sort -n|/ thirdorder_vasp.py reap x y z -c 收集力常数：

运行完成后，文件夹中出现一个新的文件： FORCE_CONSTANTS_3RD

至此，三阶力常数矩阵就算完成了。

其中需要注意的是：第一，截止半径需要测试，经验上我们需要取到10以上的截止半径才可能得到收敛的计算结果。第二，计算三阶力常数矩阵的扩胞前的POSCAR需要与计算二阶力常数保持一致。第三，如果读取到某个job时报错，请重新计算这个job。

接下来，我们就可以用FORCE_CONSTANTS_2ND  FORCE_CONSTANTS_3RD结合控制文件CONTROL来计算晶格热导率了

CONTROL文件可以参考如下设置：

&allocations

   nelements=2   #元素种类

   natoms=16      #原子数量

   ngrid(:)=21 21 1  #计算网格取点，需要进行严格测试

&end

&crystal

   lfactor=0.100000  #晶格缩放比例，ShenBTE采用nm作为单位，因此与vasp差一个数量级

   lattvec(:,1)=10.1939042571769480    0.0000000000000000   -0.0005202257779721

   lattvec(:,2)=0.0000000000000000    5.8835567629900831    0.0000000000000000

   lattvec(:,3)=-0.0016161778037250    0.0000000000000000   24.2394661807105862

   elements= "Au" "S"  #元素符号

   types= 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2   #对应不同元素

 positions(:,1)=0.0792887906745038  0.7628181349429907  0.5477417153337826

 positions(:,2)=0.5792887900793176  0.2628181570448227  0.5477417107462345

 positions(:,3)=0.3420845204951478  0.5000001343163234  0.5479707485969261

 positions(:,4)=0.8420845388897342  0.0000001380045518  0.5479707590273382

 positions(:,5)=0.0792887988980483  0.2371820787407889  0.5477416212454472

 positions(:,6)=0.5792887903691906  0.7371820856979016  0.5477416090636752

 positions(:,7)= 0.4207111952601045  0.7371818360856525  0.4522582828547815

 positions(:,8)=0.9207112108774353  0.2371818574692190  0.4522582909979565

 positions(:,9)=0.4207111998038618  0.2628179122832144  0.4522583851953175

 positions(:,10)=0.9207111996088315  0.7628179378276924  0.4522583841178008

 positions(:,11)=0.1579154828371944  0.9999998672025896  0.4520292560634112

 positions(:,12)= 0.6579154646203824  0.4999998598272342  0.4520292358640028

 positions(:,13)=0.3333765322103900  0.0000000837505425  0.3901738436731563

 positions(:,14)=0.8333765550624702  0.5000000766771890  0.3901738421682481

 positions(:,15)=0.1666234737495292  0.4999999208298789  0.6098261606528128

 positions(:,16)=0.6666234565638564  0.9999999192994125  0.6098261543990704

   scell(:)= 2 4 1     #二阶力常数计算扩胞倍数

&end

¶meters

   T=300    #温度控制

   scalebroad=0.1   #控制高斯展宽，默认取1，手册中说改小后对计算精度没有影响，但是最近的一些研究发现对计算精度影响较大。收敛测试反应，取到4才收敛

&end

&flags

   nonanalytic=.TRUE.  #控制参数具体参考手册

   nanowires=.FALSE.  #控制参数具体参考手册

&end

此时，执行ShengBTE程序计算，完成后，在文件夹中应该包括以下内容：

T300K文件夹中应该包含以下内容，

第二步，数据提取：打开BBTE.kappa_tensor文件，正确的显示为：

其中的78为使用迭代法（ShengBTE默认方法）计算晶格热导率的总迭代数。后面为3X3的9个方向的晶格热导率矩阵，顺序为xx，xy，xz，yx，yy，yz，zx，zy，zz。其中的xx，yy和zz对应正交晶胞中的abc三个方向。

注意：此时计算得到的晶格热导率不是真实值，需要考虑有效厚度。故做进一步处理：将所得值÷有效厚度×真空层厚度。有效厚度为原子层厚度+最外层原子的范德华半径，原子的范德华半径可以在网上查到。

得到准确的晶格热导率之后，我们还需要分析声子群速度、平均自由程等参数来进一步分析材料的热输运性质：用ShengBTE_TC.py生成所需要的热输运输出文件，将ShengBTE_TC.py拷贝到文件夹下，执行命令：python ShengBTE_TC.py（该文件作者为大连理工张晓亮，苑昆鹏）

得到如下输出文件:vg_freq.dat lifetime_freq.dat BTE.qpoints BTE.gruneisen thermal_conductivity_accumulated_mfpx.dat thermal_conductivity_accumulated_mfpy.dat  BTE.omega   BTE.w_anharmonic (T300K文件下)

vg_freq.dat: 声子群速度，可直接导入Origin中（如图10），第一列是频率，单位是THz；第二，三，四列分别代表x, y, z方向的group velocity（要转换成绝对值，进行abs操作）, 单位是Å/ps, 要换算成km/s，所以要除以10。

图10

lifetime_freq.dat：声子弛豫时间（不分方向， 如图11）

图11

BTE.qpoints：对应路径的点的信息（如图12）

图12

比如，Γ—X方向的路径为(0, 0, 0)—(0.5, 0, 0), 从图中可以看出对应的点为1-26，注意：如果中间有隔断的不属于此路径的点，需将这些点删去，并记住所需点的序号，以便后续分方向的热输运参数的求解。

BTE.w_anharmonic：热输运参数x轴信息（第一列，如图13）

图13

BTE.omega：声子弛豫时间（y轴数据，如图14）

图14

分方向时：找到需要路径的点对应的序号所在的行，将其余行删掉，然后将最右边一列剪切到右边倒数第二列下面，然后再整体剪切到倒数第三列下面。以此类推最后将全部数据剪切到第一列下面，组成的这列就是y轴的数据，与图12中对应序号行的x轴的数据一起plot，得到所需方向的声子弛豫时间的图形。

BTE.gruneisen：格林艾森参数（y轴数据，如图15）

图15

不分方向时：将最右边一列剪切到右边倒数第二列下面，然后再整体剪切到倒数第三列下面。以此类推最后将全部数据剪切到第一列下面，组成的这列就是y轴的数据，与图12中x轴的数据一起plot，得到所需方向的声子弛豫时间的图形。

分方向时：处理方式同声子弛豫时间。

thermal_conductivity_accumulated_mfpx.dat：x方向的MFP（如图16）

图16

thermal_conductivity_accumulated_mfpy.dat：y方向的MFP（如图17）

图17

八、电输运参数和ZT值的计算

将前面高精度自洽生成的DOSCAR,POSCAR,OUTCAR,EIGENVAL和vaspkit拷贝到新文件夹中，并将此文件夹重命名为case。在case文件夹下运行vaspkit：./vaspkit，选择92 （VASP2BoltzTraP，如图18，但注意：有的版本是选择73：VASP2BoltzTraP Interface）。

图18

得到三个case文件：case.energy  case.intrans  case.struct。case.intrans里面的参数设置如下：

主要需要修改的参数为energygrid，Tmax和temperaturegrid。将能量网格加密（一般0.0001足矣），最大温度和温度梯度按照自己的需要设置。运行boltztrap，得到case.condtens文件。接下来我们以300K为例，描述各向异性的二维材料在n型和p型掺杂下x和y方向的电输运参数（Seebeck,电导率，电子热导率，功率因子）和ZT值的计算过程：

第一步：用excel表格打开case.condtens文件，在表格中点击数据的第一列并进行数据分列，如下图：

第二步：打开Origin 9在里面建立新的4个工作表格，命名为P-X、N-X、P-Y、N-Y分别代表x方向p型掺杂，x方向n型掺杂，y方向p型掺杂和y方向n型掺杂下的数据。

第三步：在Excel中进行温度筛选，筛选出300K温度下的数据，如下图：

第四步：把Excel中的300K对应的浓度（N）复制粘贴到Origin（注意区分正负，浓度为负值是N型之后需放在对应的N-X或N-Y表格中）。

第五步：把Excel中的cond、seebeck系数x和y方向的数据依次粘贴到Origin中P-X和P-Y的表格中(cond第一列是X方向，第五列是Y方向，seebeck系数同理)，全部粘贴到Origin中后，将浓度负值相应的cond、seebeck系数依次粘贴到Origin中N-X和N-Y表格中。

第六步：为了方便整理，在Origin中重命名每列数据的名字。

第七步：在Origin中对每列数据进行单位换算，也就是把BoltzTraP输出的单位换算成国际单位或标准单位，具体参考下表：

第八步：对浓度的单位换算。因为BoltzTraP中得到的浓度单位是e/uc，即每个晶胞中的电子数目，所以我们需要得到原胞的体积，又因为我们计算的是二维材料，所以只需要知道晶胞在c方向的截面面积即可，具体操作为：选中浓度值一列并执行操作（Ctrl+Q），在弹出的命令框内输入col(a)/数值*1e16并点击‘OK’(注意a是指浓度值所在列，浓度为负时要进行绝对值化处理（在命令框最前面加abs（）），而数值是晶胞真空层方向横截面的面积a×b，从POSCAR文件中可以算出，乘1016是因为要把Å2换算成cm2)，得到的浓度单位为e/cm2。

第九步：同理，对cond执行操作（Ctrl+Q）：col（b）*数值(注意b是指电导率所在列，数值是前面用DP理论求得的弛豫时间。（注意：在这里弛豫时间单位应该换算成秒，得到的电导率单位为1/(Ω·m)。

第十步：对Seebeck执行操作（Ctrl+Q）：abs（col（c））*1e6（c是指Seebeck系数所在列），得到的单位为μV/K。

第十一步：求出功率因子pf值（在Seebeck后面插入空白列命名为pf值，选中这一列，执行操作（Ctrl+Q）：(col（c）/1e6)^2*col（b）*1e3，得到的单位为mW/mK2。

第十二步：插入空白列命名为图片，根据公式图片=LσT执行操作（Ctrl+Q）：2.45*1e-8*col(b)*300（此处300是指温度300K，L取2.45×10-8），得到的单位为W/mK。

第十三步：插入空白列，命名为ZT，根据公式 图片 执行操作：Ctrl加Q打开公式命令框，输入col(d)/1e3*300/(col(e)+数值)（此处d为pf值所在列，e为ke值所在列，数值为之前求的晶格热导率的值，300为温度），得到ZT值数据。

第十四步：以浓度为横坐标，Seebeck系数、电导率、PF因子、ZT值为纵坐标，分别画图，将横纵坐标调整到合适的大小，即可得到x或y方向，n型或p型掺杂下的最优ZT值，如下图：