计算材料学案例--浅谈表面

       表面的相关研究是材料计算中经常遇到的问题之一。这里我们通过金属、半导体表面，简单介绍结构的搭建（http://sagar.compphys.cn/sagar，无需注册，直接使用）、表面能计算、以及相关的电子结构分析，作为初学者的入门案例。在此基础上，大家可以考虑相关的表面吸附、掺杂等，以及催化反应等模拟。

       （1）通过晶体结构产生不同晶面指数、层数的表面结构。在网页的B2模块中，依次输入晶面指数（如1，1，1），真空层厚度（单位是埃，如输入6，对应真空距离是12），层数（如1，5代表一次生成1层到5层），Prac是精度控制。在网页的初始结构输入框，建议用笛卡尔坐标，这里的精度选择1e-3到1e-5通常都可以。生成结构后要注意检查，比如我们常见的FCC的111面应该是平的，发现每层原子数目不同，那结构是不合理的，需要调整Prac。

        输入的结构通常应该是惯用晶胞，这和晶面指数是对应的。如果得到的表面结构不是最小的原胞，可以在网页的A2模块生成原胞。如果需要扩胞，在B3模块可以根据晶格矢量的倍数生成超胞，也可以在B1模块设置体系倍数生成符合条件的超胞（不重不漏）。

       （2）计算不同Ag的100、111表面。先对体系截断能、k点密度进行测试（之前专栏文章有），确定ENCUT=400，kppa=2500，对应k点8x8x8。结构优化得到 Ag惯用晶胞（4个原子）能量是 -10.857 eV，而孤立的Ag原子能量为 -0.2 eV，可知 内聚能为 -2.91eV，和实验值2.95 eV符合。简单起见，我们只计算表面结构没有优化的情况。在计算平台（需要用户登录）中，可以选择多任务计算模块（文件夹提交、输入参数kpts=8,8,1,ENCUT=400,NSW=0）：

   完成总能计算后，我们得到不同层数（n）、不同晶面指数的体系总能，通过公式拟合：

      拟合得到参数分别为Es_(100)= 0.4567eV,Es_(111)= 0.6386eV。根据100、111表面原胞计算面积得到A_(100)=4.287A^2， A_(111)=7.425A^2，得到表面能γ(100)= 0.852 J/m^2 γ(111)=0.688 J/m^2。和文献Phys. Rev. B. 33, 7983(1986)结果相符，也是我们熟知的111表面能低，比100表面稳定。

       总能计算表明，111表面结构原子平均能量随层数增加逐渐降低，稳定性逐渐增强，逐渐趋向体结构。通过电子结构分析，可以推测层数增加，体系的活性将降低。在文献Adv. Catal. 45, 71 (2000); Phys. Rev. B 89, 115114 (2014)中，通过d带中心判断相对活性。我们的pdos计算也可以看到，层数增加，d带中心明显左移，活性降低。

  （3）计算不同TiO_2的001表面的电子结构。根据量子力学模型，材料的能隙通常随尺寸变小而增大，对应quantum confinement，比如之前文献（Acs Nano. 6,8203，2012)给出的不同尺寸的石墨烯团簇的发光。我们构造不同厚度的TiO_2的001表面，也能看到类似的结果。这里要注意，二维材料结构优化时，VASP有noZ的版本，对应体系的真空层出现在z方向，优化晶格常数时保证z轴长度不变，否则有时表面结构会优化成体结构。在计算平台中，勾选fix Z axis即可。

       我们用普通的PBE计算，TiO_2的带隙在2eV左右，不同厚度的表面带隙从2.8eV、2.4eV等变化，随着厚度的增加，逐渐逼近体结构数值。

    （4）半导体的表面悬挂键处理。对于Si的体结构，如果我们构建一定厚度的表面，当表面没有其他原子时，在导带、价带之间出现金属态；如果悬挂键被H饱和（保证所有Si原子都是4配位），体系是半导体。对于ZnO、GaN等体系，如果表面不是极性面，则可以不考虑悬挂键的处理；如果是极性面，需要考虑赝H的饱和，才能使体系出现半导体，参考文献详见

Phys. Rev. Lett. 88, 066103(2002)；Phys. Rev. Lett.  95, 266104 (2005);J. Appl. Phys. 114, 224304 (2013)等。