使用QuantumATK研究电子输运

本文用简单的实例介绍使用QuantumATK研究电子输运的基本概念，计算设置应该注意的事项和对体系输运性质进行分析的常用方法。

为了快速完成计算，这里用到的体系是Zn-ZnO-Zn，使用ATK-SE半经验方法进行计算。体系的结构见下图：

要了解如何构建可靠的器件模型，请参考：器件体系的建模与计算。

提示

本教程使用特定版本的QuantumATK创建，因此涉及的截图和脚本参数可能与您实际使用的版本略有区别，请在学习时务必注意。

实例简述

实例计算中使用紧束缚近似（TB）的半经验方法，结合非平衡态格林函数方法（NEGF）研究上述器件结构的电子输运特性。理论方法的详细描述见：参考手册。

这里主要介绍器件电子输运计算的核心概念与在QuantumATK中可以进行的常用计算：

器件的基本结构
选择合适的计算精度
检查块体材料是否适合用作电极
计算0偏压和有限偏压下的透射谱
使用透射分析工具分析透射谱
计算透射本征态
计算投影器件态密度
计算器件内电势分布（电压降）
计算、分析IV曲线

器件体系的几何结构

与通常的计算模拟中常见的块体材料的周期性体系（bulk）和分子体系的孤立结构（molecule）不同，用于电子输运计算的器件体系（device）比较特殊。器件结构由左、右两个电极（electrode，分别是周期性体系）和中间区域（central region，分子体系）构成。在晶体管器件中，左、右电极也称为源、漏电极。中间区域则是一个有限的区域，也称为散射区（scattering region），由于载流子从一个电极到另一个电极传输时会受到中间区域的散射，因此，中间区域是决定器件特性的关键区域。中间区域除了可能存在的不同于电极部分的材料之外，还包括左、右电极部分的一个完整的、严格的重复（electrode extension），这在使用VNL构建器件结构和QuantumATK的整个计算过程中是必须满足也是默认满足的条件。

器件结构的一些说明

电流方向规定为C方向，C方向必须与AB面垂直，以确保电子不会直接从一个电极到另一个电极；
为节约计算量起见，对于类似此实例的在AB方向为周期性的体系，应尽量减少周期性方向的重复周期数；
为研究一维、二维体系，则需要在A、B方向上留足真空层；
中间区域可以是各种形式的材料，如金属-半导体-金属的夹层结构、分子、纳米管、纳米带等；
左右电极可以是不同种材料，VNL提供功能强大的异质界面建模工具，方便构建此类模型；
除了源漏电极之外，QuantumATK计算中还可以增加静电门电极，用来模拟晶体管。

屏蔽区域

中间区域通常包含一些非均质的结构（界面、分子、缺陷等等），这导致中间区域的电子密度和自洽计算中的电子感受到的有效势发生变化，由此可能进一步导致中间区域的结构发生弛豫（当然，弛豫不包括电极扩展部分）。这种非均质的效应导致电势（电压）在中间区域发生降落。

完美周期结构的电子输运情况只能在0偏压下研究，这是因为在完美周期体系中由于没有散射中心导致没有电压降。对于此类体系可以使用块体材料来研究0偏压的输运情况。

NEGF理论方法假定中间区域的电势能够与电极部分平滑的过渡，也就是说在电极和中间区域之间不存在界面的散射。（注意：此处提到的界面指的是下图中蓝色方框的电极和红色方框的中间区域之间；在中间区域内部的Zn和ZnO之间的界面当然存在合理的散射，这部分可能是研究的主要对象）。为了确保这种平滑的过渡，中间区域要保留足够多的的电极材料层（含电极扩展部分），这部分称作屏蔽区域（下图中的中间区域的灰色部分）。如果屏蔽区域厚度不够，可能会导致自洽计算收敛变慢、结果不准确等问题。

上图中还显示了由于外加偏压（0.5V）导致的电势分布。左右电极的电势差为0.5eV=0.018Ha。图中明显看出，所有的电势降落都发生在ZnO部分，而金属Zn部分的电势则相当的平坦。

