近日，来自北京大学应用物理与技术研究中心在读博士生刘裕、陈默涵研究员和北京大学材料学院的许审镇研究员等人采用国产开源密度泛函理论软件ABACUS计算了大尺寸的（多达1620个原子）半导体表面催化反应，将前人提出的电荷外推法（Charge extrapolation method）的应用范围从之前的金属电极表面反应拓展到半导体电极表面反应，题为“A Caveat of the Charge-Extrapolation Scheme for Modeling Electrochemical Reactions on Semiconductor Surfaces: An Issue Induced by a Discontinuous Fermi Level Change”的论文已被Physical Chemistry Chemical Physics期刊接收（DOI: 10.1039/d2cp00642a）（论文链接：http://pubs.rsc.org/en/content/articlelanding/2022/cp/d2cp00642a/unauth）。

　　电荷外推法

　　电化学以及光电化学表面反应在实验里都发生在恒定电势条件下，电势在反应过程中保持恒定。但基于密度泛函理论（Density Functional Theory，简称DFT）的第一性原理计算通常是在恒定电荷条件下进行，因此在化学反应的模拟过程中，材料表面的电荷密度以及相应的电势会发生剧烈变化。这种差异实质上来源于模拟体系的尺寸效应，即使用较小的原胞会使得计算的电势在反应前后相差较大。为了在模拟中保持恒定电势，一种途径是构造足够大的超胞减弱尺寸效应的影响；另一种途径是在电子自洽迭代时额外嵌套一个循环以调整反应过程中体系的总电子数。然而，这两种方法都面临计算成本大的困难。

　　由Chan和Nørskov提出的电荷外推法在金属电极的原子尺度恒定电荷模型上实现了近似恒定电势条件下的能量修正，其关键假设是电极表面双电层可以使用恒定电容模型近似，即在一定范围内电势与表面电荷之间具有近似线性关系。该方法得到许多课题组的采用，其中甚至包括一些基于半导体电极的光电化学反应模拟。然而，如图1所示，电荷外推法用于半导体电极表面反应时，若模拟反应过程中费米能级出现跨带隙变化，那么电势与表面电荷之间的近似线性关系就会遭到破坏，导致恒定电容模型假设失效，使得该方法无法准确应用于半导体电极体系中的恒定电位反应能计算。

　　图1. 金属电极与半导体电极表面还原反应模拟的 (a) 费米能级变化示意图，CB代表Conduction Bands（导带），VB代表Valence Bands（价带），k是费米波矢，E(k)是k点的能量；(b) 表面电荷Q与电势U（DFT计算中电势即为功函数）的变化示意图。金属电极表面电化学反应中Q与U呈近似线性关系，即电容近似恒定。由于能隙的存在，半导体电极表面电化学反应中Q与U的线性关系会遭到破坏

　　电荷外推法的推广

　　为了验证电荷外推法在应用于半导体电极表面反应时存在的问题并给出相应的修正方法，如图2所示，我们采用面心立方结构Pt最稳定的(111)面以及金刚石结构GaP最稳定的(110)面作为电极分别构造了1个和3个表面模型反应。采用电荷外推法，我们可以预测近似恒定电势反应能（E_chg-ex）。为了评估E_chg-ex的精度以及电荷外推法的有效性，我们通过扩胞的方式构造了足够大的超胞使得反应前后电势变化可以忽略不计，则可以直接得到近似恒定电势反应能（E_ref）。通过比较E_chg-ex与E_ref，即可判断电荷外推法是否有效。 

