本教程专门介绍分子器件的计算和分析，研究的经典器件由一个对苯二硫酚（DTB）与两个金电极表面构成，这个结构在很多文献里都有报道。

本教程中的计算时间很长，作为入门的ATK量子输运计算教程，请参考：

• 使用QuantumATK研究电子输运

本教程中使用的分子导电结的建模请参见：

• 构建金-DTB-金分子结

本教程中使用的ATK-DFT方法参见ATK的使用手册。

参考这里构建结构模型，或直接下载脚本builder_device.zip，并将脚本拖动到Builder中，如下图。

选中器件结构，点击 Send to 将结构发送至 Script Generator。

• 在 Script Generator 中添加三个计算单元：
  • New Calculator
  • Analysis – DeviceDensityOfStates
  • Analysis – TransmissionSpectrum
• 将默认输出文件名改为 au-dtb-au.nc。
• 双击打开 New Calculator 单元，进行如下设置：
  • 设置 k 点为3×3×50；
  • 在 Basis set/exchange correlation 里，将金元素的基组设置为 SingleZetaPolarized 以加快计算速度。

• 双击打开 DeviceDensityOfStates 单元，进行如下设置：
  • 设置 k 点为 3×3；
  • 设置能量范围为 (-3,3) eV
  • 将能量点数从 101 改为 401
• 打开 TransmissionSpectrum 单元，进行相同的设置：
  • 设置 k 点为 3×3；
  • 设置能量范围为 (-3,3) eV
  • 将能量点数从 101 改为 401
• 最后，将脚本保存为 au_dtb_au.py。这个脚本可以直接在这里下载：au-dtb-au.zip

点击 Script Generator 右下角的 Send To 将脚本传送至 Job Manager。在 Job Manager 里点击 Run Queue 启动计算。此计算在单机上大约耗时一小时。

这一节主要介绍如何分析并提取零偏压下的数据：

• 将透射谱与苯环的分波态密度（PDDOS）进行对比；
• 计算分子投影自洽哈密顿量（MPSH），并用MPSH态理解透射的机制；
• 计算透射本征态；
• 计算能量分解的局域态密度（LDOS），得到能级的空间分布。

回到 VNL 主窗口，应该可以在 Project Files 里看到 au-dtb-au.nc 文件，并在 LabFloor 里看到文件内容：

• DeviceConfiguration 包含双电极结构和使用设定的计算方法自洽得到的电子态；
• TransmissionSpectrum 包含了透射谱计算的全部信息（计算设置和结果）；
• DeviceDensityOfState 包含了所有器件区域态密度的信息。

选中透射谱数据对象，点击右侧 Transmission Analyzer 打开透射谱分析工具。左侧题为Spectrum的部分显示的是综合k点的透射强度，为能量的函数。图中的虚线为能量零点（左右点击费米能的平均值）。

图中用红色圆点标出的能量和透射概率值在图的下方给出。你可以用鼠标在曲线上移动圆点，读取想要的数据。Curves下拉菜单可以选择不同的自旋分量。要指出的是只有对于非共线自旋计算，自旋的X、Y、Z分量才有意义。

在任意的能量处（红色圆点标记），总透射是通过在k空间取样并计算透射系数得到的。右侧 Coefficients 图即画出了透射系数与倒空间矢量 $k_A$ 和 $k_B$ 的关系。例如，在能量为 -1.4 eV 的透射峰处，$\Gamma$ 点贡献却最小，如下图。

在 -3、-1.4、2.5 eV处有三个宽峰；2.4 eV附近还有一个窄峰，T(E) 超过了 1。后面我们会研究这些峰的来源。

接下来选中 DeviceDensityOfStates 数据对象，点击 2D Plot…按钮，打开PDDOS分析工具窗口。要将所有PDDOS线画在一起，点击左上的 Option 按钮，关闭 flipping。

选中所有的碳和氢原子（鼠标左键拖动方框选中，或者按住Ctrl键用鼠标单击逐一选中）

MPSH可以帮助理解在有电极存在时的分子能级的情况。MPSH是通过对全自洽的哈密顿量的分子部分进行对角化得到的。

打开 Script Generator 工具：

• 双击添加 Analysis from File 单元；
• 选择 au-dtb-au.nc，选择 DeviceConfiguration 数据对象（gID000）；
• 添加 Analysis – MolecularEnergySpectrum；
• 将输出文件名改为 dtb_analysis.nc；
• 打开 MolecularEnergySpectrum 单元，输入所有碳和氢原子的序号（37,38,39,40,41,42,43,44,45,46）；

• 将脚本发送至 Job Manager 并执行计算（几秒钟即可完成）。MPSH能量将出现在log文件里。

要计算MPSH某个能级的本征态，可以在Windows系统的任务栏或者 VNL 的主窗口 “Windows” 菜单里重新打开 Script Generator：

• 将输出文件名改为 molecular_eigenstate.nc；
• 添加四个 Analysis – Eigenstate单元；
• 每个单元都设置好碳和氢原子的序号，四个单元的Quantum Number 分别设置为 13、14、15、16。

• 将 MolecularEnergySpectrum 删去（已经计算过了）；
• 将脚本提交至 Job Manager 并执行计算。

