这是本文档旧的修订版!
版本:2017.dev
在本教程中,您将学习两个主题:(i)通过密度泛函理论(DFT)结合元广义梯度近似(meta-GGA)交换关联(XC)函数计算能带结构(ii)受限结构的计算,重点专注于如何钝化表面的悬空键。完成本教程后,您将能够:
如果您还不熟悉 QuantumATK,请先阅读基本教程。本教程所用基础计算引擎是 ATK-DFT,您可以在 QuantumATK 参考手册中找到其详细信息。
当由 Kohn-Sham 能量本征值计算带隙时,DFT 因预测值过小而不受欢迎。然而,meta-GGA 泛函的研究进展使我们能够计算得到准确的带隙。
Meta-GGA 泛函属于所谓 XC 函数的 Jacob 阶梯的第三阶,它不仅包含局域密度 $\rho(\mathbf{r})$ (如在 LDA 中,第一阶)和密度梯度 $\nabla\rho(\mathbf{r})$ (如在 GGA 中,第二阶),而且还有动能密度 $\tau(\mathbf{r})$。在 2009 年,Tran 和 Blaha[1] 表明可以获得各种材料的准确带隙。在他们的公式中,交换式可由以下给出:
$$v_x^{TB}(\mathbf{r}) = cv_x^{BR}(\mathbf{r}) + \frac{3c-2}{\pi}\sqrt{\frac{4\tau(\mathbf{r})}{6\rho(\mathbf{r})}},$$
在这里,$\tau(\mathbf{r})=1/2\sum_{i=1}^N|\nabla\psi_i(\mathbf{r})|^2$ 是动能密度,$\psi_i(\mathbf{r})$ 是 Kohn-Sham 轨道,$v_x^{BR}(\mathbf{r})$ 是 Becke-Roussel 交换势Blaha[2]。c 参数 可以由下式计算得到:
$$c = \alpha +\beta\left[\frac{1}{\Omega}\int_\Omega\frac{|\nabla\rho(\mathbf{r})|}{\rho(\mathbf{r})}d\mathbf{r}\right]^{1/2},$$
此处,$\alpha=-0.012$ 和 $\beta=1.023\ \mathrm{Bohr}^{1/2}$ 是通过拟合复制大量半导体和绝缘体的实验带隙确定的,${\Omega}$ 为晶胞体积。Tran 和 Blaha 采用的关联势是一个普通的 LDA 关联。
在 AT K中,您可以通过两种方式使用 Tran和Blaha 采用的 XC 泛函。您可以将 XC 指定为
exchange_correlation = MGGA.TB09LDA
在这种情况下,c 参数可以基于以上表达式自洽地地确定。这是默认设置,脚本的方式是通过 QuantumATK 里的 Script Generator 设置。
您也可以手动设置 c参数的值:
exchange_correlation = MGGA.TB09LDA(c=1.0)
电子和动能密度仍然完全自洽计算,但使用固定的 c 值。
如果以自洽方式计算 c 参数,则该值将被记录在日志文件中。它也可以在脚本中按照如下方法被提取
c = calculateTB09C(bulk_configuration)
在自洽计算之后。
使用 c 参数的自洽确定应仅针对在所有方向上具有周期性的块体构型。受限的系统如平板或纳米线将包含较大真空区域。由于 c 参数被作为整个晶胞体积的积分来计算,因此来自真空区域的贡献将会导致不正确甚至不同的结果。如果尝试此操作,您很可能会在日志文件中看到警告消息。
对于受限系统,必须首先为相应的块体系统确定适当的 c值。可以通过自洽或通过将 c 参数拟合到例如实验带隙(参见 Fitting the meta-GGA c-parameter)来实现。然后将由此确定的 c 参数运用于受限结构。
现在您应该设置一个 InAs 块状晶体。在 Builder,点击 Add Database…,在搜索区域输入 InAs,添加结构到 Stash。
下一步是更改晶格常数到实验上室温时的晶格常数 $a=6.0583 Å$ [3]。
接下来,将结构发送到 Script Generator。
在下文中,您将使用 meta-GGA 泛函和自洽及非自洽计算得到的 c参数计算 InAs 块体的能带结构。
为实现此目的:
打开 New Calculator 模块,做出如下更改:
保存脚本,然后使用 “Send To” 按钮将其发送到 Job Manager,运行作业。该计算将会花费 1 或 2 分钟。
对于 c 参数的自洽计算,为得到更精确的能带结构通常使用包含半芯态的 HGH 赝势。这是针对铟 HGH [Z = 13] 的情况,它包括 10个 4d 电子,而 HGH [Z = 3] 仅包括 3 个价电子(5s 态和 5p 态)。然而,在 ATK 中仅提供一些元素的半芯态赝势;对于砷,只有一个 HGH [Z = 5] 赝势,没有伪芯态。
当计算结束后,找到 LabFloor 上的能带结构数据块,点击 Bandstructure Analyzer 画出如下图的能带结构。
您会发现计算的带隙非常接近实验值 0.354 eV。当使用 HGH 赝势时,这实际上是一个相当普遍的趋势。下图显示了一系列半导体和绝缘体的计算带隙与实验带隙。图中展示了 LDA 和 meta-GGA 两种的计算结果(其他实例可参见 Ref[1]),也显示了 LDA 总是低估带隙(在 DFT 中众所周知的带隙问题),而 meta-GGA 的结果通常非常接近实验值,且通常与代价更高的混合泛函或甚至 GW 计算的结果一样准确。
图 75 计算带隙与实验带隙。采用 LDA(蓝色方块)和 meta-GGA(红色圆圈)XC 泛函进行的计算。对于 meta-GGA 的计算,HGH 赝势常与 Tier 4 基组一起使用,并且 c 参数是自洽计算得到的。
虽然 HGH 赝势通过结合自洽计算的 c 参数通常能给出相当准确的带隙,但我们想要得到能够与实验值更加一致的结果。为了实现这一点,可以手动设置 meta-GGA 的 c参数使得计算的带隙与实验带隙一致。这也可以作为如本节所示补偿稍微不准确的赝势或基组的方法。
返回到 Script Generator,打开 New Calculator 模块并做如下设置:
点击 “Send To” 按钮将脚本发送到 Editor。
在 Editor 里,找到指定 XC 函数的所在行,将 c 参数这个自变量添加到 MGGA.TB09LDA 泛函。
#---------------------------------------- # Exchange-Correlation #---------------------------------------- exchange_correlation = MGGA.TB09LDA(c=1.0)
然后将脚本发送到 Job Manager(再次使用 “Send To” 按钮)
计算结束后,画出能带结构图,测量带隙值。如果带隙比实验值(0.354 eV)大,就返回到 Editor 将 c 值降低一点,再次运行计算。如果带隙值较小,就必须增大 c 值。可以尝试一些迭代观察是否能得到与实验值一致的结果。
下图展示了在 4 个不同 c 值(0.9,0.95,1.0,1.05)下计算得到的带隙(左图)。由于特有的线性行为,我们可以安全地在计算带隙的插值线和实验带隙的交点处找到最佳 c 参数。于是就产生了我们将在接下来使用的 c 值 0.936。
右图为同样这 4 个 c 参数下计算得到的导带有效质量。还显示了实验的导带有效质量值(0.023 $m_e$)。线之间的交点在 c = 0.9356 处,所以非常重要的是,近似相同的 c 值可以优化带隙和有效质量,因此至少可以给 InAs 的最低 $\Gamma$ 凹谷提供一个非常好的描述。
图 76 计算的带隙(左)和导带有效质量(右)与meta-GGA c 参数的关系图。增加 c 值将始终导致带隙逐渐变大。带隙和有效质量都通过近似相同的 c 参数(0.936)得到优化。
在继续研究受限 InAs 结构之前,您将首先使用优化的 c 参数更详细地计算和分析块体能带结构。
返回到 Script Generator,作如下调整:
将脚本转移到 Editor,插入优化的 c 参数值:
#---------------------------------------- # Exchange-Correlation #---------------------------------------- exchange_correlation = MGGA.TB09LDA(c=0.936)
最后,发送脚本到 Job Manager,运行计算。
为了分析结果,可以采用 Local Bandstructure Analyzer。
LocalBandstructure 数据块计算有效质量 $m^*$,并由函数 $E(1+\alpha E)=\frac{\hbar^2 k^2}{2m^*}$ 进一步拟合能带结构,该式中 $\alpha$ 是要拟合的参数。您可以使用适合的点的数量来控制拟合的范围。如下图所示,这种非抛物线色散关系比通常的抛物线模型 $E=\frac{\hbar^2 k^2}{2m^*}$ 能更好地拟合 DFT-meta-GGA 能带结构。红色的点表示用于拟合 $\alpha$ 的范围。
图 77局部能带结构分析。黑点显示计算的能带结构数据,而蓝色和红色曲线分别显示抛物线和非抛物线模型的结果。拟合的 $\alpha$ 值和有效质量在图下方显示。
从 LocalBandstructure Analyzer,我们可以看出 InAs 在 [1,0,0] 方向上的导带可由 $m^*=0.028\,m_e$ 和 $\alpha=2.869\,eV^{-1}$ 精确描述,而抛物线模型只能较好地描述导带最小值邻近的能带。
尝试更改方向为 $(0.5, 0.5, 0)\,Å^{-1}$,计算此方向上的局部能带结构。注意,有效质量和之前的相同,即它是各向同性,而抛物线参数 $\alpha$ 为各向异性。