　　基于以上验证方法，我们需要处理最大1620个原子、11214个电子的大体系，然而传统的基于平面波的第一性原理计算软件例如VASP或者Quantum Espreesso难以处理上千原子的大体系计算。ABACUS采用了空间上可以严格截断的数值原子轨道，利用其局域性来构造体系的哈密顿矩阵再加上并行化处理，使得ABACUS有能力进行几千原子的大体系高效模拟。我们利用ABACUS进行了DFT计算，使用SG15模守恒赝势和PBE泛函，并且采用了ABACUS兼具效率和精度的Double Zeta plus Polarization (DZP) 级别的数值原子轨道基组：H、C、O、P、Ga和Pt的轨道截断半径分别是8，8，7，8，9和9玻尔；我们采用了共轭梯度算法做原子结构弛豫，力的收敛标准是0.04 eV/Angstrom。我们采用了Gamma only算法（数值原子基组下针对Gamma点，即布里渊区下的（0，0，0）k点采样的加速算法）。通过以上参数设置，ABACUS可以高效地进行大体系的计算。例如，对于1620个原子的体系（加真空），采用120个E5-2690v3的CPU进程上进行并行计算，单步电子迭代的时间在80秒左右。

　　图2. (a) GaP slab和 (b) Pt slab结构俯视图和正视图。正视图中红色虚线框内的原子冻结在体相位置

　　具有金刚石结构的GaP最稳定的(110)表面对许多表面反应都具有催化能力，如图3(e)，我们发现对于半导体GaP表面的析氢反应，即在GaP表面存在一个三个水分子和一个水合氢离子组合而成的基团，其中水合氢离子的一个质子会和衬底材料的一个电子以及一个质子发生反应，生成氢分子（化学表达式为H⁺_sol + e^- + H^* → H₂, 这里H⁺_sol 表示H₃O⁺-3H₂O基团，“H^*”表示吸附在电极表面的H原子），反应前后半导体衬底由于失去一个电子，费米能级从反应前的导带底变到反应后的价带顶，电势的不连续变化导致恒定电容模型失效，电荷外推法得到的ΔE_chg-ex相较于ΔE_ref严重偏小，说明在这个半导体表面反应中电荷外推法不再适用。为了消除费米能级跨带隙所引起的误差，我们在电极表面模型中额外增加了electron-donating基团，可以避免反应前后费米能级跨带隙行为的发生，如图2(f)，费米能级在反应前后始终位于导带底，由此电荷外推法预测的E_chg-ex与E_ref比较一致。此外，如图3(g)，1-H⁺_sol的ΔE_ref与2-H⁺_sol的ΔE_chg-ex之间的差值甚至小于其与1-H⁺_sol的ΔE_chg-ex之间的差值，这说明由费米能级跨带隙引起的电荷外推法的误差甚至要比结构微小差异引起的误差还要严重。

　　图3. (a)(b)(c)(d) 半导体GaP表面析氢反应初态（IS）和末态（FS）的结构图；(e)(f) IS和FS的态密度（DOS）；(g) 电荷外推法预测的反应能（ΔE）与功函数（Φ）的拟合曲线，根据拟合曲线得到的E_chg-ex与大胞反应能Eref的比较图

　　除了图3所示的GaP表面析氢反应之外，我们还对其他半导体表面电化学反应进行了计算模拟。如图4所示，若体系计算出来的费米能级反应前后没有跨带隙（即从导带底到价带顶），则电荷外推法可以得到较为合理的反应能，反之，其反应能则会严重偏低。因此电荷外推法能否用于半导体表面电化学反应的关键在于费米能级在反应前后是否会跨带隙。根据我们提出的方法，反应前后费米能级跨带隙行为不再出现，则电荷外推法在此情况下可以重新用于半导体表面的光电化学反应。

　　此工作从恒定电容假设失效的问题出发，利用ABACUS处理了最大1620个原子的体系，通过两种方法修正恒定电荷反应能以获取近似恒定电势反应能并加以比较，揭示了电荷外推法应用于半导体电极表面反应模拟时可能引起的误差，给出了电荷外推法应用于半导体表面的适用条件，并提出了当电荷外推法不再适用时如何调整表面结构使其重新适用的解决方案，为拓展电荷外推法在（光）电催化反应模拟中的应用提供了理论基础。

　　图4. 金属Pt以及半导体GaP表面（光）电化学反应的E_chg-ex与E_ref的比较图。红色虚线框表示费米能级没有跨带隙的反应，黑色虚线框表示费米能级跨带隙的反应

　　原文链接：http://mp.weixin.qq.com/s/M0R3WdrpcPtVV68rlsNZuQ