计算结束后，文件 molecular_eigenstate.nc 应该已经默认勾选，其中的数据对象则可以在 LabFloor 里看到。选择第一个本征态对象，点击 Viewer按钮，选择 Isosurface。在 Viewer 里，打开 Property 菜单，选择 Isosurface，进行如下设置：

• 设置 Isovalue 0.16；
• 设置 Color map 为 BlueRed；
• 勾选 $\pm$Isovalue

要同时将结构显示在图上，只需要将 DeviceConfiguration 数据对象拖动到上面已经显示了本征态的Viewer窗口里。

注意到，所有的态都是关于苯分子平面反对称的，这说明这几个态都是苯分子的 $\pi$-轨道，由与平面垂直的碳原子 6 个 pz 轨道线性组合而成的，因此这样的 $\pi$-轨道共有 6 个。

观察不同本征态的对称性可以确定还有两个 $\pi$-轨道分别是 10 和 18。用 Editor 将之前的脚本打开， 删去两个 Eigenstate 单元，将剩下的两个计算单元的 quantum number 改为 10 和18。计算后即可得到这两个轨道，如下图（isovalue=0.16）。

接下来可以研究透射谱上峰位置的透射本征态，即-3、-1.5、2.4 V处。

• 重新打开 Script Generator，删去 Eigenstate 单元；
• 添加三个 Analysis → TransmissionEigenvalues 单元；
• 分别将三个 TransmissionEigenvalues 的能量设置成 -3、-1.5、2.4；
• 将脚本发送到 Job Manager 并执行计算，查看log文件。

透射本征值是通过对透射矩阵进行对角化得到的。透射本征值的个数表明了通过分子的输运通道的个数，数值的大小则是该通道的透射强度。这些透射本征值是真正的透射概率，因此应位于 [0,1] 之间。如果某一能量处有不止一个透射通道，那么这些通道的透射总和（即该能量处的总透射系数）可能大于1。

多数的透射本征值都很小，可以忽略。只有两个比较大的透射本征值最为重要。

进一步，可以对前两个透射能量峰处的透射本征态进行可视化。

• 重新打开 Script Generator，删去 TransmissionEigenvalues单元；
• 添加两个 Analysis → TransmissionEigenstates 单元。
• 将两个 TransmissionEigenstates 单元的能量分别设置为 -3 和 -1.5；quantum number保持为 0，对应最上面的一个eigenvalue。
• 通过 Job Manager 执行脚本。

所得的数据对象将保存在 dtb_analysis.nc 中， 可以用 与上面分析本征态的方法一样用 Viewer 分析（见下图）。

我们继续分析 2.4eV 处的透射。从前面分析可以看出，这个能量处有一个透射概率为 0.618 的本征态，这与我们从总透射系数图上读出的数值（1.45）有很大的差别。这是由于前面的总透射图给出的是k平均透射。现在我们继续计算此能量处的k分辨的透射系数。

• 重新打开 Script Generator，删去 TransmissionEigenstates 单元；
• 添加 Analysis → TransmissionSpectrum 单元，设置
  • 能量起始和终点同为 2.4 eV；
  • 能量点数为 1；
  • k点取样为 21×21；
  • 通过 Job Manager 执行脚本。

计算结束后，可以使用 Transmission Analyzer 分析、查看结果。

此k分辨的透射散点图中在 $(k_A,k_B)=(0.0,0.38)$ 处有显著的峰，透射系数接近 1。在前面的计算中，由于k点取样密度不够，因此透射被高估了（1.45），此处k点为 21×21 时，透射约为 1.1。

现在，我们来计算最强透射k点处的透射本征态，可以是 (0,-0.38) 或 (0,0.38)，这两个点是时间反演对称性下等价的。

• 在右侧散点图上点击 (0,-0.38) 或 (0, 0.38)；
• 在左下方点击 Eigenvalues，得到一组本征值列表。

• 选中第一个本征值，点击 Eigenstates 按钮。稍等片刻，对应的本征态就计算好了，Viewer 窗口自动弹出。适当调整 isovalue 数值，得到漂亮的图。

这一节，我们来计算整个器件区域的局域器件态密度并在实空间的 Z 轴上进行投影。

• 重新打开 Script Generator，删去所有之前的分析单元；
• 添加一个 Analysis → ProjectedLocalDensityOfStates 单元；
• 双击打开这个单元，进行如下设置：能量窗口设为 -3 eV 到 3 eV，k点为 3×3；
• 将输出文件设为 pldos_0.0V.nc，脚本保存为pldos_0.0V.py。

执行脚本（计算可能耗时1小时）结束后，PLDOS 即出现在 VNL LabFloor 里。使用右侧 Projected Local Density Of States 工具即可给出 PLDOS如下图：

PLDOS 图显示了金电极部分很高的态密度（粉色），但是在中间真空区域里 DOS 几乎为零，只有在分子处有三个局域化的峰，分别位于 -3、-1.4 和 + 2.5 eV 处，这与透射谱的峰是吻合的。

这一章，我们正式计算 Au-DTB-Au 分子器件的 伏安特性（IV曲线），并分析 1.6V 偏压处的器件特性。