本征半导体材料电极中的屏蔽长度可能非常长，因此需要非常长的中间区域才能使电势平滑过渡到块体电极部分。不过在实际情况中，屏蔽距离会随掺杂浓度的提高而变短。

这里计算模拟的电子输运性质都是相干、弹性散射，不包含声子对电子的散射效应。但是，这种弹性散射则考虑了整个中间区域（含电极扩展部分）的特性。
NEGF理论方法假定电极处于热平衡，这要求器件的电流比较低（电压比较低，或者电子透射概率低）。

开始计算

打开VNL，创建一个新的Project。一般情况下，可以使用Builder构建新的结构，这里可以直接下载结构文件，放在打开的Project文件夹下。将脚本拖动到Builder图标上即可打开Builder窗口看到结构。

设置电极参数的精度

QuantumATK的器件计算通常分为多个步骤。首先是将电极和中间区域作为周期性的块体材料体系进行自洽电子密度计算。随后开始的器件计算首先将中间区域的电子密度作为初始猜测，使两个电极和中间区域在接触的边界匹配。最后计算得到中间区域的自洽电子结构。

上述步骤中，电极块体区域的收敛性很重要，既要保证电极长度足够，也要保证数值精度足够。这一节介绍两个重要的数值精度控制参数：实空间密度网格截断能（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界面反射。

提示

如果采用等值面图的方式对透射本征态作图，等值面对应了波函数的波幅，颜色则对应了波函数的相位。

透射谱计算的对k点的收敛性测试

如前所述，通常情况下要得到收敛很好的透射谱需要比SCF自洽更多的k点。例如，下图显示了三个能量下的总透射（所有k点平均）对横向（A、B方向）的k点数的关系图。得到下图的脚本也附在后面（transmission_kpts.py）。

transmission_kpts.py

# Restart from the previous calculation.
device_configuration = nlread('zero-bias.nc', object_id='gID000')[0]
 
# Array with k-values.
numbers_of_kpoints = range(3,40,2)
 
# Energies at which to calculate the transmission function.
energies = [-3,0,3]
 
# Array with transmission values.
T = numpy.zeros([len(energies), len(numbers_of_kpoints)])
 
# Loop over numbers of k-points.
for i, num_k in enumerate(numbers_of_kpoints):        
    # -------------------------------------------------------------
    # Transmission spectrum
    # -------------------------------------------------------------
    transmission_spectrum = TransmissionSpectrum(
        configuration=device_configuration,
        energies=numpy.array(energies)*eV,
        kpoints=MonkhorstPackGrid(num_k,num_k),
        energy_zero_parameter=AverageFermiLevel,
        infinitesimal=1e-06*eV,
        self_energy_calculator=KrylovSelfEnergy(),
        )
 
    # Put the k-point averaged transmission into the array.
    T[:,i] = transmission_spectrum.evaluate()
 
# Convert to numpy array.
numbers_of_kpoints  = numpy.array(numbers_of_kpoints)
 
# Plot the results.
import pylab
 
pylab.subplot(311)
pylab.plot(numbers_of_kpoints , T[0,:] ,'ro-',label='E=-3 eV')
pylab.xlabel('# k-points',fontsize=16)
pylab.ylabel('Transmission',fontsize=16)
pylab.legend()
 
pylab.subplot(312)
pylab.plot(numbers_of_kpoints , T[1,:] ,'ro-',label='E=0 eV')
pylab.xlabel('# k-points',fontsize=16)
pylab.ylabel('Transmission',fontsize=16)
pylab.legend()
 
pylab.subplot(313)
pylab.plot(numbers_of_kpoints , T[2,:] ,'ro-',label='E=3 eV')
pylab.xlabel('# k-points',fontsize=16)
pylab.ylabel('Transmission',fontsize=16)
pylab.legend()
 
# Show the plot.
pylab.show()

图中可以看出，不同能量处的透射收敛情况非常不同。费米能量处大约 5~7 个k点即收敛，但在 -3eV 和 3eV 处，看上去要取13个k点才获得收敛的结果。

