计算材料学案例--表面吸附

       我们知道，在元素种类、原子数目相同情况时，通过第一性原理计算，通常认为总能最低的结构是最稳定的。在研究表面吸附时，体系包含原子数目不同，需要根据形成能随化学势的变化来判断。结合理想气体模型，我们可以计算不同随温度、压强下的化学势，得到稳定结构随温度、压强的演化。

      （1）形成能随化学势的变化

       以乙烷、乙烯、乙炔为例，我们得到的总能，并计算形成能。其中\Delta mu_H 是氢原子的化学势（相对于氢分子），这里选择在[-1.5,-0.5] eV之间。随着H化学势的增加，氢原子数目多的体系稳定性逐步增加。这里形成能有两个交点，分别在 -1.2和 -0.9 左右，当化学势低于 -1.2 eV，乙炔是比另外两个分子稳定；当化学势高于 -0.9 eV，乙烷最稳定；当化学势在 -1.2 和 -0.9 之间，乙烯最稳定。

       在文献Phys. Rev. B.33.7683,1986给出的实验中，Cs在石墨表面有不同浓度的有序相，分别对应根号7x根号7，2x2，根号3x根号3（高温亚稳相）。根据实验选择对应的超胞搭建吸附结构，分别对应C14Cs,C8Cs,C6Cs，总能分别为：-130.215 eV，-74.826 eV，-56.047 eV.

       其中E_0是这三个吸附的总能，n_Cs和n_C分别为Cs和C原子数目，E_Cs是Cs孤立原子的能量（-0.174 eV），E_C是石墨烯中C原子的平均能量（-18.449/2 eV），\Delta mu_Cs 是Cs原子的化学势（相对于Cs的孤立原子），这里选择在[-1.2,0] eV之间。

       需要注意，形成能小于零才有利于吸附。定性的结论：化学势大于-0.9 eV（附近）后，根号7x根号7的相将开始形成，超过-0.8 eV（附近）后，2x2的相变得更稳定，而根号3x根号3的形成能一直比另外两个相高，在基态相图上是不稳定的，和实验结果符合。这里的蓝色虚线对应的是Cs体结构的化学势（-0.678 eV），当化学势高于该数值，吸附结构也不稳定，Cs倾向于在石墨表面形成体结构。

（2）吸附有序相随温度、压强的演化

       在综述文献Surface Science Reports 23，43 (1996) 中，提到了K在Ir（001）表面吸附不同覆盖度的的实验，对应文献Surf. Sci. 152/153 (1985) 303，在不同条件下得到覆盖度从1/6到1/2的有序结构。我们选择了包含4层Ir原子的（001）表面结构，根据实验文献给出的超胞构造吸附结构，包含Ir4(无吸附)，Ir32K,Ir20K,Ir16K,Ir24K2,Ir8K，分别对应总能为：-32.751 eV, -264.931 eV, -166.551 eV, -133.670 eV, -201.339 eV，-67.595 eV。

       其中E_0是这三个吸附的总能，n_K和n_Ir分别为K和Ir原子数目，E_K是K原子在BCC结构中的平均能量（-1.04 eV），E_Ir是(001)表面上Ir原子的平均能量（-8.187 eV），\Delta mu_K 是K原子的化学势（相对于K的体结构），这里选择在[-2,0] eV之间。可以看到，随着化学势的增加，稳定吸附结构对应的覆盖度从0到1/6,1/5,1/4,1/3,1/2，和实验结果符合。

       在文献PRL 99, 066101 (2007)中，考虑O在金属表面的吸附，通过形成能确定稳定结构，进一步通过化学势随温度、压强的变化，得到吸附相图。在文献Computer Physics Communications 169， 20 (2005)中，K的化学势由理想气体模型得到（如下图左所示）。

        根据计算的覆盖度和化学势的关系，可以得到有序相随温度、压强的演化。温度越高，压强越低，覆盖度越低；反之，覆盖度越高。温度相图，覆盖度随压强增加而增大；压强相同，覆盖度随温度升高而减小。

（3）高通量计算方法对吸附体系的系统研究

      在文献Surface Science 171， L379 (1986)中，实验给出了H在Pd（111）表面的吸附有序相，在文献PhysRevB.58.10890（1998）通过第一原理对相应结构进行了计算。

        对于一般的吸附问题，我们可以选择不同的超胞，从覆盖度为1的结构出发，通过删除不同数目的吸附原子，得到可能的候选结构。我们选择包含4层Pd原子的（111）表面，采用12倍面积的超胞（这样可以包含2，3，4，6等有序相）。

      这里可以看成是H和空位的二元合金，从覆盖度为1的结构出发，选择加入空位，删除的是H原子，原子数目缺省（对应产生从0-12个空位的所有结构），可以得到105个不等价结构。

       公式中的x是H数目除以12（最多是12个H），E(PdH)是表面上12个H的体系总能，E(Pd)是表面上没有H的体系总能。下图右边是形成能随覆盖度的变化，蓝色的点对应105个结构的形成能和浓度的关系。红色的凸包络给出了可能稳定的结构。查看凸点对应的结构，和实验中看到的高覆盖度有序相一致。