这是本文档旧的修订版!
本文用简单的实例介绍使用ATK研究电子输运的基本概念,计算设置应该注意的事项和对体系输运性质进行的常用分析方法。
为了快速完成计算,这里用到的体系是Zn-ZnO-Zn,使用ATK-SE半经验方法进行计算。体系的结构见下图:
实例计算中使用紧束缚近似(TB)的半经验方法,结合非平衡态格林函数方法(NEGF)研究上述器件结构的电子输运特性。理论方法的详细描述见:参考手册。
这里主要介绍器件电子输运计算的核心概念与在VNL-ATK中可以进行的常用计算:
与通常的计算模拟中常见的块体材料的周期性体系(bulk)和分子体系的孤立结构(molecule)不同,用于电子输运计算的器件体系(device)比较特殊。器件结构由左、右两个电极(electrode,分别是周期性体系)和中间区域(central region,分子体系)构成。在晶体管器件中,左、右电极也称为源、漏电极。中间区域则是一个有限的区域,也称为散射区(scattering region),由于载流子从一个电极到另一个电极传输时会受到中间区域的散射,因此,中间区域是决定器件特性的关键区域。 中间区域除了可能存在的不同于电极部分的材料之外,还包括左、右电极部分的一个完整的、严格的重复(electrode extension),这在使用VNL构建器件结构和ATK的整个计算过程中是必须满足也是默认满足的条件。
中间区域通常包含一些非均质的结构(界面、分子、缺陷等等),这导致中间区域的电子密度和自洽计算中的电子感受到的有效势发生变化,由此可能进一步导致中间区域的结构发生弛豫(当然,弛豫不包括电极扩展部分)。这种非均质的效应导致电势(电压)在中间区域发生降落。
完美周期结构的电子输运情况只能在0偏压下研究,这是因为在完美周期体系中由于没有散射中心导致没有电压降。对于此类体系可以使用块体材料来研究0偏压的输运情况。
NEGF理论方法假定中间区域的电势能够与电极部分平滑的过渡,也就是说在电极和中间区域之间不存在界面的散射。(注意:此处提到的界面指的是下图中蓝色方框的电极和红色方框的中间区域之间;在中间区域内部的Zn和ZnO之间的界面当然存在合理的散射,这部分可能是研究的主要对象)。为了确保这种平滑的过渡,中间区域要保留足够多的的电极材料层(含电极扩展部分),这部分称作屏蔽区域(下图中的中间区域的灰色部分)。如果屏蔽区域厚度不够,可能会导致自洽计算收敛变慢、结果不准确等问题。
上图中还显示了由于外加偏压(0.5V)导致的电势分布。左右电极的电势差为0.5eV=0.018Ha。图中明显看出,所有的电势降落都发生在ZnO部分,而金属Zn部分的电势则相当的平坦。
本征半导体材料电极中的屏蔽长度可能非常长,因此需要非常长的中间区域才能使电势平滑过渡到块体电极部分。不过在实际情况中,屏蔽距离会随掺杂浓度的提高而变短。
打开VNL,创建一个新的Project。一般情况下,可以使用Builder构建新的结构,这里可以直接下载结构文件,放在打开的Project文件夹下。 将脚本拖动到Builder图标上即可打开Builder窗口看到结构。
ATK的器件计算通常分为多个步骤。 首先是将电极和中间区域作为周期性的块体材料体系进行自洽电子密度计算。 随后开始的器件计算首先将中间区域的电子密度作为初始猜测,使两个电极和中间区域在接触的边界匹配。 最后计算得到中间区域的自洽电子结构。
上述步骤中,电极块体区域的收敛性很重要,既要保证电极长度足够,也要保证数值精度足够。 这一节介绍两个重要的数值精度控制参数:实空间密度网格截断能(mesh cut-off)和横向(A向和B向)的k点数。
下面的收敛性测试只针对左电极,因为这个实例体系的左右电极完全相同。 对于左右电极不全同的情况,要对左右电极分别进行测试,并在最后的器件计算里选取精度更高的一组参数作为最终计算的输入。
使用上方工具条中的器件分割工具(splitter tool)将器件结构分割成两个电极和中间区域,三个产生的结构出现在下部的结构浏览区域(Stash)。
接下来选中左电极结构,点击窗口右下角箭头发送到Script Generator,用不同的mesh cut-off设置一系列块体材料计算,以便研究数值收敛性。
为方便起见,对mesh cut-off和k点进行测试的脚本可以直接在这里下载。接下来只要运行两个脚本即可。
第一个要测试收敛性的参数是实空间网格截断能,这个参数控制实空间网格(用来投影静电势和电子密度)的精细程度。这个精细程度之所以表述为一个能量是因为网格格点的间距可以由下面公式给出: \[ \Delta x = \frac{\pi \hbar}{\sqrt{2mE}} \]
$m$是电子的质量。
脚本 mesh_loop.py 计算电极块体的基态,mesh cut-off取值从 2.5 Ha 到 80 Ha。所有控制数值精度的其他参数都保持不变。电极的费米能级(化学势)作为收敛标准的判据量。
# ------------------------------------------------------------- # Bulk Configuration # ------------------------------------------------------------- # Set up lattice vector_a = [1.3325, -2.30795770109, 0.0]*Angstrom vector_b = [1.3325, 2.30795770109, 0.0]*Angstrom vector_c = [0.0, 0.0, 4.947]*Angstrom lattice = UnitCell(vector_a, vector_b, vector_c) # Define elements elements = [Zinc, Zinc] # Define coordinates fractional_coordinates = [[ 0.333333329801, 0.666666670199, 0.25 ], [ 0.666666670199, 0.333333329801, 0.75 ]] # Set up configuration bulk_configuration = BulkConfiguration( bravais_lattice=lattice, elements=elements, fractional_coordinates=fractional_coordinates ) # ------------------------------------------------------------- # Calculator # ------------------------------------------------------------- #---------------------------------------- # Basis Set #---------------------------------------- basis_set = DFTBDirectory("cp2k/scc/") #---------------------------------------- # Pair Potentials #---------------------------------------- pair_potentials = DFTBDirectory("cp2k/scc/") # List of values for the grid mesh cutoff cutoffs = [2.5, 5, 10, 20, 30, 40, 60, 80] # List to hold the chemical potential for each calculation. chemical_potentials = [] # Loop over the grid mesh cutoffs for cutoff in cutoffs: numerical_accuracy_parameters = NumericalAccuracyParameters( k_point_sampling=(7, 7, 51), density_mesh_cutoff=cutoff*Hartree ) iteration_control_parameters = IterationControlParameters() calculator = SlaterKosterCalculator( basis_set=basis_set, pair_potentials=pair_potentials, numerical_accuracy_parameters=numerical_accuracy_parameters, iteration_control_parameters=iteration_control_parameters, ) bulk_configuration.setCalculator(calculator) bulk_configuration.update() # ------------------------------------------------------------- # Chemical Potential # ------------------------------------------------------------- chemical_potential = ChemicalPotential(bulk_configuration) value = chemical_potential.evaluate().inUnitsOf(eV) chemical_potentials.append(value) # Plot data. import pylab pylab.plot(cutoffs, chemical_potentials, 'ro-') pylab.xlabel('Grid mesh cut-off (Ha)') pylab.ylabel('Fermi level (eV)') # Show the plot. pylab.show()
脚本的内容很简单,对mesh cut-off的循环中的每个值进行一次计算,并最终得对得到的数据作图。下载脚本并在Job Manager上运行,或用命令行运行:
$ atkpython mesh_loop.py > mesh_loop.log
计算只需要很短的时间,结果作图就会显示出来(如下图)。很明显,默认的10Ha的mesh cut-off对于紧束缚近似计算是足够的,因为继续增加mesh cut-off对费米能级的影响小于10meV,但是却会大大增加计算量。因此后续计算就采用10Ha作为理想的设置。
注意:此处费米能级对mesh cut-off的变化并不敏感,这是紧束缚近似的特点。使用DFT计算时情况并不相同,计算量最小的赝势对应的默认mesh cut-off为75Ha。
与上面的方式相同,可以继续研究电极费米能级对横向(A向和B向)的k点数的依赖关系。这里A向和B向与电流方向(C向)一定是垂直的。
在器件计算中,输运方向的k点数总是必须很大($n_{C}=51$)。其他计算精度参数都保持不变,mesh cut-off保持为10Ha。
使用脚本kpts_loop.py
可以直接进行收敛性测试,此脚本和mesh cut-off测试脚本几乎一样,只是循环变量编程横向k点数,从1到11。
# ------------------------------------------------------------- # Bulk Configuration # ------------------------------------------------------------- # Set up lattice vector_a = [1.3325, -2.30795770109, 0.0]*Angstrom vector_b = [1.3325, 2.30795770109, 0.0]*Angstrom vector_c = [0.0, 0.0, 4.947]*Angstrom lattice = UnitCell(vector_a, vector_b, vector_c) # Define elements elements = [Zinc, Zinc] # Define coordinates fractional_coordinates = [[ 0.333333329801, 0.666666670199, 0.25 ], [ 0.666666670199, 0.333333329801, 0.75 ]] # Set up configuration bulk_configuration = BulkConfiguration( bravais_lattice=lattice, elements=elements, fractional_coordinates=fractional_coordinates ) # ------------------------------------------------------------- # Calculator # ------------------------------------------------------------- #---------------------------------------- # Basis Set #---------------------------------------- basis_set = DFTBDirectory("cp2k/scc/") #---------------------------------------- # Pair Potentials #---------------------------------------- pair_potentials = DFTBDirectory("cp2k/scc/") # List of values for the k-point grid. k_list = numpy.arange(1,13,2) # List to hold the chemical potential for each calculation. chemical_potentials = [] # Loop over the grid mesh cutoffs for k in k_list: numerical_accuracy_parameters = NumericalAccuracyParameters( k_point_sampling=(k, k, 51), ) iteration_control_parameters = IterationControlParameters() calculator = SlaterKosterCalculator( basis_set=basis_set, pair_potentials=pair_potentials, numerical_accuracy_parameters=numerical_accuracy_parameters, iteration_control_parameters=iteration_control_parameters, ) bulk_configuration.setCalculator(calculator) bulk_configuration.update() # ------------------------------------------------------------- # Chemical Potential # ------------------------------------------------------------- chemical_potential = ChemicalPotential(bulk_configuration) value = chemical_potential.evaluate().inUnitsOf(eV) chemical_potentials.append(value) # Plot data. import pylab pylab.plot(k_list, chemical_potentials, 'ro-') pylab.xlabel('k-points along A and B directions') pylab.ylabel('Fermi level (eV)') pylab.xticks(k_list) # Show the plot. pylab.show()
下载并运行上述脚本,就可以得到费米能级对k点的依赖关系(下图)。从图中很容易看出,A、B方向1个k点是远远不够描述Zn电极块体的,至少要5个以上的横向k点就可以得到很好的结果。
在确定了电极计算的数值精度参数(mesh cut-off和k点)之后,准备计算的最后一步是用Electrode Validator插件工具检查电极结构。这个插件可以帮你检查电极结构是否合适(晶轴C与坐标轴Z平行、AB平面与C轴垂直、电极在C方向是否足够长)。
首先需要对电极进行基态计算,步骤如下:
electrode.nc
。脚本应如下:
# ------------------------------------------------------------- # Bulk Configuration # ------------------------------------------------------------- # Set up lattice vector_a = [1.3325, -2.30795770109, 0.0]*Angstrom vector_b = [1.3325, 2.30795770109, 0.0]*Angstrom vector_c = [0.0, 0.0, 4.947]*Angstrom lattice = UnitCell(vector_a, vector_b, vector_c) # Define elements elements = [Zinc, Zinc] # Define coordinates fractional_coordinates = [[ 0.333333329801, 0.666666670199, 0.25 ], [ 0.666666670199, 0.333333329801, 0.75 ]] # Set up configuration bulk_configuration = BulkConfiguration( bravais_lattice=lattice, elements=elements, fractional_coordinates=fractional_coordinates ) # ------------------------------------------------------------- # Calculator # ------------------------------------------------------------- #---------------------------------------- # Basis Set #---------------------------------------- basis_set = DFTBDirectory("cp2k/scc/") #---------------------------------------- # Pair Potentials #---------------------------------------- pair_potentials = DFTBDirectory("cp2k/scc/") k_point_sampling = MonkhorstPackGrid( na=7, nb=7, nc=7, ) numerical_accuracy_parameters = NumericalAccuracyParameters( k_point_sampling=k_point_sampling, ) iteration_control_parameters = IterationControlParameters() calculator = SlaterKosterCalculator( basis_set=basis_set, pair_potentials=pair_potentials, numerical_accuracy_parameters=numerical_accuracy_parameters, iteration_control_parameters=iteration_control_parameters, ) bulk_configuration.setCalculator(calculator) nlprint(bulk_configuration) bulk_configuration.update() nlsave('electrode.nc', bulk_configuration)
将脚本保存为electrode.py
(或直接从上面下载),运行脚本。结果很快出现在LabFloor里。
选中electrode.nc
文件里的块体结构图标,展开右侧Electrode Validator,点击Check。
说明:
C方向的电极长度是否足够取决于基组的选择。Electrode Validator的判断标准是哈密顿量的残余矩阵元不能超过$10^{-10}$ eV。
这个标准的设置是为了确保电极在器件中的C方向只有最近邻的相互作用。为了将电极描述为半无限电极必须如此,否则电极在C方向将成为完整周期性体系。
如果电极过短,也可以继续进行器件计算,但是很有可能导致计算无法收敛或结果不准确,这要看残余矩阵元的大小。
在确定了电极结构适合进行器件计算之后,就可以开始真正的输运计算。在Builder里,选中Zn-ZnO-Zn器件结构(zn-zno-zn.py
),发送至Script Generator。
在Script Generator里,按以下步骤设置Zn-ZnO-Zn体系的电子输运谱计算。
zero-bias.nc
接下来双极添加的New Calculator,设置计算的方法:
下一步,编辑添加的TransmissionSpectrum:
说明:
至此,计算脚本设置完毕,将其保存为zero-bias.py
,或者在这里下载。
使用Job Manager或命令行方法可以运行脚本。这个计算只需要几分钟。
计算结束后,要浏览一下计算作业的log文件,可以看到器件计算的整个过程。大体上计算是按以下顺序完成的:
如果log没有错误和警告信息,说明成功的进行了一个器件的电子输运计算。接下来就可以进行结果分析和可视化了。
计算输出的结果会出现在LabFloor里,找到结果文件zero-bias.nc
,可以看到里面包含两个数据对象:
选中TransmissionSpectrum,点击右侧工具栏的Transmission Analyzer。
在打开的透射分析工具窗口的左侧是透射谱(Spectrum),即k点平均透射系数随能量的变化。左右两电极的费米能级平均值设为零,在图中用水平虚线标出。图中可以看出,在-2.5eV到2.5eV很大的区间里的透射概率值都很小,这与半导体ZnO的带隙对应。
在透射曲线上用红色圆点标出的透射系数和能量在图的下部给出,用鼠标可以在$T(E)$曲线上移动、重新选择要分析的点。左上角的Curves菜单可以用来选择不同的自旋方向的透射。对于非共线自旋情况,还可以选择自旋的X、Y、Z分量。
在(圆点标出的)特定能量处的总透射由取样的k点上的透射系数的平均求得。右侧的彩色等值图(Coefficients)显示了透射系数与倒空间矢量$k_A$、$k_B$的关系。
在上图中可以看出,Zn-ZnO-Zn体系在费米能级处的透射主要由$\Gamma$-点($k_A=k_B=0$)附近贡献。然而在其他能量处(比如下图的E=-3.16eV处),情况则完全不同。
增加透射谱计算设置时的k点数,可以提高彩色等值图的精度,也就可以提高透射谱计算的准确性。很多器件计算中都需要设置较大的k点数,以获得准确收敛的透射谱。
点击彩色等值图中的某一k点,所选点数值就显示在下方的选择(Selection)框里。对指定的自旋和k点处的透射系数,可以求算本征值和本征态:
经过短时间的计算后,所选能量、k点和自旋方向的透射本征值将会列出在白色框中。
提示:
如果透射来自于几个透射通道的贡献,则都会显示在此。每个本征值都在[0,1]之间,但是本征值的总和,即k分辨透射值$T(E,k)$,可以大于1,因此总透射$T(E)$当然也一样。
接下来可以计算并显示每个本征值对应的透射本征态:
Viewer将显示Zn-ZnO-Zn器件体系结构和中间区域的透射本征态的截面图。
接下来可以进一步调整本征态截面显示属性,以便选择更好的色调更清楚的显示透射本征态的细节:
截面图显示情况应该如下图所示。本征态对应了从左电极到右电极的散射态。因此,本征态在左侧总是振幅较大。较高透射本征值对应的本征态在散射区的右侧也有比较大的振幅,这反映了左侧的入射态有较大的透射概率穿过中间区域进入右电极。
在散射区左侧的透射本征态除了包括入射的部分(右行),还包括了反射态(左行)。这两部分会有一定的迭加干涉效应。从上图中可以看到左侧第一个Zn原子上的本征态波幅较小,这是由于入射和反射的态的相消干涉造成的。在ZnO的右侧,本征态只有透射(右行)部分,因此观察不到类似的干涉效应。
你可以再尝试画出费米能级(能量$E = 0 eV$)、$k_A=k_B=0$处的透射本征态。调整作图属性之后可以得到如下图所示。这时的本征态的波幅在散射区的中部到右部几乎消失为零,对应于很小的透射值,也就是说只有很小的散射本征态透过,而大部分都被Zn-ZnO界面反射。
提示: 如果采用等值面图的方式对透射本征态作图,等值面对应了波函数的波幅,颜色则对应了波函数的相位。