回顾透射谱的形状，可以注意到在透射谱的 -3 eV 处有许多尖峰。通常存在尖峰的透射谱需要更多的k点才能得到收敛的结果。然而，有时候可能由于尝试计算时k点过于粗糙导致某些透射谱上某些孤立尖峰被忽略，为了避免这种情况，就需要仔细研究k分辨透射图，以大致确定到底需要多少的k点取样。

此外，区域里存在大量（尖锐的）透射峰时，也需要增加能量取样点数，以改善为了计算电流时对能量空间进行积分的准确性。

进一步分析

这一节介绍几种进一步分析零偏压透射的方法。ATK 的一个重要特性是，不用从一开始就定义所有要进行的分析计算；你可以随时进行各种分析计算得到有用的物理量，而且不必从头进行费时的自洽计算。

接下来，我们计算Zn–ZnO–Zn体系的器件态密度（DDOS）和差别静电势（EDP）：

打开 Script Generator 窗口；
将输出文件改为analysis.nc，添加以下几个步骤：
- Analysis from File；
- DeviceDensityOfStates；
- ElectrostaticDifferencePotential。

打开“Analysis from File”。这个功能决定QuantumATK脚本从何处读取已经自洽的零偏压计算得到的态。你需要指定保存了DeviceConfiguration的文件名，以在此基础上继续计算。
选择zero-bias.nc，数据对象ID “gID000”。从2016版本开始，可以点击Detect获得可用的数据对象ID。点击OK。

提示

每个数据对象的ID号显示在LabFloor里对应的图标下方。为区别同一个文件里的不同Configuration对象，你可以点击右侧的General information插件查看最重要的计算参数等信息。

接下来，编辑DeviceDensityOfStates的设置，与计算透射谱的设置类似：
- 能量从 -4 eV 到 4 eV，共201个点；
- k点7×7；
- Krylov自能计算器
不需要修改 ElectrostaticDifferencePotential。

分析脚本就准备好了。保存为analysis.py，并用Job Manager或命令行执行，这个计算只需要几分钟。完整的脚本如下：

analysis.py

# -------------------------------------------------------------
# Analysis from File
# -------------------------------------------------------------
configuration = nlread(u'zero-bias.nc', object_id='gID000')[0]
 
# -------------------------------------------------------------
# Device Density Of States
# -------------------------------------------------------------
device_density_of_states = DeviceDensityOfStates(
    configuration=configuration,
    energies=numpy.linspace(-4,4,201)*eV,
    kpoints=MonkhorstPackGrid(7,7),
    contributions=All,
    energy_zero_parameter=AverageFermiLevel,
    infinitesimal=1e-06*eV,
    self_energy_calculator=KrylovSelfEnergy(),
    )
nlsave('analysis.nc', device_density_of_states)
nlprint(device_density_of_states)
 
# -------------------------------------------------------------
# Electrostatic Difference Potential
# -------------------------------------------------------------
electrostatic_difference_potential = ElectrostaticDifferencePotential(configuration)
nlsave('analysis.nc', electrostatic_difference_potential)

计算结束后，会LabFloor里会出现两个新的结果：

DeviceDensityOfStates

ElectrostaticDifferencePotential

器件态密度

选中DeviceDensityOfStates对象，点击右侧的 2D Plot 工具可以作图可视化。

打开的 2D Plot 窗口左侧显示了DDOS图，右侧则用以选择投影DOS的原子序号和轨道。

选择zero-bias.nc即可看到器件体系结构，选择8、9、10、11号原子用于DOS投影；
放大DDOS如下图，显示ZnO中间部分的投影器件态密度。

注意到，DDOS的形状与透射谱形状十分类似：两者在费米能级附近都有相当宽的区域接近零，表明了ZnO的半导体特性；能量低于 -2 eV 处有很多尖峰；还有，两者在 +3 eV处都有一个宽峰。透射谱与DDOS的对应关系在器件计算中是十分常见的特性。

差别静电势

