版本：2016.04

在本教程中，您将学习如何计算中性和带电点缺陷的形成能。具体地，我们将以 GaAs 为例，因为它的结构中包含了 As 和 Ga 空位的各种电荷态。

您可能还需要参考形成能的入门教程，该教程侧重于块体系统和中性点缺陷：Calculation of Formation Energies。

下载 & 链接：

• ↓ PDF version
• ↓ Slides (PDF)
• ↓ bulk-As.py
• ↓ bulk-GaAs.py
• ↓ big-bulk-GaAs.py
• ↓ vacancy-calc.py
• ↓ formation-energies.py
• Basic QuantumATK Tutorial
• ATK Reference Manual

点缺陷 $X$ 的形成能 $E^f$ 被定义为研究系统和参考状态组分的能量差。类似于块状材料的形成能，但对于带电缺陷，如何选择参考状态则不是太清楚。在这里，我们遵循 C. Freysoldt 等人在综述中概述的程序 [FGH+14]。

$$E^f[X^q] = E_{tot}[X^q]-E_{tot}[bulk]-\sum_i n_i \mu_i + q \left(E_{VBM} + \mu_e \right) + E_{corr} .$$

$E_{tot}[X^q]$ 是带电荷 $q$ 的缺陷态总能量，其中所带电荷以基本电荷（$e>0$）的单位给出。$E_{tot}[bulk]$ 是与带缺陷结构大小相同的块体晶胞的总能量，$n_i \mu_i$ 是 $n_i$ 的参考能量与元素 $i$ 化学势 $\mu_i$ 的乘积。括号中的术语说明了对涉及到使缺陷带电的电子的化学势。$E_{VBM}$ 是由 ATK 计算研究块体（本教程中为无缺陷的块体 GaAs）的能带结构所给出的价带最大值，$\mu_e$ 在此处定义的是对应价带顶部的电子化学势。考虑到例如掺杂等原因引起费米能级的移位，$\mu_e$ 参数可以被当做一个自由参数。注意 $\mu_e=E_{gap}/2$ 对应的是未掺杂半导体的情况，$E_{gap}$ 是本征半导体带隙。最后的 $E_{corr}$ 是相关修正项之和。本例中，它校正了带电点缺陷周期性区域间的静电相互作用。由于库仑相互作用的长程特性，很难使用足够大的系统而不需要这种校正。想要了解有关修正方案的更多细节，请参阅下文的附录。

为了计算形成能，我们需要计算上面列出的每一项。流程如下：

• 选择要研究的化合物和缺陷，以及要使用的计算模型（交换关联、赝势等）。
• 采用所选计算模型弛豫原始的块体系统。
• 确定缺陷中元素的化学势。本例中，我们使用块体的总能量为移除的元素创造空位。
• 设置无缺陷的超胞计算 $E_{tot}[bulk]$ 和 $E_{VBM}$。
• 引入（带电的）缺陷，弛豫系统获得 $E_{tot}[X^q]$。
• 计算 $E_{corr}$ 中的关联项，最终得到形成能。

作为点缺陷的一个示例，我们将计算 GaAs 中 As 空位的形成能，这种空位能够支持许多不同的电荷态。但是我们在给缺陷加上电荷之前，首先关注一下中性空位点缺陷。在这种情况下，形成能方程简化为

$$E^f[V^0] = E_{tot}[V^0] - E_{tot}[As] -(-1) \mu _{As} .$$

前面两项很容易计算，但可能会耗时久。相反地，第三项是 As 的参考化学势，我们在这里的选择会直接影响点缺陷计算形成能的绝对值。这里最明显的选择块体 As，也是我们将要采用的，但其他选择可能在特定情况下可以合理的表征。此处我们就不深入研究了。

总能量的计算将采用 DFT 与 PBESol 交换关联函数，设计 GGA 近似适用于固体 [PRC + 08]。

