计算材料学案例--浅谈缺陷形成能计算

       在半导体材料中掺入杂质可以有效调控体系的输运性质，此时体系的缺陷形成能主要有来自不同原子化学势的贡献和不同价态的电子化学势的贡献。缺陷形成能（formation energy）的大小决定了杂质掺杂的难易程度，不同价态的转换能级（transition levels）决定了缺陷的电离能（通常浅能级对电导的贡献要大于深能级）。详细理论见综述：Rev. Mod. Phys. 86，253（2014）；计算案例见综述：Computational Materials Science 130，1 (2017) 。

      （1）缺陷形成能计算

        主要思路和之前吸附、合金体系类似，原子数目、种类不同，要考虑原子化学势的贡献，价态不同（电子数目有变化），要考虑电子化学势的贡献。具体公式如下：

       其中D代表掺杂元素，q对应价态（得失电子数目）。在等式的右边，第一项是带电缺陷体系的总能（不同价态是通过VASP的电子总数NELECT来调整），第二项是掺杂前体结构总能；第三项是原子的化学势贡献，相对完整的体结构，n_i是原子的数目差别，n_i小于0对应的是增加原子，u_i是原子化学势；第四项是电子数目不同的贡献，其中q大于零，说明体系电子数变少，小于零对应电子增加；第5项是带电体系的修正。（这里n_i和q的符号是一致的，类似失去电子，就是正号，得到就是负号。）

     （2）中性体系情况（考虑原子化学势的约束）

         在文献Phys. Rev. Mater. 3, 074604 (2019)中，作者选择80个原子的In2O3超胞，计算了体系可能的缺陷形成能。我们这里先看中性氧空位的缺陷形成能（红色五角星位置，大约为1eV和4eV，分别对应O-poor和O-rich条件）。我们分别计算得到In 体结构 -2.56 eV， O2分子总能 -9.87 eV，体结构In32O48 总能为 -448.889 eV，出现一个氧空位的结构 In32O47 总能是-440.492 eV。

       在O-poor条件下，In的化学势是 -2.56 eV（就是金属In体结构能量，使得O的化学势最低），O的化学势是(-448.889 +32*2.56 )/48 = -7.6452 eV。根据公式计算氧空位形成能，可得：-440.492 -7.6452+448.889 = 0.7518 eV。

       在O-rich条件下，O 的化学势是 -4.935 eV（就是O2分子的O平均能量，使得In化学势最低）, In的化学势是 (-448.889 +48*4.935 )/32 = -6.6253 eV。对应的氧空位形成能为：-440.492 -4.935 +448.889 = 3.4620 eV。

        因此，O空位的形成能将在0.7518eV到3.462eV之间变化，和文献定性一致。

     （3）带电体系（transition levels）

        在文献Computational Materials Science 130，1 (2017) 中，不同价态的Si空位稳定性将随电子化学势变化（下图左边的蓝色线）。我们计算用的体结构VBM数据为5.61 eV, gap计算值为0.6eV，决定了电子的化学势范围。从文献结果看，体系将依次经历2+，1+，0，1-等价态，右边是我们计算的结果，这里q=0的价态的稳定性偏差，可能是体结构能量数值的影响。注意Si的体结构能量会影响q=0价态的位置，文献中单质的化学势用了修正，拟合不同的化合物形成能 Phys. Rev. B 85, 115104 (2012)，这个情况在氧化物中更明显。

      这里可以看到transition levels，比如电子化学势从0开始，Si空位的2+是最稳定的，在化学势达到0.1eV附近，1+价态就更稳定。如下图所示，我们在512个Si原子体系中掺入一个Al原子，并计算不同价态的形成能。可以看到，q=0和q=-1是可能稳定的，在0.1eV附近，出现q=0到q=-1的电离过程，体系从0价态逐渐变成-1价态，约束住一个电子，然后产生空穴，对应p型导电。实验的电离能是57meV，定性也是合理的。

     （4）可能的修正方法

       需要注意的是，实验中的缺陷浓度通常是非常低的，而计算中无法采用包含太多原子的超胞。尤其在计算带电体系时，缺陷形成能随体积变化收敛缓慢，通常需要考虑修正（详见文献Phys. Rev. Lett. 102, 016402 ，2009）。下图是2011年的文献，通过马德龙能的拟合，可以在原子数目较少时得到较好的结果。