差别静电势（ElectrostaticDifferencePotential）数据对象包含了一组三维网格数据，代表了自洽得到的价电子密度的静电势（自洽价电子密度代入泊松方程解得的电势）与原子价电子密度直接重叠得到的静电势的差别。后面会用到这个物理量来计算电压降。

使用 Viewer 工具可以可视化这种三维网格数据，效果与上面显示透射本征态十分类似。此外，还有一种分析三维网格数据的方法是使用一维投影工具 1D Projector。这个工具可以将电势或电子密度在某一个方向作图，而在另外两个垂直方向上进行积分或平均。对于器件体系，最有意义的是沿输运方向（C方向）进行电势分析，而对A和B方向进行平均。

选中 LabFloor 中的 ElectrostaticDifferencePotential 数据对象，打开右侧的1D Projector插件。
将Projection Type改为 Average，点击 Add Line。

你现在可以看到如上图的一条曲线，代表了差别静电势$\Delta V_H$ 随 C 方向坐标（0~1）的变化关系。

图中可以看到由于Zn–ZnO界面上显著的电荷转移导致了自洽的静电势相对于原子价电子电势发生了显著的变化。在金属Zn区域，振荡相当小，表明这个区域的的自洽电子密度分布与孤立原子价电子密度分布十分接近。

有限偏压计算

以上所有计算都仅限于零偏压的情况（左右电极上电压都为零）。接下来这一节介绍如何使用零偏压的结果作为起点进行非零偏压下的计算，此时左右电极电压将不再相同。这将导致左右电极费米能级的相对移动，这也是电子器件的运作的基本模式。

从零偏压结果出发逐步增加电压进行非零偏压计算是快速、准确的获得收敛结果的最佳技术方案。这里对器件进行0.5V偏压的计算。

计算设置

在VNL 项目文件列表里，勾选zero-bias.nc文件，在LabFloor里查看并将DeviceConfiugration（gID000）拖动到Script Generator上。这会直接设置好一个Slater-Koster的计算，所选参数也与之前的零偏压计算完全相同。

添加三个计算步骤：Initial State，Analysis → ElectrostaticDifferencePotential，以及 Analysis → TransmissionSpectrum.

将输出文件改为finite-bias.nc，打开New Calculator来设置电压：

将左电极电压设为 -0.5V。

接下来，打开“Initial State”告诉QuantumATK使用零偏压的计算得到的态作为有限偏压的初始猜测：

选中 Use old calculation，选择zero-bias.nc文件作为Old Calculation 的文件名。

最后，设置透射谱计算，与零偏压完全相同：

能量范围 -4eV 到 4eV，共201个点；
7×7 横向k点；
Krylov自能计算器。

将脚本保存为finite-bias.py并运行。这个计算比零偏压计算耗时略长。完整的脚本如下。

finite-bias.py

# -------------------------------------------------------------
# Two-probe Configuration
# -------------------------------------------------------------
 
# -------------------------------------------------------------
# Left Electrode
# -------------------------------------------------------------
 
# 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
left_electrode_lattice = UnitCell(vector_a, vector_b, vector_c)
 
# Define elements
left_electrode_elements = [Zinc, Zinc]
 
# Define coordinates
left_electrode_coordinates = [[ 1.3325    ,  0.76931925,  1.23675   ],
                              [ 1.3325    , -0.76931925,  3.71025   ]]*Angstrom
 
# Set up configuration
left_electrode = BulkConfiguration(
    bravais_lattice=left_electrode_lattice,
    elements=left_electrode_elements,
    cartesian_coordinates=left_electrode_coordinates
    )
 
# -------------------------------------------------------------
# Right Electrode
# -------------------------------------------------------------
 
# 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
right_electrode_lattice = UnitCell(vector_a, vector_b, vector_c)
 
# Define elements
right_electrode_elements = [Zinc, Zinc]
 
# Define coordinates
right_electrode_coordinates = [[ 1.0280478 ,  0.59354366,  1.23675   ],
                               [ 1.0280478 , -0.94509484,  3.71025   ]]*Angstrom
 
