对于周期性固体,薛定谔方程 $H \psi_{n{\bf k}} = E_{n{\bf k}} S \psi_{n{\bf k}}$ ($S$ 为重叠矩阵)中的 $\psi_{n{\bf k}}$ 可以写为 $\psi_{n{\bf k}}({\bf r}) = e^{-i {\bf k}\cdot {\bf r}} U_{n{\bf k}}({\bf r})$,这里的 $U_{n{\bf k}}({\bf r})$ 是与晶体自身周期性相同的周期函数。在一般的能带结构计算中,波矢量 ${\bf k}$ 为实数,通过求解上面的薛定谔方程得到不同 ${\bf k}$ (通常位于第一布里渊区的高对称线上) 值上的本征矢量,由此确定本征能量 $E_{n{\bf k}}$ (即能带结构)。
计算复数能带的方法可参考 [CS82]。在固定能量 $E$ 的情况下,要求解的是满足薛定谔方程的 ${\bf k}$ 值。这种解法可以得到实数和复数的 ${\bf k}$,实数 ${\bf k}$ 的解是通常的 Bloch 态,而带有虚部的解是在一个方向上呈指数递减,相反方向上递增。这样的解通常不能稳定存在于块体材料中,因此在能带结构计算中通常被忽略。然而,它们可以存在于表面或界面处,并且可以提供关于电子态如何在材料中衰减的信息,例如电子态如何通过一个薄的隧道势垒。
更多关于复数能带结构的概念和在导电结中的应用,请参考:
您将主要使用图形用户界面 QuantumATK 来设置和分析结果。如果您是 QuantumATK 的新手,建议您阅读 Basic QuantumATK Tutorial。
本教程中的计算将使用 QuantumATK 的半经验模型。所有参数的完整说明和有关其物理相关性在很多情况下的详细讨论,可参阅 ATK Reference Manual。特别是,参考手册条目中与复数能带计算相关的:ComplexBandstructure。
本教程使用特定版本的QuantumATK创建,因此涉及的截图和脚本参数可能与您实际使用的版本略有区别,请在学习时务必注意。
启动 QuantumATK,创建一个名为“ComplexBandstructure”的新项目。然后点击 Open,点击工具栏的按钮 开启 Builder。
在 Builder 里,点击 Add From Database,搜索“Silicon(alpha)”结构。双击结果所在行或者点击右下角的 按钮将结构添加到 Stash。
在面板中展开 Builders,打开 Surface(Cleave)工具,按照以下步骤设置 Si(100) 表面:
为缩短计算时间,避免区域折叠,我们通过以上设置和最小的层数表示表面。
切割平面即本例中的(100)面总是跨越两个第一原胞矢量和 ${\bf B}$。在“electrode”模式下,平面的法线与第三原胞矢量 ${\bf C}$ 重合,但在“bulk-like”模式下则不同。在 QuantumATK 中,复数能带始终沿第三倒易矢量 ${\bf g}_C$ 计算得到,其显然平行于 ${\bf A} \times {\bf B}$。因此使用当前结构您将获得沿(100)面的复数能带结构。
点击 Builder 右下角的 按钮将结构发送到 Script Generator ,执行以下步骤。
si_100_cbs.nc
。接下来,打开嵌入的 New Calculator 模块(双击),设置一个紧束缚计算。
下一步,打开 ComplexBandstructure 分析模块,然后编辑:
在复数能带计算中, ${\bf k}$ 在切割平面的投影保持不变,且投影值从 $(k_A,k_B)$ 参数获得,该参数用两个第一倒易晶格矢量 ${\bf g}_A$ 和 ${\bf g}_B$ 表示。因此,所得解将位于与切割平面法线平行的第三倒易晶格矢量 ${\bf g}_C$。对应于每一个 $(k_A,k_B)$ 值可以获得一组新解 $k_C+i\kappa_C$。
请注意,$k_C$ 是复数能带机构的实部,而 $kappa_C$ 为虚部。
您现在已经完成了 Python 脚本,保存为 si_100_cbs.py
并发送到 Job Manager 执行。像这种小的体系仅需要几分钟。如有需要,您也可以在此处下载脚本:↓si_100_cbs.py。
这种类型的计算并行得非常好,因此对于较大的结构,强烈建议并行执行脚本。
文件 si_100_cbs.nc
应该已经出现在了 Quantum 的 LabFloor。它包含了保存的 Hückel 计算和 ComplexBandstructure 分析数据块。
选择那个分析数据块,选择右边插件栏的 Show 2D Plot 工具,显示复数能带图的窗口弹出。放大一些可以过滤掉那些复数部分 $\kappa_C$ 的较大值,因为这些都不是真正相关的。
图51 Si(100) 的能带结构。图右侧显示的是实数能带——注意 ~1.2 eV 的间接带隙。图左侧显示的是复数能带随虚部 $\kappa_C$ 变化。具有一个小虚部的能带是最重要的,因为在您试图沿(100)方向上用电子穿过硅薄片时,它们会衰减较少。
在上图中,$\kappa_C$ 的解可从倒易的笛卡尔坐标中获取,以便更容易地比较不同结构。另一方面,对于能带结构的实部,$\kappa_C$ 的解由 $L$ 归一化后得到,分层垂直于切割平面。在 fcc 晶体沿(100)面切割的案例中,分层为 $L=a/2$,$a$ 为晶格常数。分层可用 layerSeparation() 方法从 ComplexBandstructure 数据块中提取获得,具体可参见 ComplexBandstructure。
提示:在较新的版本中,复数能带可以直接在分析工具中用Complex Band Structure Analyzer可视化,无需再使用脚本。此处保留脚本供大家学习如何使用Python语言进行复杂的作图。
在能带结构的复数部分中,解也包含了实数部分。因此解可以在三维图中可视化,$(x,y,z)=(k_C, \kappa_C, E)$。使用以下脚本可实现此目的。打开 QuantumATK 的 Editor ,复制粘贴脚本并保存为 3D-plot.py
。确保缩进正确。或者,您也可以在此处下载:↓ 3D-plot.py。
1 from QuantumATK import * 2 from mpl_toolkits.mplot3d import Axes3D 3 import matplotlib.pyplot as plt 4 import math 5 6 # Read the complex bandstructure object from the NC file 7 cbs = nlread('si_100_cbs.nc', ComplexBandstructure)[-1] 8 energies = cbs.energies().inUnitsOf(eV) 9 k_real, k_complex = cbs.evaluate() 10 L = cbs.layerSeparation() 11 12 fig = plt.figure() 13 ax = fig.add_subplot(111, projection='3d') 14 15 # First plot the real bands 16 kvr = numpy.array([]) 17 e = numpy.array([]) 18 19 # Loop over energies, and pick those where we have solutions 20 for (j, energy) in enumerate(energies): 21 k = k_real[j]*L/math.pi 22 if len(k)>0: 23 e = numpy.append(e,[energy,]*len(k)) 24 kvr = numpy.append(kvr,k) 25 26 # Plot the bands with red 27 ax.scatter([0.0,]*len(kvr), kvr, e, c='r', marker='o', linewidths=0, s=10) 28 29 # Next plot the complex bands 30 kvr = [] 31 kvi = [] 32 e = [] 33 34 # Again loop over energies and pick solutions 35 for (j, energy) in enumerate(energies): 36 if len(k_complex[j])>0: 37 for x in numpy.array(k_complex[j]*L/math.pi): 38 kr = numpy.abs(x.real) 39 ki = -numpy.abs(x.imag) 40 # Discard rapidly decaying modes which clutter the plot 41 if ki>-0.3: 42 e += [energy] 43 kvr += [kr] 44 kvi += [ki] 45 46 # Plot the complex bands with blue 47 ax.scatter(kvi, kvr, e, c='b', marker='o', linewidths=0, s=10) 48 49 # Put on labels 50 ax.set_xlabel('$\kappa$ (1/Ang)') 51 ax.set_ylabel('$kL/\pi$') 52 ax.set_zlabel('Energy / eV') 53 54 plt.show()
使用作业管理器 Job Manager 或从终端执行脚本。输出图将类似于下图,但您的点将更加分散,因为下图是由 Hückel 计算产生的,包含 10001 个能量点而不是本例的 501 个。
图 52 Si(100) 面复数能带的 3D 可视化图,实数能带由红色圆点表示,复数部分由蓝色圆点表示。请注意复数能带是怎样与实数能带连接的。
另一种将复数能带结构实数部分可视化的方式就是利用颜色。脚本 ↓ 2D_plot.py 可以实现(您可以看到之后的脚本)。它的最后可能会有些复杂,但那部分仅用于放置颜色条,可以省略。
1 from QuantumATK import * 2 import matplotlib.pyplot as plt 3 import math 4 5 # Read the complex bandstructure object from the NC file 6 cbs = nlread('si_100_cbs.nc', ComplexBandstructure)[-1] 7 energies = cbs.energies().inUnitsOf(eV) 8 k_real, k_complex = cbs.evaluate() 9 10 ax = plt.axes() 11 cmap="Spectral" 12 13 # First plot the real bands 14 kvr = numpy.array([]) 15 e = numpy.array([]) 16 for (j, energy) in enumerate(energies): 17 k = k_real[j]*cbs.layerSeparation()/math.pi 18 if len(k)>0: 19 e = numpy.append(e,[energy,]*len(k)) 20 kvr = numpy.append(kvr,k) 21 22 # Plot 23 ax.scatter(kvr, e, 24 c=numpy.abs(kvr), 25 cmap=cmap, marker='o', linewidths=0, s=10) 26 27 # Next plot the complex bands 28 kvr = numpy.array([]) 29 kvi = numpy.array([]) 30 e = numpy.array([]) 31 32 for (j, energy) in enumerate(energies): 33 if len(k_complex[j])>0: 34 kr = [numpy.abs(x.real) for x in k_complex[j]] 35 ki = [numpy.abs(x.imag) for x in k_complex[j]] 36 e = numpy.append(e,[energy,]*len(kr)) 37 kvr = numpy.append(kvr,kr) 38 kvi = numpy.append(kvi,ki) 39 40 # Plot with color depending on the imaginary part (corresponding to real k-points) 41 sc = ax.scatter(-kvi, e, 42 c=kvr, 43 cmap=cmap, marker='o', linewidths=0, s=10) 44 45 # Put on labels and decorations 46 ax.axvline(0,color='b') 47 ax.grid(True, which='major') 48 ax.set_xlim(-1, 1) 49 ax.set_ylim(-15, 10) 50 ax.annotate('$\kappa$ (1/Ang)', xy=(0.25,-0.07), xycoords="axes fraction", ha="center") 51 ax.annotate('$kL / \pi$', xy=(0.75,-0.07), xycoords="axes fraction", ha="center") 52 ax.set_ylabel('Energy / eV') 53 54 # Add a colorbar 55 fig = plt.gcf() 56 x1, x2, y1, y2 = 0., 1, ax.get_ylim()[0], ax.get_ylim()[0]+1 57 trans = ax.transData + fig.transFigure.inverted() 58 ax_x1, ax_y1 = trans.transform_point([x1, y1]) 59 ax_x2, ax_y2 = trans.transform_point([x2, y2]) 60 ax_dx, ax_dy = ax_x2 - ax_x1, ax_y2 - ax_y1 61 cmap_axes = plt.axes([ax_x1, ax_y1, ax_dx, ax_dy]) 62 a = numpy.outer(numpy.arange(0,1,0.01),numpy.ones(10)).transpose() 63 cmap_plt = plt.imshow(a,aspect='auto',cmap=plt.get_cmap(cmap),origin=(0,0)) 64 ax = plt.gca() 65 ax.set_xticks([]) 66 ax.set_yticks([]) 67 ax.set_xticklabels([]) 68 ax.set_yticklabels([]) 69 70 plt.show()
图 53 复数能带结构的 2D 可视化图。注意 k 值的颜色编码如何应用于能带结构的实数和复数部分,这使得能更容易分辨出复数能带附着于实数能带的位置。
复数能带的“森林”中有一个相当大的 $\kappa$ 值,在紧束缚模拟中并不常见。然而,对于 DFT 和 Hückel 理论中,基组更大,因此与含有各种价带的未占用能级连接的复数能带也就更多。
您可能注意到了可视化中存在着“间隙”。原因在于,与正常的能带图不同,不是将数据绘制成线而是圆点。在标准的能带图中,您可以用一些置信水平定义“能带”,在对称点(只有在能带交叉处会产生一些小问题)间连续分布。在当前案例中,首先解的数量在每个能量处(尤其是复数一侧)不同,且依赖于能量取样的密度,您可能不会找到一条特别的能带接近于某些点,这些点位于能带非常平缓(重的有效质量)的位置。通过增加能量点的数量可以在一定程度上缓解这种情况。