这是本文档旧的修订版！
目录
Meta-GGA 和二维受限的砷化铟
版本：2017.dev
在本教程中，您将学习两个主题：（i）通过密度泛函理论（DFT）结合元广义梯度近似（meta-GGA）交换关联（XC）函数计算能带结构（ii）受限结构的计算，重点专注于如何钝化表面的悬空键。完成本教程后，您将能够：
- 使用 meta-GGA XC 函数执行计算，获得精确的半导体带隙；
- 设置并钝化受限的 III-V 体系（InAs，GaAs，InP，……）;
- 采用具有非整数电荷的特殊赝势钝化氢，即钝化 III-V 和 II-VI 体系表面所需的人工氢原子或赝氢原子，这是因为正常的氢钝化会导致金属表面态;
- 对钝化（正常）的氢原子加上补偿电荷可以替代赝氢原子。
注意
如果您还不熟悉 QuantumATK，请先阅读基本教程。本教程所用基础计算引擎是 ATK-DFT，您可以在 QuantumATK 参考手册中找到其详细信息。
Meta-GGA
meta-GGA 简介
当由 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 关联。
在 QuantumATK 里的 Meta-GGA
在 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 参数运用于受限结构。
采用 meta-GGA 计算块体 InAs 的能带结构
设置块体 InAs
现在您应该设置一个 InAs 块状晶体。在 Builder，点击 Add 
Database…，在搜索区域输入 InAs，添加结构到 Stash。
下一步是更改晶格常数到实验上室温时的晶格常数 $a=6.0583 Å$ [3]。
- 在控制板的右侧，打开 Bulk tools
Lattice Parameters… - 重要的是，您必须选择保持 fractional 坐标不变（如下图所示），否则晶体的对称性将在晶格常数改变时会被破坏。
- 设置 $a=6.0583 Å$。
接下来，将结构发送到 Script Generator。
设置脚本
在下文中，您将使用 meta-GGA 泛函和自洽及非自洽计算得到的 c参数计算 InAs 块体的能带结构。
为实现此目的：
- 添加 New Calculator 和 Bandstructure 模块。
打开 New Calculator 模块，做出如下更改：
- 设置 Exchange correlation 为 meta-GGA。
- 设置 k 点取样为 7,7,7。
- 在 Basis set/exchange correlation 下，更改砷的赝势为“HGH[Z=5], LDA.PZ”，铟的赝势为“HGH[Z=13], LDA.PZ”。为两种元素都选择 Tier 3 基组。
- 最后设置输出文件名为 InAs_hgh.nc。
保存脚本，然后使用 “Send To” 按钮将其发送到 Job Manager，运行作业。该计算将会花费 1 或 2 分钟。
注意
对于 c 参数的自洽计算，为得到更精确的能带结构通常使用包含半芯态的 HGH 赝势。这是针对铟 HGH [Z = 13] 的情况，它包括 10个 4d 电子，而 HGH [Z = 3] 仅包括 3 个价电子（5s 态和 5p 态）。然而，在 ATK 中仅提供一些元素的半芯态赝势；对于砷，只有一个 HGH [Z = 5] 赝势，没有伪芯态。
自洽确定 c 参数的能带结构
当计算结束后，找到 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 参数是自洽计算得到的。
拟合 meta-GGA c 参数
虽然 HGH 赝势通过结合自洽计算的 c 参数通常能给出相当准确的带隙，但我们想要得到能够与实验值更加一致的结果。为了实现这一点，可以手动设置 meta-GGA 的 c参数使得计算的带隙与实验带隙一致。这也可以作为如本节所示补偿稍微不准确的赝势或基组的方法。
设置脚本
返回到 Script Generator，打开 New Calculator 模块并做如下设置：
- 在菜单选项下：Basis set exchange correlation，更改砷的赝势为“FHI[Z=5], LDA.PZ”，铟的赝势为“FHI[Z=3], LDA.PZ”。也就是说选择没有半芯态的赝势。
- 给两种元素都设置基组为 DoubleZetaPolarized。
- 更改输出文件名称为“InAs_c_fit.nc”。
点击 “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）得到优化。
采用优化后 c 参数的能带结构
在继续研究受限 InAs 结构之前，您将首先使用优化的 c 参数更详细地计算和分析块体能带结构。
返回到 Script Generator，作如下调整：
- 删除 Bandstructure 数据块。
- 插入一个 LocalBandstructure 模块，打开它（如下所示）。
- 设置“Bands below Fermi level”为 0。
- 设置“Points per segment”为 50。
- 增加第二个方向：0.5，0.5，0。
- 其他保持为默认值。
- 更改输出文件名称为 “InAs_bulk_c_optimized.nc”。
将脚本转移到 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$ 为各向异性。
设置并钝化 InAs 平板
平板结构
现在您将设置一个在 [100] 方向上受限的 InAs 平板结构。
返回至 Builder 窗口，假设您仍拥有之前的块状 InAs 晶体。
下一步，点击 Builders Surface (Cleave)…，沿 (100) 方向切割晶体。
- 保持默认的密勒指数 (100)，按下“Next”。
- 保持默认的表面晶格矢量，按下“Next”。
- 保持默认的超胞——这将确保线的方向垂直于表面。
- 选择“Non-periodic and normal (Slab)”，厚度为 8 层，顶部和底部的真空为 5 Å 的构型为最终结构。
- 按下“Finish”将切割后得结构添加到 Stash。
加氢
由于未配对的 In 原子和 As 原子，所谓裸露的 InAs 平板表面会有悬空键。这将导致在带隙中带有能量的金属表面态。为了钝化表面，您将现在添加氢原子。
打开 Coordinate Tools Custom Passivator。
设置 Hybridization（杂化）为“4 (SP3)”。
更改键长为 1.3 Å（否则相邻的氢原子会距离太近）。
按下“Passivate”按钮。
注意
Custom Passivator 为添加的氢原子加了标签 —— 这对本教程至关重要，因此请勿更改或删除它们。您可以通过打开 Select Tags 来检查标签。在指定补偿电荷和使用特定的赝势时，您将会稍后使用这些标签辨认分别与 In 和 As 结合的H原子。
该结构现在应该类似于下图。
发送结构到 Script Generator。
采用默认氢原子的能带结构
在下文中，您将对钝化氢原子采用默认设置计能算带结构。为实现此目的：
- 添加一个 New Calculator。
- 添加 Analysis/Bandstructure。
打开 New Calculator 模块，设置交换关联为 meta-GGA，k 点取样为 7，7，1。注意 FHI 赝势是默认的选项。
最后，打开 Bandstructure 模块，设置布里渊区路径为“Y，G，X”，增加每段的点数至 100。
正如在 meta-GGA 的简介里所讨论的，c 参数的自洽计算仅适用于没有真空区域的块体系统。然而对于平板结构的计算，我们应该采用一个固定值。因此将脚本发送到 Editor，然后指定您在优化 c 参数的能带结构中找到的 c 参数使 FHI 赝势重现实验带隙。
#---------------------------------------- # Exchange-Correlation #---------------------------------------- exchange_correlation = MGGA.TB09LDA(c=0.936)
保存脚本，将其发送到 Job Manager（再次使用 “Send To” 按钮）并运行计算。本次耗时略长，约 5 分钟。
注意
运行脚本时，您可以能会在 log 文件中发现一个警告：
################################################################################ # WARNING # # # # The computed TB09 meta-GGA XC potential # # diverged in 0.01 % of the simulation volume, and was # # truncated to be in the range [-10.0000, 10.0000] Hartree # # # ################################################################################
您无需为此担心，因为您已将 c 参数设置为了固定值。另一方面，如果您试图通过自洽计算它，则警告将是相关的。然而，如在 meta-GGA 的简介中所解释的，c 参数的自洽评估涉及包括真空区域的整个计算单胞的积分，其中 meta-GGA XC 势可以分离。因此，对于平板和纳米线等受限系统，您始终需要手动指定 c 参数。