我们首先需要弛豫 GaAs 和 As 块体的最小单胞和内部坐标。对于后者，我们还需要每个原子的总能量。使用含有 PBESol 的 DFT 计算器和 OptimizeGeometry 模块弛豫这两个块体上的压力和应力。这将在其他地方详细描述，如 Calculation of Formation Energies。为方便起见，我们还提供了以下脚本：

↓ bulk-GaAs.py

↓ bulk-As.py

在 LabFloor 上找到优化的 BulkConfiguration，它应该有个名为 gID001 的号码，并将其拖放至 Builder，设置一个优化了晶格常数的新 BulkConfiguration。您还可以在输出文本中找到优化的晶格参数。但我们仍然需要创建更大的晶胞，而且还希望是一个正交晶胞，因为它会简化静电校正。

• 为实现以上，请在 Stash 选中 GaAs 系统，然后转至右侧的 BulkTools Supercell，选择 Conventional 并点击 Transform。GaAs 现在由一个更大的立方晶胞所代表。
• 为使缺陷获得一个合理的结果，我们仍然需要将晶胞设置得更大。点击 BulkTools Repeat，将 3 个数字都增加到 2，点击 Apply。
• 现在将构型发送到 Script Generator 并添加 New calculator、Density of States 模块和 TotalEnergy 分析模块。
• 将交换关联函数更改为 PBES（即 PBESol），并将 k 点取样设置为 4x4x4，如此计算晶胞已足够大。将 Density of States 模块的 k 点取样更改为 8x8x8。
• 现在执行计算，大约需要 10 分钟。

↓ big-bulk-GaAs.py

为设置中性点缺陷的计算，您可以使用脚本 ↓ vacancy-calc.py 或自己构建脚本。若是后一种情况，需要返回 Builder 并在 Stash 中选择 2x2x2 的正交晶胞。然后执行以下操作：

• 选择一个 As 原子，最好是最接近原点的原子，然后点击顶部栏的 图标将其转换为鬼原子。
• 将构型发送到 Script Generator 并添加 New calculator、OptimizeGeometry 模块和 TotalEnergy 分析模块。
• 更改计算器设置以适合原始的 2x2x2 GaAs。其他两个模块采用默认设置就可以了。运行计算，可能需要几个小时进行结构优化。

我们得到无缺陷块体结构的总能 $E_{tot}[GaAs]$ = -7460.13 eV，含有中性 As 空位的构型总能 $E_{tot}[GaAs]$ = -7285.59 eV，然后结合上文 As 块体的结果，得知中性 As 空位的形成能是 3.22 eV。您也可以使用脚本 ↓formation-energies.py 计算 $E^f[V^0]$。

现在我们已经计算了中性点缺陷的形成能，接下来将把注意力转移到在建模中电荷的并入。如果我们再次聚焦 As 空位，则根据以下公式计算形成能：

$$E^f[V^q] = E_{tot}[V^q]-E_{tot}[GaAs]+ \mu_{As} + q \left(E_{VBM} + \mu_e \right) + E_{corr} .$$

第一项与我们刚刚对中性点缺陷的计算类似，而第二项和第三项与之前完全相同，不需要重新计算。价带最大值由无缺陷块体的态密度（DOS）计算确定，$\mu_e$ 将被视为自由参数。最后一项 $E_{corr}$ 根据 Freysoldt 等人 [FNVanDWalle09] 的方法计算，其中如附录所述，为说明因缺陷引入产生的偏移，也会对 VBM 做一个小的修正。请注意，校正项取决于带电缺陷周围块状材料的屏蔽性质。因此，如附录中进一步解释的那样，我们需要从实验或计算中获得无缺陷块体 GaAs 的介电常数值。

As 空位总共可以维持 5 个电荷态，包括上文研究过的中性状态。我们现在将计算这些点缺陷的形成能。对于非常大的系统，使用弛豫过的中性点缺陷作为电荷态的初始猜测可能是值得的，但对于本例使用的系统尺寸来说并不重要。

