这是本文档旧的修订版！
使用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自能计算方法，但是在其他一些情况下，这种方法的近似可能会得到负的透射系数。