这是本文档旧的修订版！
使用ATK研究电子输运
本文用简单的实例介绍使用ATK研究电子输运的基本概念，计算设置应该注意的事项和对体系输运性质进行的常用分析方法。
为了快速完成计算，这里用到的体系是Zn-ZnO-Zn，使用ATK-SE半经验方法进行计算。体系的结构见下图：
实例简述
实例计算中使用紧束缚近似（TB）的半经验方法，结合非平衡态格林函数方法（NEGF）研究上述器件结构的电子输运特性。理论方法的详细描述见：参考手册。
这里主要介绍器件电子输运计算的核心概念与在VNL-ATK中可以进行的常用计算：
- 器件的基本结构
- 选择合适的计算精度
- 检查块体材料是否适合用作电极
- 计算0偏压和有限偏压下的透射谱
- 使用透射分析工具分析透射谱
- 计算透射本征态
- 计算投影器件态密度
- 计算器件内电势分布（电压降）
- 计算、分析IV曲线
器件体系的几何结构
与通常的计算模拟中常见的块体材料的周期性体系（bulk）和分子体系的孤立结构（molecule）不同，用于电子输运计算的器件体系（device）比较特殊。器件结构由左、右两个电极（electrode，分别是周期性体系）和中间区域（central region，分子体系）构成。在晶体管器件中，左、右电极也称为源、漏电极。中间区域则是一个有限的区域，也称为散射区（scattering region），由于载流子从一个电极到另一个电极传输时会受到中间区域的散射，因此，中间区域是决定器件特性的关键区域。 中间区域除了可能存在的不同于电极部分的材料之外，还包括左、右电极部分的一个完整的、严格的重复（electrode extension），这在使用VNL构建器件结构和ATK的整个计算过程中是必须满足也是默认满足的条件。
器件结构的一些说明
- 电流方向规定为C方向，C方向必须与AB面垂直，以确保电子不会直接从一个电极到另一个电极；
- 为节约计算量起见，对于类似此实例的在AB方向为周期性的体系，应尽量减少周期性方向的重复周期数；
- 为研究一维、二维体系，则需要在A、B方向上留足真空层；
- 中间区域可以是各种形式的材料，如金属-半导体-金属的夹层结构、分子、纳米管、纳米带等；
- 左右电极可以是不同种材料，VNL提供功能强大的异质界面建模工具，方便构建此类模型；
- 除了源漏电极之外，ATK计算中还可以增加静电门电极，用来模拟晶体管。
屏蔽区域
中间区域通常包含一些非均质的结构（界面、分子、缺陷等等），这导致中间区域的电子密度和自洽计算中的电子感受到的有效势发生变化，由此可能进一步导致中间区域的结构发生弛豫（当然，弛豫不包括电极扩展部分）。这种非均质的效应导致电势（电压）在中间区域发生降落。
完美周期结构的电子输运情况只能在0偏压下研究，这是因为在完美周期体系中由于没有散射中心导致没有电压降。对于此类体系可以使用块体材料来研究0偏压的输运情况。
NEGF理论方法假定中间区域的电势能够与电极部分平滑的过渡，也就是说在电极和中间区域之间不存在界面的散射。（注意：此处提到的界面指的是下图中蓝色方框的电极和红色方框的中间区域之间；在中间区域内部的Zn和ZnO之间的界面当然存在合理的散射，这部分可能是研究的主要对象）。为了确保这种平滑的过渡，中间区域要保留足够多的的电极材料层（含电极扩展部分），这部分称作屏蔽区域（下图中的中间区域的灰色部分）。如果屏蔽区域厚度不够，可能会导致自洽计算收敛变慢、结果不准确等问题。
上图中还显示了由于外加偏压（0.5V）导致的电势分布。左右电极的电势差为0.5eV=0.018Ha。图中明显看出，所有的电势降落都发生在ZnO部分，而金属Zn部分的电势则相当的平坦。
本征半导体材料电极中的屏蔽长度可能非常长，因此需要非常长的中间区域才能使电势平滑过渡到块体电极部分。不过在实际情况中，屏蔽距离会随掺杂浓度的提高而变短。
- 这里计算模拟的电子输运性质都是相干、弹性散射，不包含声子对电子的散射效应。但是，这种弹性散射则考虑了整个中间区域（含电极扩展部分）的特性。
- NEGF理论方法假定电极处于热平衡，这要求器件的电流比较低（电压比较低，或者电子透射概率低）。
开始计算
打开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点进行测试的脚本可以直接在这里下载。接下来只要运行两个脚本即可。
实空间网格截断能：Mesh cut-off
第一个要测试收敛性的参数是实空间网格截断能，这个参数控制实空间网格（用来投影静电势和电子密度）的精细程度。这个精细程度之所以表述为一个能量是因为网格格点的间距可以由下面公式给出： \[ \Delta x = \frac{\pi \hbar}{\sqrt{2mE}} \]
$m$是电子的质量。
脚本 mesh_loop.py 计算电极块体的基态，mesh cut-off取值从 2.5 Ha 到 80 Ha。所有控制数值精度的其他参数都保持不变。电极的费米能级（化学势）作为收敛标准的判据量。
- mesh_loop.py
# ------------------------------------------------------------- # 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。
横向的k点数
与上面的方式相同，可以继续研究电极费米能级对横向（A向和B向）的k点数的依赖关系。这里A向和B向与电流方向（C向）一定是垂直的。
在器件计算中，输运方向的k点数总是必须很大（$n_{C}=51$）。其他计算精度参数都保持不变，mesh cut-off保持为10Ha。
使用脚本kpts_loop.py可以直接进行收敛性测试，此脚本和mesh cut-off测试脚本几乎一样，只是循环变量编程横向k点数，从1到11。
- kpts_loop.py
# ------------------------------------------------------------- # 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方向是否足够长）。
首先需要对电极进行基态计算，步骤如下：
- 打开Builder，选中左电极，使用右下箭头将结构发送至Script Generator，将默认输出文件名（output file）改为
electrode.nc。
- 添加New Calculator，双击打开进行详细设置。选择Slater-Koster紧束缚近似方法计算引擎，确认“No SCF”没有勾选，将k点设为7x7x7。
- 最后在Slater-Koster basis set里选择“DFTB(CP2K,self-consistent)”。
脚本应如下：
- electrode.py
# ------------------------------------------------------------- # 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方向将成为完整周期性体系。
如果电极过短，也可以继续进行器件计算，但是很有可能导致计算无法收敛或结果不准确，这要看残余矩阵元的大小。
提示
这里的屏幕截图中主界面的插件工具的显示方式为“不显示不可用的插件”。设置的方法是打开菜单“File”–“Preference”，去掉“Show unavailable plugins”前面的勾。这会使得VNL主界面动态显示在LabFloor里选中的数据对象可用的插件，不可用的插件将不会显示。
零偏压分析
在确定了电极结构适合进行器件计算之后，就可以开始真正的输运计算。在Builder里，选中Zn-ZnO-Zn器件结构（zn-zno-zn.py），发送至Script Generator。
透射谱分析
在Script Generator里，按以下步骤设置Zn-ZnO-Zn体系的电子输运谱计算。
- 将默认输出文件名改为
zero-bias.nc - 添加New Calculator
- 添加Analysis–TransmissionSpectrum。
接下来双极添加的New Calculator，设置计算的方法：
- 选择ATK-SE：Slater-Koster(Device)紧束缚近似计算引擎。
- 将横向的k点数（A、B方向）设为5。
- 去掉No SCF Iteration前面的勾。
- 确定Slater-Koster basis set里选择的是DFTB CP2K。
下一步，编辑添加的TransmissionSpectrum：
- 将能量窗口改为从-4eV到4eV。
- 将取样点数增加为201。
- 将横向k点增加到7×7。
- 将自能（self-energy）计算方法改为Krylov。
说明
- 能量点的密度决定了计算有限偏压下电流的精度。
- 在零偏压计算时，建议选取一个很宽的窗口，以便对透射峰出现的位置有一个大概的了解。但在非零有限偏压计算时，建议测试计算窗口大小（能量点密度）的收敛性，因为能量窗口太宽（点密度太低）可能导致有些很窄的透射峰未被采样或峰型不准确。
- 在进行透射分析（或其他分析计算）时进行的k点设置不必与New Calculator里的设置一致。通常透射分析需要比自洽更多k点。
- 这里选择了速度较快的Krylov自能计算方法，但是在其他一些情况下，这种方法的近似可能会得到负的透射系数。
至此，计算脚本设置完毕，将其保存为zero-bias.py，或者在这里下载。
提交计算
使用Job Manager或命令行方法可以运行脚本。这个计算只需要几分钟。
计算结束后，要浏览一下计算作业的log文件，可以看到器件计算的整个过程。大体上计算是按以下顺序完成的：
- 左电极自洽（bulk）
- 右电极自洽（bulk）
- 器件自洽（NEGF）
- 透射谱分析。
如果log没有错误和警告信息，说明成功的进行了一个器件的电子输运计算。接下来就可以进行结果分析和可视化了。
透射分析工具
计算输出的结果会出现在LabFloor里，找到结果文件zero-bias.nc，可以看到里面包含两个数据对象：
- DeviceConfiguration包含了双电极器件结构和所使用计算方法得到的自洽电子态。
- TransmissionSpectrum包含了计算得到的透射谱的全部信息，包括计算设置和结果。
选中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点;$E=-3.16eV$处的$(k_A,k_B)=(0.15,0.15)$。
- 点击左下方的Eigenvalues。
经过短时间的计算后，所选能量、k点和自旋方向的透射本征值将会列出在白色框中。
提示
如果透射来自于几个透射通道的贡献，则都会显示在此。每个本征值都在[0,1]之间，但是本征值的总和，即k分辨透射值$T(E,k)$,可以大于1，因此总透射$T(E)$当然也一样。
接下来可以计算并显示每个本征值对应的透射本征态：
- 选择一个本征值，点击Eigenstates，过几秒钟将显示Viewer窗口。
- 选择Cut plane类型作图，点击OK。
Viewer将显示Zn-ZnO-Zn器件体系结构和中间区域的透射本征态的截面图。
接下来可以进一步调整本征态截面显示属性，以便选择更好的色调更清楚的显示透射本征态的细节：
- 在Viewer窗口的右侧，打开Properties编辑器。
- 在Cut plane标签下，将色条改为“Hot”；点击Fit，输入0.3作为本征值要显示的最大值。
截面图显示情况应该如下图所示。本征态对应了从左电极到右电极的散射态。因此，本征态在左侧总是振幅较大。较高透射本征值对应的本征态在散射区的右侧也有比较大的振幅，这反映了左侧的入射态有较大的透射概率穿过中间区域进入右电极。
在散射区左侧的透射本征态除了包括入射的部分（右行），还包括了反射态（左行）。这两部分会有一定的迭加干涉效应。从上图中可以看到左侧第一个Zn原子上的本征态波幅较小，这是由于入射和反射的态的相消干涉造成的。在ZnO的右侧，本征态只有透射（右行）部分，因此观察不到类似的干涉效应。
你可以再尝试画出费米能级（能量$E = 0 eV$）、$k_A=k_B=0$处的透射本征态。调整作图属性之后可以得到如下图所示。这时的本征态的波幅在散射区的中部到右部几乎消失为零，对应于很小的透射值，也就是说只有很小的散射本征态透过，而大部分都被Zn-ZnO界面反射。
提示
如果采用等值面图的方式对透射本征态作图，等值面对应了波函数的波幅，颜色则对应了波函数的相位。