可以按照上述步骤在计算器中添加电荷，也可以直接复制脚本并在计算器模块中添加 charge = q。您还可以使用 ↓ vacancy-calc.py 脚本，它将大块体超胞的 .nc 文件作为输入，并在结构优化和计算总能量之前添加缺陷。

现在我们已经准备好了所有需要的部分，可以通过从保存的数据文件中提取总能量和计算修正项得到完整的形成能。为此，您可以使用脚本 ↓ formation-energies.py。它将把弛豫过的缺陷态 .nc 文件和包含有效 DOS 计算的同等大小的无缺陷超胞作为输入文件。 您还需要如教程：Optical Properties of Silicon 所示，利用您所选的计算模型估算介电常数。我们还在附录中展示了如何计算 GaAs 的介电常数。

表 7 含有 As 空位的 2x2x2 GaAs 超胞形成能（包括周期性和能带偏移校正）。

上面给出的形成能是根据 VBM 处的电子化学势计算的，即 $\mu_e = 0$ eV。考虑到 $\mu_e$ 作为自由参数能使费米能级移位，我们也可以为每个缺陷画出形成能与电子化学势的函数关系图。这也相当于随材料的掺杂水平而变化。结果如下图所示，在给定 $\mu_e$ 下，最稳定的电荷态即最低的线。由图可知，对于 p 型掺杂的 GaAs，在 $\mu_e = 0$ eV 时，带正电荷的 As 空位是最稳定的。随着电子化学势的逐渐增大，最稳定的电荷态从 $\mu_e > 0.35$ 开始变为中性，然后是负电荷态。尤其是本征（无掺杂） GaAs 在 $\mu_e = E_{gap}/2\sim 0.75$ eV 时，带有电荷 $q=-2$ 的 As 空位是最稳定的点缺陷。

图97 GaAs 中 As 空位各种电荷态的形成能随 $\mu_e$ 函数的变化（相对于 VBM 的定义），包含 VBM 为零的情况，计算采用了 2x2x2 正交晶系超胞。

如果我们将超胞尺寸增加到 3x3x3，k 点网格更改为 2x2x2，我们会得到下表列出的形成能，以及使用相同大小超胞的的文献值[EMM05]。

表 8 计算形成能（包括周期性和能带偏移校正）与文献值的比较

我们可以看到形成能有怎样的细微变化，中性点缺陷的偏移约为 0.04 eV。带电缺陷的差异略大，这是由于 k 点密度的小更改引起 VBM 的轻微偏移。本例没有与参考文献完全一致的，但应该注意的是，参考文献中没有使用鬼原子，采用了较陈旧的静电校正方法和 LDA 交换关联。

这些变化也会影响形成能和 $\mu_e$ 的函数关系图，尽管线的斜率相同，但截距略有不同，导致缺陷稳定性的整体图像发生显着变化。具体而言，中性状态现在仅在非常窄的区域稳定，且从 $\mu_e = 0.4$ 起带负电的状态最稳定。

图98 GaAs 中 As 空位各种电荷态的形成能随 $\mu_e$ 函数的变化（相对于 VBM 的定义），包含 VBM 为零的情况，计算采用了 3x3x3 正交晶系超胞。

正如上文已经提到的，我们会采用由 Freysoldt 等人 [FNVanDWalle09] 提出的校正方法修正因点缺陷引入造成的电荷周期图像间的静电相互作用和能带偏移。校正包含以下两项：

$$E_{corr} = E_{lat} - q \Delta V_{q/b} ,$$

$E_{lat}$ 修正了电荷周期图像间的静电相互作用，$q \Delta V_{q/b}$ 则修正了能带偏移。$E_{lat}$ 通过考虑模型电荷分布的静电作用和对比周期晶胞中边界条件为非周期性的势计算得到：

$$E_{lat} = \int_\Omega \left[\frac{1}{2} q (\tilde{V}_q^{lr}-V_q^{lr}) \right] d^3 r .$$

