用户工具

站点工具


atk:分子器件模拟

分子器件模拟

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

本教程中的计算时间很长,作为入门的QuantumATK量子输运计算教程,请参考:

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

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

零偏压计算

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

注意

这个几何结构未经优化。这对本教程中要讨论的概念没有影响,但对于实际研究中的计算,结构优化至关重要。要了解如何弛豫器件几何结构,请参见:如何弛豫器件体系的几何结构

使用Script Generator设置计算

选中器件结构,点击 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 启动计算。此计算在单机上大约耗时一小时。

并行计算可以大幅度提升计算速度,要了解如何并行计算,请参考:QuantumATK并行计算

零偏压结果分析

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

  • 将透射谱与苯环的分波态密度(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键用鼠标单击逐一选中)

将Atom Indices拷贝下来,后面步骤会用到。

PDDOS的峰和透射谱峰的相似性提示了苯分子能级和透射谱的对应关系。

分子投影自洽哈密顿量(MPSH)

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)。

说明

$\pi$-能级共有6个电子,下方三个为占据轨道。因此我们可以确定 14 号轨道为 HOMO,15 号轨道为 LUMO。

透射本征值和本征态

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

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

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

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

前面两个能量处只有一个透射本征态主导,最后一个能量处则有两个主导的透射本征值。后面会讨论为什么他们加在一起没有得到总的透射谱 T(E)。

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

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

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

透射本征态是复数的波函数。等值面图(isosurface)显示了绝对值,颜色则表明了相位。由于相位是周期的,因此最好使用周期的color map。

k分辨透射和本征通道

我们继续分析 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点采样点数对透射谱计算的影响的重要性。透射谱计算使用的k点数与进行自洽电子密度计算使用的k点数也没有关系。换句话说,某k点数可以收敛电子密度自洽计算,但同样的k点数可能无法给出准确的透射结果。

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

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

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

投影局域态密度 (PLDOS)

这一节,我们来计算整个器件区域的局域器件态密度并在实空间的 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 处,这与透射谱的峰是吻合的。

IV特性

这一章,我们正式计算 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特性分析

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

提示

IVCurve数据对象除了保存了IV曲线的信息外,还保存了每个偏压点的透射谱信息。

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

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

每条透射谱曲线都对应了一个偏压点。某一偏压下的能量窗口在透射谱曲线上用蓝色标识。用鼠标指向某一偏压点时(红色圆点),对应的透射谱和谱电流曲线将被点亮显示。

注意

要特别注意谱电流曲线:所有的对总电流的贡献都必须包含在能量窗口内,此时IV计算才是可靠的。如果实际情况并非如此,则需要使用更大的能量窗口重新进行IV计算。

dI/dV 在1.6V 处有峰,这是因为透射谱上的共振峰在这个偏压时开始进入偏压窗口。

1.6 V 偏压时的 PLDOS

要研究 偏压为 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_16Vvoltage_drop_0V 中的 $\delta V_E(\vec(r))$ 图标。
  • 点击打开 GridOperations 工具;
  • 设置函数为 A1 - B1;
  • 保存为 EDP.nc

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

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

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

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

如果设置如上图,则得到下面的图:

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

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

说明

两个图都显示了器件区域 1.6V 的电压降。很明显,电压降主要发生在分子附近,而不是电极区域。电压降对应了量子电阻。分子区域的电压降也解释了非零偏压下的共振峰的分裂。

参考

  • 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.
atk/分子器件模拟.txt · 最后更改: 2018/03/20 21:51 由 liu.jun

© 2014-2020 费米科技(京ICP备14023855号