# Set up configuration
right_electrode = BulkConfiguration(
    bravais_lattice=right_electrode_lattice,
    elements=right_electrode_elements,
    cartesian_coordinates=right_electrode_coordinates
    )
 
# -------------------------------------------------------------
# Central Region
# -------------------------------------------------------------
 
# 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, 40.5741587317]*Angstrom
central_region_lattice = UnitCell(vector_a, vector_b, vector_c)
 
# Define elements
central_region_elements = [Zinc, Zinc, Zinc, Zinc, Zinc, Zinc, Oxygen, Zinc, Oxygen, Zinc,
                           Oxygen, Zinc, Oxygen, Zinc, Oxygen, Zinc, Zinc, Zinc, Zinc, Zinc,
                           Zinc]
 
# Define coordinates
central_region_coordinates = [[  1.3325    ,   0.76931925,   1.23675   ],
                              [  1.3325    ,  -0.76931925,   3.71025   ],
                              [  1.3325    ,   0.76931925,   6.18375   ],
                              [  1.3325    ,  -0.76931925,   8.65725   ],
                              [  1.3325    ,   0.76931925,  11.13075   ],
                              [  1.3325    ,  -0.76931925,  13.60425   ],
                              [  1.3325    ,  -0.76931925,  15.60909881],
                              [  2.3605478 ,   0.17577554,  16.27738175],
                              [  2.3605478 ,   0.17577554,  18.28223056],
                              [  1.0280478 ,   0.59354366,  18.95051349],
                              [  1.0280478 ,   0.59354366,  20.9553623 ],
                              [  1.3325    ,  -0.76931925,  21.62364524],
                              [  1.3325    ,  -0.76931925,  23.62849405],
                              [  2.3605478 ,   0.17577554,  24.29677699],
                              [  2.3605478 ,   0.17577554,  26.3016258 ],
                              [  1.0280478 ,   0.59354366,  26.96990873],
                              [  1.0280478 ,  -0.94509484,  29.44340873],
                              [  1.0280478 ,   0.59354366,  31.91690873],
                              [  1.0280478 ,  -0.94509484,  34.39040873],
                              [  1.0280478 ,   0.59354366,  36.86390873],
                              [  1.0280478 ,  -0.94509484,  39.33740873]]*Angstrom
 
# Set up configuration
central_region = BulkConfiguration(
    bravais_lattice=central_region_lattice,
    elements=central_region_elements,
    cartesian_coordinates=central_region_coordinates
    )
 
device_configuration = DeviceConfiguration(
    central_region,
    [left_electrode, right_electrode]
    )
 
# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = DFTBDirectory("cp2k/scc/")
 
#----------------------------------------
# Pair Potentials
#----------------------------------------
pair_potentials = DFTBDirectory("cp2k/scc/")
 
#----------------------------------------
# Numerical Accuracy Settings
#----------------------------------------
left_electrode_k_point_sampling = MonkhorstPackGrid(
    na=5,
    nb=5,
    nc=100,
    )
left_electrode_numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=left_electrode_k_point_sampling,
    )
 
right_electrode_k_point_sampling = MonkhorstPackGrid(
    na=5,
    nb=5,
    nc=100,
    )
right_electrode_numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=right_electrode_k_point_sampling,
    )
 
device_k_point_sampling = MonkhorstPackGrid(
    na=5,
    nb=5,
    nc=100,
    )
device_numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=device_k_point_sampling,
    )
 
#----------------------------------------
# Iteration Control Settings
#----------------------------------------
left_electrode_iteration_control_parameters = IterationControlParameters(
    tolerance=4e-05,
    )
 
right_electrode_iteration_control_parameters = IterationControlParameters(
    tolerance=4e-05,
    )
 
device_iteration_control_parameters = IterationControlParameters(
    tolerance=4e-05,
    )
 