$V_q^{lr}$ 是模型电荷分布的长程势，$\tilde{V}_q^{lr}$ 周期晶胞的修正量。$\Delta V_{q/b}$ 由引入点缺陷造成的静电势移位计算得来:

$$ \Delta V_{q/b} = \tilde{V}_{q/b} - \tilde{V}_q^{lr} , $$

这里的 $\tilde{V}_{q/b}$ 是原始块体构型和添加缺陷间静电势的差值。为了得到零点参考能量的正确偏移量，应在远离点缺陷且是常数的位置计算 $\Delta V_{q/b}$。在我们的实际操作中，在离缺陷最远的五个网格点上取平均值。我们会在下图中画出这些值，$\Delta V_{q/b}$ 由黑线表示：

图 99 能带偏移校正图。绿色是缺陷的模型电位，红色是由于引入缺陷引起的变化，而黑色则是绝对能量零点的偏移。

注意，所获得的差异是在两个屏蔽势间定义的，参见文献 [FNVanDWalle09] 和脚本 ↓ forma-energies.py，其中屏蔽是因为带电缺陷周围的半导体材料引起的。描述电介质屏蔽的静电介电常数的计算在下文将以 GaAs 为例介绍。

为了说明半导体材料（围绕带电缺陷）对能带偏移校正的屏蔽影响，应该在脚本 ↓ formation-energies.py. 中计算介电常数。可以通过对原始 GaAs 块体晶胞进行 OpticalSpectrum 分析计算介电常数。重要的是，计算对于 k 点能够很好地收敛，且 OpticalSpectrum 中的能带数量足够对所有相关跃迁进行采样。本例中，20x20x20 的 k 点网格，费米能级以下 10 条能带和以上 20 条能带足以确保这一点。

将优化的 GaAs 块体结构拖放到 Builder 并发送到 Script Generator。添加一个 New Calculator，除 k 点外，参数设置同上文计算相同，此处 k 点必须为 20x20x20。添加一个 OpticalSpectrum 分析模块，设置 k 点为 20x20x20，费米能级以下 10 条能带和以上 20 条能带。现在的窗口应如下图所示。

运行计算，应该不会超过几分钟。您现在应该在 LabFloor 上看到 OpticalSpectrum，用右侧的 Optical Spectrum Analyzer 打开它，将会出现一个如下图片的窗口：

本例中，我们只对上边那张图的截距感兴趣，看起来大约为 13。要获得精确计算的值，请打开带有 Text Representation 的 OpticalSpectrum，找到 “Real part of the dielectric tensor” 部分。在这里，您可以找到零频率下介电常数的计算值为 12.9，如下图所示：

• 英文原文：https://docs.quantumatk.com/tutorials/charged_defects/charged_defects.html#charged-defects
• [EMM05] Fedwa El-Mellouhi and Normand Mousseau. Self-vacancies in gallium arsenide: An ab initio calculation. Phys. Rev. B, 71(12):125207, 2005. doi:10.1103/PhysRevB.71.125207.
• [FGH+14] Christoph Freysoldt, Blazej Grabowski, Tilmann Hickel, Jörg Neugebauer, Georg Kresse, Anderson Janotti, and Chris G. Van De Walle. First-principles calculations for point defects in solids. Rev. Mod. Phys., 86(1):253–305, 2014. doi:10.1103/RevModPhys.86.253.
• [FNVanDWalle09] (1, 2, 3) Christoph Freysoldt, Jörg Neugebauer, and Chris G. Van De Walle. Fully Ab Initio finite-size corrections for charged-defect supercell calculations. Phys. Rev. Lett., 102:016402, 2009. doi:10.1103/PhysRevLett.102.016402.
• [PRC+08] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke. Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces. Phys. Rev. Lett., 100(13):136406, 2008. doi:10.1103/PhysRevLett.100.136406.