打开一个新的 Script Generator 窗口，进行如下操作：

• 添加 Analysis from File 单元；
• 添加 Analysis → IVCurve 单元；
• 将输出文件名改为 IV-curve.nc。

编辑 Analysis from File 单元：

• 设置 NetCDF 文件为 au-dtb-au.nc，使用 gID000 从零偏压计算结果继续计算。

编辑 IVCurve 单元：

• 设置偏压扫面范围从 V0=-2V 到 V1=2V；
• 设置电压点数为 11；
• 设置能量范围为 (-3,3)；
• 横向k点设置为 3×3；
• 选择 Krylov 自能计算器（为了加快计算速度）。

脚本设置结束，将其保存为 iv-curve.py 并通过 Job Manager 或命令行执行。这个计算耗时比较久，因为对每个偏压点（本计算共 11 个偏压点）都要进行自洽计算并分析透射谱。使用 12 核并行时，计算耗时约6小时。

计算结束后，文件夹里将出现 IV-curve.nc 文件，此外还有 11 个前缀为 ivcurve 的输出文件，它们保存了每个偏压点得自洽结果。

选中 IV-curve，点击右侧工具栏 IV-Plot 分析工具，I-V 曲线分析工具即弹出。

勾选 Additional plots，将显示三个新增面板，分别是微分电导（dI/dV），透射谱，谱电流。

要研究 偏压为 1.6 V 时分子器件的电子态性质，可以计算投影局域态密度（PLDOS），和零偏压下类似。

• 打开一个新的 Script Generator 窗口，添加两个单元：
  • Analysis from file
  • Analysis → PorjectedLocalDensityOfStates
• 在 Analysis from file 单元里选择 ivcurve_selfconsistent_configuration_1.60000V.nc；
• 在 PLDOS 分析单元里，将能量窗口 设为 -3 到 3 eV，设置 k 点为 3×3；
• 将输出文件名设置 为 pldos_1.6V.nc，保存脚本为 pldos_1.6V.py。

运行脚本后，使用 Projected Local Density of States 工具显示结果，如下图。

和零偏压相比，-1.5 eV 处的共振分裂成了两个，分别位于 -0.8 和 -2 eV，这个现象在透射谱上也可以看到。此外，我们还可以注意到 -0.8 的峰接近右侧电极，-2.0 eV 处的峰接近左侧电极。这种分裂来源于分子区域的电势降落（有外加偏压时）。

现在我们来分析偏压 1.6 V 时的电压降和诱导电荷密度。首先我们计算 1.6 V 和 0 V 下的 ElectrostaticDifferencePotential 和 ElectronDifferenceDensity。

• 在 Script Generator 中，添加 三个计算单元：
  • Analysis from File
  • Analysis → ElectrostaticDifferencePotential
  • Analysis → ElectronDifferenceDensity
• 在 Analysis from File 里设置 ivcurve_selfconsistent_configuration_1.60000V.nc；
• 将默认输出文件名改为 voltage_drop_16V.nc；
• 将脚本发送至 Job Manager

重复一次这个计算，但使用 0V 的自洽结果，保存文件为 voltage_drop_0V.nc。

为使用 1.6V 和 0V 的 ElectrostaticDifferencePotential 和 ElectronDifferenceDensity 结果计算电压降，我们使用右侧的GridOperations工具。

• 按住 ctrl 键同时选中 voltage_drop_16V 和 voltage_drop_0V 中的 $\delta V_E(\vec(r))$ 图标。
• 点击打开 GridOperations 工具；
• 设置函数为 A1 - B1；
• 保存为 EDP.nc。

重复上述步骤，计算电压诱导的电荷密度。

现在计算结果分别保存于 EDP.nc 和 EDD.nc。要图示这两个结果：

• 选择 EDP.nc 文件，打开 Viewer；
• 设置作图模式为 CutPlanes；
• 将结构（DeviceConfiguration）拖动到此打开的窗口里；
• 选择 EDD.nc 文件，也拖动到这个窗口里，使用默认的isosurface作图模式。

在 Property 里详细调节图示的细节。

如果设置如上图，则得到下面的图：

将电压降进行一维作图通常比较有用：选中 EDP 中的 $f(x)$ 图标，点击右侧的 1D Projector 工具。投影方法改为 Average，点击 Add Line。

此处的一维投影作图，C轴取值范围为 0 ~ 1，对应了真实的 z 方向长度。

• K. Stokbro, J. Taylor, M. Brandbyge, J.-L. Mozos, and P. Ordejón. Theoretical study of the nonlinear conductance of di-thiol benzene coupled to au(1 1 1) surfaces via thiol and thiolate bonds. Computational Materials Science, 27(1–2):151 – 160, 2003. doi:10.1016/S0927-0256(02)00439-1.
• K.S. Thygesen and K.W. Jacobsen. Molecular transport calculations with wannier functions. Chemical Physics, 319(1–3):111 – 125, 2005. doi:10.1016/j.chemphys.2005.05.032.
• 英文原文：http://docs.quantumwise.com/tutorials/molecular_device/molecular_device.html