#----------------------------------------
# Poisson Solver Settings
#----------------------------------------
left_electrode_poisson_solver = FastFourier2DSolver(
    boundary_conditions=[[PeriodicBoundaryCondition,PeriodicBoundaryCondition],
                         [PeriodicBoundaryCondition,PeriodicBoundaryCondition],
                         [PeriodicBoundaryCondition,PeriodicBoundaryCondition]]
    )
 
right_electrode_poisson_solver = FastFourier2DSolver(
    boundary_conditions=[[PeriodicBoundaryCondition,PeriodicBoundaryCondition],
                         [PeriodicBoundaryCondition,PeriodicBoundaryCondition],
                         [PeriodicBoundaryCondition,PeriodicBoundaryCondition]]
    )
 
#----------------------------------------
# Electrode Calculators
#----------------------------------------
left_electrode_calculator = SlaterKosterCalculator(
    basis_set=basis_set,
    pair_potentials=pair_potentials,
    numerical_accuracy_parameters=left_electrode_numerical_accuracy_parameters,
    iteration_control_parameters=left_electrode_iteration_control_parameters,
    poisson_solver=left_electrode_poisson_solver,
    )
 
right_electrode_calculator = SlaterKosterCalculator(
    basis_set=basis_set,
    pair_potentials=pair_potentials,
    numerical_accuracy_parameters=right_electrode_numerical_accuracy_parameters,
    iteration_control_parameters=right_electrode_iteration_control_parameters,
    poisson_solver=right_electrode_poisson_solver,
    )
 
#----------------------------------------
# Device Calculator
#----------------------------------------
calculator = DeviceSlaterKosterCalculator(
    basis_set=basis_set,
    pair_potentials=pair_potentials,
    numerical_accuracy_parameters=device_numerical_accuracy_parameters,
    iteration_control_parameters=device_iteration_control_parameters,
    electrode_calculators=
        [left_electrode_calculator, right_electrode_calculator],
    electrode_voltages=( -0.5*Volt, 0.0*Volt)
    )
 
device_configuration.setCalculator(calculator)
 
# -------------------------------------------------------------
# Initial State
# -------------------------------------------------------------
old_calculation = nlread('zero-bias.nc', DeviceConfiguration)[0]
device_configuration.setCalculator(
    calculator,
    initial_state=old_calculation,
)
device_configuration.update()
nlsave('finite-bias.nc', device_configuration)
nlprint(device_configuration)
 
# -------------------------------------------------------------
# Electrostatic Difference Potential
# -------------------------------------------------------------
electrostatic_difference_potential = ElectrostaticDifferencePotential(device_configuration)
nlsave('finite-bias.nc', electrostatic_difference_potential)
 
# -------------------------------------------------------------
# Transmission Spectrum
# -------------------------------------------------------------
transmission_spectrum = TransmissionSpectrum(
    configuration=device_configuration,
    energies=numpy.linspace(-4,4,201)*eV,
    kpoints=MonkhorstPackGrid(7,7),
    energy_zero_parameter=AverageFermiLevel,
    infinitesimal=1e-06*eV,
    self_energy_calculator=KrylovSelfEnergy(),
    )
nlsave('finite-bias.nc', transmission_spectrum)
nlprint(transmission_spectrum)

关于电极电压和电流

对于不含门电极的体系，只有左右电极的电压差会影响计算结果。也就是说，电压设为 -0.25/0.25V 和 -0.5/0.0V设置完全一样。

偏压的定义是左右电压的差别（右电压减去左电压）。正电流对应于电荷从左移动到右，也就是电子从右移动到侧（电子电荷为负）。因此，当左电极电压高于右电极时，电流为正；而对于上面的计算的偏压设置，将得到负电流。

NEGF 环路积分

对于非零偏压计算，一个影响收敛的重要的参数设置是“Contour Integral Parameters”里的 Real axis point density。这个参数决定了在进行左、右谱密度矩阵沿能量实轴方向上在左、右电极化学势区间里求算积分时使用的能量点的密度。对于当前计算，默认值即可；但更一般的计算，需要测试这个参数对收敛性的影响。

