声子性质、热输运和热电性质

本计算实例演示如何使用QuantumATK的声子模块。计算中使用的模型很小，这样保证计算实例能够快速完成。第一部分计算声子的能带和态密度，并简单介绍QuantumATK声子模块的工作原理。第二部分对一个器件体系进行了电子和声子输运计算，并由此得到热电品质因数ZT。

本实例主要使用Virtual NanoLab (VNL) 进行计算设置和结果分析，用户可以参考VNL的教程熟悉相关的软件操作。全部计算主要使用ATK-Classical和ATK-SE两个模块，具体参数设置方法需要参考QuantumATK的参考手册。

本实例计算需要ATK-SE的license。用户可以运行有时间限制的试用版或联系费米科技购买正式版权。

启动VNL，建立一个新的Project，打开新建的Project。从主窗口左侧工具栏启动Builder。在Builder窗口中点击Add→From Plugin→Naoribbon。

• 去掉Passivate dangling bonds with hydrogen前的勾选
• 选择Armchair edge
• 选择8原子宽度

点击Build按钮，创建结构。

为了让声子模块能够确定这是个一维的结构，结构的两边都要有最少7Å的真空。为此，打开Bulk Tools→Lattice Parameters…。

• A-x方向增加到20Å
• B-y方向增加到30Å

打开Coordinate Tools→Center，点击Apply，将纳米带移到胞的中心。

使用窗口右下角的箭头图标中的Send To选项将结构传送到Script Generator。

接下来设置计算，用ATK-Classical模块中的Tersoff势计算结构优化和声子能带和态密度的计算。

在Script Generator窗口中

• 添加New Calculator
• 添加Optimization→OptimizeGeometry
• 添加Analysis→PhononBandstructure
• 添加Analysis→PhononDensityOfStates
• 将输出文件（output file）改成armchair.nc

• 双击New Calculator，选择ATK-Classical和Tersoff势
• 双击OptmizeGeometry，设置maximum force=0.01eV/Å
• 双击PhononBandstructure，设置Points pr. segment=100
• 双击PhononDensityOfStates，不勾选sync，设置nC=51

将脚本传送至Job Manager，提交计算。

计算生成的结果文件armchair.nc可以在VNL主窗口中的Project Files下看到，文件中包含了PhononBandstructure对象。选中该图标，点击右侧分析工具栏中的Show 2D Plot…。将图的低频部分放大如下图。

负能带说明该结构并不在能量曲面的稳定点上，进一步的结构弛豫可以进一步的降低能量。前面的操作中仅优化了原子受力但是没有优化张力，因此结构有沿垂直石墨烯方向弛豫的倾向。

重新打开Script Generator，打开OptmizeGeometry。

• 将Maximum Stress设置为0.0001 eV/ų
• 去掉constrain cell下的z前面的勾选

点击OK，将脚本提交到Job Manager，运行计算。计算结束重新画出的能带结构中负能带消失了。

armchair.nc文件还包括声子态密度（DOS）结果，包含优化张力前后的两个结果。下图的结构弛豫之后的结果。

从声子态密度还可以使用entropy()函数计算体系的熵。下面脚本用于计算体系的熵并得到自由能。

entropy.py

# Read the last configuration in the output file
configuration = nlread('armchair.nc', BulkConfiguration)[-1]
# Read the last phonon DOS in the output file
phonon_dos = nlread('armchair.nc', PhononDensityOfStates)[-1]
# Calculate the total energy
total_energy = TotalEnergy(configuration)
e = total_energy.evaluate()
 
# Define the temperature
t = 300*Kelvin
 
# Calculate the entropy
s = phonon_dos.entropy(t)
 
# Print results
print 'Energy  = ', e
print 'Entropy = ', s
print 'Free energy at 300K = ', e - t*s

将脚本拖到Job Manager，得到结果。

Energy  =  -124.297968287 eV
Entropy =  3.33141143492 meV/K
Free energy at 300K =  -125.297391718 eV

声子从体系的动力学矩阵计算得到。动力学矩阵是体系能量的二阶导数，也是力的一阶导数。力的导数采用有限位移的方法进行计算，这种方法需要将体系沿每个自由度进行微小位移，也叫冻结声子法。

在进行声子计算时，得到的动力学矩阵会保存在体系对应的configuration数据对象中。

动力学矩阵的计算由DynamicalMatrixParameters对象控制。

