版本：2017.dev

在本教程中，您将学习两个主题：（i）通过密度泛函理论（DFT）结合元广义梯度近似（meta-GGA）交换关联（XC）函数计算能带结构（ii）受限结构的计算，重点专注于如何钝化表面的悬空键。完成本教程后，您将能够：

• 使用 meta-GGA XC 函数执行计算，获得精确的半导体带隙；
• 设置并钝化受限的 III-V 体系（InAs，GaAs，InP，……）;
• 采用具有非整数电荷的特殊赝势钝化氢，即钝化 III-V 和 II-VI 体系表面所需的人工氢原子或赝氢原子，这是因为正常的氢钝化会导致金属表面态;
• 对钝化（正常）的氢原子加上补偿电荷可以替代赝氢原子。

当由 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)

在自洽计算之后。

现在您应该设置一个 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 分钟。

当计算结束后，找到 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 模块并做如下设置：

• 在菜单选项下：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）得到优化。

在继续研究受限 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$ 为各向异性。

现在您将设置一个在 [100] 方向上受限的 InAs 平板结构。

返回至 Builder 窗口，假设您仍拥有之前的块状 InAs 晶体。

下一步，点击 Builders Surface (Cleave)…，沿 (100) 方向切割晶体。

• 保持默认的密勒指数 (100)，按下“Next”。
• 保持默认的表面晶格矢量，按下“Next”。
• 保持默认的超胞——这将确保线的方向垂直于表面。
• 选择“Non-periodic and normal (Slab)”，厚度为 8 层，顶部和底部的真空为 5 Å 的构型为最终结构。
• 按下“Finish”将切割后得结构添加到 Stash。