另外两个重要的参数是 Integral lower bound 和 Circle points。NEGF环路积分需要包含所有的能量轴上的本征值，但是本征值分布的区间大小无法预先得知。默认的积分限是费米面以下1.5 Ha，但对于使用“较硬”的赝势（大量的价电子、半芯态）进行的计算，很可能在这个能量以下还有本征态。这时需要增加 Integral lower bound，并同时增加 Circle points。器件收敛困难多数是由这个问题造成。

结果分析

回到VNL主窗口，打开finite-bias.nc文件，同样使用 Transmission Analyzer 分析透射谱。很容易注意到左右电极的化学势分开了0.5eV，对称的位于平均费米能级（参考能量零点）上下。还可以看到体系中有非常小的负电流，这与我们增加的负偏压一致。

非零偏压得到透射谱与零偏压基本一致。回顾之前关于DDOS的分析也可以发现，透射函数的形状基本上是由ZnO得态密度决定，与电极的电子态几乎没有什么关系。由于偏压窗口中的透射非常小，体系电流也就很小，ZnO的DDOS变化也很小。更一般的情况并非如此，尤其是如果偏压大到能够跨越散射区域的带隙的话。

电压降

接下来的任务是分析器件输运方向上的电势分布，或者叫电压降（由于外加偏压导致的电势重新分布）。这可以由非零偏压和零偏压下的差别静电势作差得到： $$ \Delta V_H(\vec{r})=\delta V_H^{\mathtt{finite-bias}}-\delta V_H^{\mathtt{zero-bias}} $$

差别静电势是自洽得到的价电子密度对应的电势与孤立原子价电子密度叠加对应的电势的差别。由于原子价电子密度在零偏压和非零偏压下一样，因此上述得到的$\Delta V_H(\vec{r})$就是由于外加电压导致的电势重新分布，即电压降。

VNL中的网格数据运算（Gird Operation）工具在这里非常方便。这个工具可以将两套网格数据相加或作差，甚至乘法或除法。下面用这个工具计算$\Delta V_H(\vec{r})$，并保存文件，之后用Viewer工具可视化。

LabFloor里默认按文件名分组。为方便起见，可以改为Group by Item Type。
之后选择两个ElectrostaticDifferencePotential数据对象，打开Grid Operation工具。

在这个工具窗口中，每套数据都给定了一个别名，根据选择数据时的先后顺序，有限偏压和零偏压的数据分别被命名为A1和B1，或者相反。为计算电压降（有限偏压的数据减去零偏压数据），在上面输入B1- A1 （或者相反），回车。

你马上就可以在3D 显示窗口里画出这个电压降数据的切面图或轮廓图，或者将结果保存到一个新的NetCDF文件，点击Save即可。

将网格数据差别保存为voltage_drop.nc，在LabFloor里可以看到出现了一个新的GridValue对象。

$f(r)$ GridValue 是表示3D网格数据的通用符号。

选择DeviceConfiguration，用Viewer打开这个器件结构。然后将GridValue对象拖动到Viewer上，即可在器件结构上显示电压降。选择 Cut Plane模式，并将色条改为RGB。

显示结果如下图所示。电压降明显的出现在中间区域的ZnO附近，表明电子散射发生的主要区域。

为了定量分析，打开 1D Projector，选择投影类型为Average，点击Add Line。应该可以看到如下图的显示的曲线。电势降落的单位是Ha。因此总电势降落为0.0183 Ha * 27.2 eV/Ha = 0.5 eV，即外加偏压。

也应注意到接近电极部分的电势是平坦的，说明中间区域内两端的Zn层提供了足够的屏蔽。

IV曲线

本教程的最后一部分器件非零偏压分析是IV曲线计算。这可以通过使用IVCurve分析计算工具得到，并用IV-Plot可视化工具分析。

计算设置

打开一个新的Script Generator窗口，按下面操作：

添加“Analysis from File”；
添加Analysis→IVCurve;
将输出文件名改为finite-bias.nc（把IV曲线数据追加写入现有文件）。

编辑 Analysis from File：

将NetCDF文件设置为zero-bias.nc，使用gID000作为计算的起点。