dynamical_matrix_parameters = DynamicalMatrixParameters(repeats = (1,1,5),
atomic_displacement = 0.01*Ang)
calculator =
Tersoff(dynamical_matrix_parameters=dynamical_matrix_parameters)

为了考虑临近晶胞对动力学矩阵的贡献，计算将对体系进行重复，参数可控。如果重复参数设为automatic，体系自动确定重复单元数，输出文件中会给出

| Phonon: Automatically detected repetitions = [1 1 5]

atomic_displacement参数控制冻结声子时的原子位移大小。

使用ATK-DFT和ATK-SE也可以计算声子。此时要注意k点设置会用在经过重复的晶胞上。例如上面体系计算使用(1,1,21)的k点进行结构优化，而用于动力学计算的5倍重复单胞则可以设置较少的k点，这样更可以节约计算时间，下面脚本举例说明如何修改k点。

# Change k-point sampling of a calculator
calculator = bulk_configuration.calculator()
numerical_accuracy_parameters = calculator.numericalAccuracyParameters()
new_numerical_accuracy_parameters = numerical_accuracy_parameters(k_point_sampling = (1, 1, 5))
new_calculator = calculator(numerical_accuracy_parameters=new_numerical_accuracy_parameters)
bulk_configuration.setCalculator(new_calculator)

接下来我们将计算armchair纳米带的声子透射系数，进而得到声子的热导。将这部分结果与电子透射系数（可以得到电导、Peltier系数和电子热导）结合起来，可以计算热电的品质因数（figure of merit）ZT。计算方法在文献[1]中描述。

第一步是从上一部分得到的经过弛豫的armchair纳米带结构建立一个器件模型。在VNL主窗口中选择armchair.nc中的最后一个configuration，拖动到Builder里。 在建模工具窗口：

• Bulk Tools→Repeat:在C方向重复5次；
• Device Tools→Device From Bulk:使用默认设置，点击OK

将结构传送到Script Generator。

接下来设置声子和电子的透射计算。声子计算使用Tersoff势，电子透射计算则采用使用$\pi$紧束缚模型。

• 添加New Calculator
• 添加Analysis→PhononTransmissionSpectrum
• 添加New Calculator
• 添加Analysis→TransmissionSpectrum
• 将输出文件（output file）改成armchair_device.nc

按以下步骤设置计算参数：

• 双击打开New Calculator，选择ATK-Classical计算器和Tersoff势。
• 双击打开PhononTransmissionSepctrum，将energy range设置为0 eV到0.5 eV
• 双击打开第二个New Calculator，选择ATK-SE：Slater-Koster计算器，选择Hancock.C ppPi基组
• 不用设置TransmissionSpectrum

将脚本传送至Job Manager，运行计算。

计算结束后，在主窗口找到armchair_device.nc文件，你可以选择声子或是电子的透射系数谱，用Transmission Analyzer…或Show 2D plot进行分析。

接下来你可以计算在两电极温度差（$dT$）或电压差（$dU$）下的输运性质的线性响应系数： 电导 $$ G_e =\left.\frac{dI}{dU}\right|_{dT=0} $$ Peltier系数 $$ \Pi=\left.\frac{I_Q}{I}\right|_{dT=0} $$ Seebeck系数 $$ S=-\left.\frac{dU}{dT}\right|_{I=0}=\frac{\Pi}{T} $$ 电子和声子热输运系数$\kappa_e$和$\kappa_{ph}$ $$ \kappa=\left.\frac{dI_Q}{dT}\right|_{I=0} $$ 其中$I_Q=dQ/dT$是热流。

有了这些系数，就可以计算热电品质因数，热电品质因数描述了热电转换的效率。 $$ ZT=\frac{G_eS^2T}{\kappa} $$ 这些线性响应系数可以用Transport Coefficient插件计算。同时选中PhononTransmissionSpectrum和ElectronTransmissionSpectrum，点击Calculate按钮即可。

未掺杂的纳米带，费米能级位于能隙中间。为了得到大的Peltier系数，费米能级必须位于能带的边缘，这可以用掺杂实现。这个计算插件可以手动移动费米能级，下图显示将费米能级移动了$\Delta E_F=0.08 eV$的结果，对应了纳米带的n-型掺杂。

1. T. Markussen, A.-P. Jauho , M. Brandbyge , Phys. Rev. Lett., 103, 055502, 2009.