本文只是笔者对软件操作的一些总结与记录，如有错误欢迎指正。所有更详细的内容都可以在文章最后的引用找到。

0.ZT值的浅析（thermoelectric figure of merit）

ZT值是衡量热电材料热电性能的指标，它决定了在特定温度下热电材料能量转化的最大效率

$ZT%3D%5Csigma%20S%5E2%20T%2F%5Ckappa%20$

其中σ为电导率，κ为热导率，S为塞贝克系数，T 为温度。而热导率κ又可写为电子热导 $%20%5Ckappa%20_e$ 与晶格热导 $%5Ckappa_c$ 之和：

$%5Ckappa%20%3D%20%5Ckappa%20_e%20%2B%5Ckappa_c$

因此要计算ZT值，我们需要分别计算电导率σ，塞贝克系数S，电子热导率 $%20%5Ckappa%20_e$ ，晶格热导率 $%5Ckappa_c$

本文主要介绍使用vaspkit计算电导率σ，塞贝克系数S，及使用ShengBTE计算电子热导率 $%20%5Ckappa%20_e$ ，晶格热导率 $%5Ckappa_c$ 。

1.电导率σ，塞贝克系数S的计算

这两个物理量是通过半经典的Boltzmann输运理论进行计算的，详细的理论推导可以参考文章[1][4]。

一般情况下我们可以通过vasp+BoltzTrap 程序包来计算，可参考[2][3]。但在不久前vapskit更新了1.31版本[4][5]后，我们可以完全脱离BoltzTrap 程序包，仅使用vaspkit计算处对应的物理量。

具体流程为：

准备好计算的材料对应的POSCAR。如果是二维材料可以使用vaspkit 的921或923功能对二维材料POSCAR进行标准化。
进行结构优化。
使用 vaspkit-681命令生成高密度的KPOINTS，然后进行静态计算（注意只有使用这项功能生成KPOINTS计算的结果才能继续使用vaspkit命令计算下一步，使用M-S方法自动生成K点的计算结果无法进行下一步）。
准备对应的INPUT.in文件用于输运性质计算。

#INPUT.in  仅供参考
1           # 1 for gradient method with constant relaxation time approximation
2D          # 2D for slab, 3D for bulk
-2.0  2.0   # Minimum and maximum values of chemical potential (float type, in units of eV), referenced to the Fermi energy
1000        # Number of intervals between the minimum and maximum chemical potential values
300.0       # Temperature (float type, in units of K)
1.0         # Relaxation time (float type, in units of s), Set 1.0 if you don't know the exact relaxation time
0.0         # Desired gap value in scissor correction (float type, in units of V), Zero means make no correction, available only non-spin-polarization calculation
1           # Finite difference method to calculate electron group velocity (integer type)

#关于INPUT文件中的参数，一般只需要调整Temperature 以及Relaxation time。其中Relaxation time 一般可以通过查找文献得到当前研究的结构的载流子Relaxation time/scattering time 。当我们查找不到时就需要自己计算出Relaxation time。

#基于形变势理论(DP)的Relaxation time计算,参考[10]

##原版的BoltzTrap程序中输入参数并不包括relaxation time,vaspkit中INPUT文件的relaxation time到底是什么意思我不太清楚。

##但是将relaxation time设置为1计算得到的结果除以通过其他方式得到的relaxation time后，最终的不带时间量纲的结果才是正常。这一步类似于传统的BoltzTrap程序输出结果除以通过其他方式得到的relaxation time。

5.使用vaspkit-682计算输运性质，对应的结果将以.dat文件形式保存在当前目录。

2.晶格热导率 $%5Ckappa_c$ 的计算

计算晶格热导率我们需要用到的软件包括Phonopy,Thirdorder,ShengBTE[6]。

其中Phonopy用于计算声子谱及二阶力常数，Thirdorder用于计算三阶力常数，ShengBTE用于结合前面两者的结果计算晶格热导率。可参考王宁博士在B站的视频[7][8][9]。

2.1 Phonopy计算声子谱及二阶力常数

计算声子谱及二阶力常数的具体流程如下：

对初始结构进行高精度的结构优化
这一步中INCAR的主要参数是EDIFFG，一般情况下应达到EDIFFG=-1E-8的标准。考虑优化速度，可以通过优化多次，每次优化时逐步减小EDIFFG直到EDIFFG=-1E-8的方法进行优化。高精度优化中IBRION建议设置为1，且当EDIFFG较小时建议设置ISIF=2。
使用Phonopy进行扩胞
一般情况下，扩胞后的超胞中的原子数达到100就可以了。扩胞后会生成SPOSCAR 与 POSCAR-* 等文件。前者可用于DFPT（密度泛函微扰）方法，后者应用于有限位移法。两种方法计算的结果没有区别。

phonopy -d --dim='4 4 1'  #表示x，y方向扩胞4倍，z方向不扩胞

3.1 DFPT法

DFPT法需要使用SPOSCAR进行计算（单个任务）。可参考王宁博士在B站的视频[7]。笔者在这里贴出自己的代码仅供参考。

计算前需要将原高精度优化的POSCAR 重命名为POSCAR-unitcell SPOSCAR命名为POSCAR。

#INCAR for DFPT

Global Parameters
ISTART =  0            (Read existing wavefunction; if there)
ICHARG =  2         (Non-self-consistent: GGA/LDA band structures)
LREAL  = .FALSE.       (Projection operators: automatic)
ENCUT  =  400        (Cut-off energy for plane wave basis set, in eV)
PREC   =  High       (Precision level)
ADDGRID= .TRUE.        (Increase grid; helps GGA convergence)

IALGO=38
NSW    =  1          (number of ionic steps)
ISMEAR =  0            (gaussian smearing method )
SIGMA  =  0.05         (please check the width of the smearing)
IBRION =  8            (Algorithm: 0-MD; 1-Quasi-New; 2-CG)
#最近发现对于二维材料应该取IBRION=7，因为IBRION=8考虑了对称性，而二维材料缺少了Z方向的对称性
#如果发现自己的二维体系IBION=8算出来虚频，可以试试IBRION=7
EDIFFG = -1E-07      (Ionic convergence; eV/AA)
EDIFF = 1E-07
PREC   =  Accurate     (Precision level)
# DFPT 需要IBRION =  8
#考虑范德华相互作用时IBRION = 或5或6
#考虑SOC时 需要加上以下参数
{
#LSORBIT = .TRUE.
#SAXIS = 0 0 1
#GGA_COMPAT = .F.
#ISPIN = 2
}

vasp计算完成后编写band.conf 并运行以下命令

phonopy --fc vasprun.xml
phonopy --dim='5 5 1' -c POSCAR-unitcell -p -s band.conf
bandplot  --gnuplot >a.dat

#band.conf for DFPT
ATOM_NAME = P
DIM = 5 5 1
BAND =   0.0 0.0 0.0 0.5 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0 0.0 0.0 0.0 #高对称点可由vaspkit得到
BAND_LABELS = G X S Y G
FORCE_CONSTANTS = READ
EIGENVECTORS .TRUE.
BAND_CONNECTION = .TURE.
#FC_SYMMETRY = .TRUE.