编辑IVCurve:

设置偏压结束点V₁=1V;
设置偏压点数为5；
设置横向k点：5×5；
选择Krylov自能计算器。

计算设置到此完成。保存脚本为iv-curve.py，用Job Manager或命令行运行。这个计算需要一些时间，因为在每个偏压点都需要进行自洽计算和透射谱分析（大约10-20分钟）。完整的计算脚本如下。计算结果可以直接下载：finite-bias.nc。

iv-curve.py

# -------------------------------------------------------------
# Analysis from File
# -------------------------------------------------------------
configuration = nlread(u'zero-bias.nc', object_id='gID000')[0]
 
# -------------------------------------------------------------
# IV Curve
# -------------------------------------------------------------
biases = [0.000000, 0.250000, 0.500000, 0.750000, 1.000000]*Volt
 
iv_curve = IVCurve(
    configuration=configuration,
    biases=biases,
    energies=numpy.linspace(-2,2,101)*eV,
    kpoints=MonkhorstPackGrid(5,5),
    self_energy_calculator=KrylovSelfEnergy(),
    energy_zero_parameter=AverageFermiLevel,
    infinitesimal=1e-06*eV,
    selfconsistent_configurations_filename_prefix="ivcurve_selfconsistent_configuration_",
    log_filename_prefix="ivcurve_"
    )
nlsave('finite-bias.nc', iv_curve)
nlprint(iv_curve)

可视化

计算完成后，IV曲线数据对象出现在finite-bias.nc文件下。此外还出现了五个新的输出文件，文件名开头为ivcurve。这些文件里保存了每个偏压点下的SCF自洽结果。

IVCurve数据对象包含了计算得到IV曲线以及每个偏压下的透射谱。

选择IVCurve，点击右侧的 IV-Plot 分析工具，IV曲线分析窗口出现。点击additional plots，可以显示微分电导（dI/dV），计算得到的透射谱以及谱电流。

每个透射谱和谱电流对应了一个偏压点。偏压窗口对应的能量积分区间在曲线上用蓝色显示。将鼠标沿着IV曲线或dI/dV上的点移动，对应的透射谱和谱电流将被高亮显示。

需特别注意谱电流：为得到准确度的电流，所有可能对总电流有贡献的部分（偏压窗口）都必须包含在能量范围以内。如果不是，应该用更大的能量范围重新进行计算。

使用更精细的偏压点可以得到更精确的微分电导（dI/dV曲线）。为此，你可以将偏压点数增加（例如21点），重新进行IVCurve计算分析，或在此下载结果文件：iv-curve_fine.nc

更精细的IV曲线如下图所示。与上面的粗略采样相比，可以发现IV曲线大致形状一致，但dI/dV曲线却大不相同。

总结

上面介绍了如何用QuantumATK进行电子输运的计算。这篇教程要达到的主要目的有：

了解了器件基本的几何结构，可以继续学习如何构建其他的器件结构：VNL基本操作；
了解最重要的收敛性参数，例如横向k点，密度网格截断，NEGF环路积分参数设置等；
学会进行零偏压和非零偏压的透射谱计算，并使用Transmission Analyzer分析计算结果；
掌握如何进行分析计算，如IV曲线、投影器件态密度、差别静电势、透射本征值与本征态；
学会如何用3D Viewer可视化三维数据对象，并作一维投影。

参考

英文教程原文：英文

目录

使用QuantumATK研究电子输运

提示

实例简述

器件体系的几何结构

器件结构的一些说明

屏蔽区域

开始计算

设置电极参数的精度

将器件结构分割开

实空间网格截断能：Mesh cut-off

横向的k点数

检查电极的几何结构和大小

说明

提示

零偏压分析

透射谱分析

说明

提交计算

透射分析工具

透射本征态

提示

提示

透射谱计算的对k点的收敛性测试

进一步分析

提示

器件态密度

差别静电势

有限偏压计算

计算设置

关于电极电压和电流

NEGF 环路积分

结果分析

电压降

IV曲线

计算设置

可视化

总结

参考