如何利用 DFT 计算研究催化反应

在实际的电化学反应环境中，需要建立催化模型去研究电极与电解液之间的界面环境。在高真空的环境下，实验显示吸附在 Pt(111) 表面的 H₂O 分子是平行于电极表面层的六边形的分子层，计算同时也在氧的还原反应（ORR）中发现了相似的结构类型。

在这种情况下，当氧还原反应中间物的覆盖度为 1/3ML 时，每两个中间反应物就会间隔一个 H₂O 分子，如图 1 所示。当然过于真实地模拟水环境将会增加计算的负担，增加运算时间。假如能用一个或者少量几个表面代替“真实”的表面来计算吸附自由能的话将会大大的节省计算时间，减少计算机的负担。

图 1. Pt(111) 表面层的双水分子层结构

目前为止，模拟电极中金属／溶液界面之间相互关系的理论方法主要有两个，它们分别是 Neurock 的双参考模型以及 Anderson 和 Liu 的的溶剂化计算方法。Neurock 的双参考模型主要是通过在晶胞中加入电荷以及直接的水溶液描述了极化金属／溶液界面的性质；而 Anderson 和 Liu 的溶剂化计算方法则直接通过连续介质溶剂化模型模拟了溶液中离子的分布以及溶剂对电极的作用。

如果我们仅仅通过图 1 来理解化学反应，那么我们对化学反应的理解将是非常肤浅。然而，通过图 2 和图 3 分别描述的电极 Pt(111) 表面在不同电化学势下的氧还原反应自由能图以及 Au，Pt，Ni 和 Mo 在平衡条件下的析氢反应自由能图却能让人深刻理解反应的实质。

通过查看反应自由能图中 X 轴所对应的数值变化可以使人清晰地了解在反应中有多少质子-电子对转移，反应的势垒有多高，反应路径中的最大正电荷是多少，最大过电位有多高，哪种反应中间产物更关键等，这些数据对研究催化反应都是非常重要的。

图 2. Pt(111) 表面在不同电化学势下的氧还原反应自由能图

例如，从图 2 中的反应自由能图可知，决定反应关键步骤的是 H₂O 中 OH 基的氧化过程△G_max(1.23V)=0.45eV，因此 OH 基被认为是一种毒化反应的种类，占据了表面的活性位置。

电势 U_ORR=1.23V-0.45V=0.78V 是反应自由能没有上升的步骤中的最高值。然而，在图 2 所示的氧还原反应自由能图中，由于用密度泛函理论计算 O₂ 是不精确的，所以能量为 4.92eV 的热力学值被全部应用于反应的吉布斯自由能交换，以上的方法也同样适用于分解机制。

图 3. Au、Pt、Ni 和 Mo 在平衡条件下的析氢反应自由能图

在实际的氧还原过程中，可用于建构反应自由能图的中间产物种类有很多。然而为了更加方便地查看反应自由能图，反应自由能图却没有过多描述中间产物，结果如图 2 和图 3 所示。为了保持计算模型尽可能简单的同时又能保持过电位有足够的精度，反应自由能图中所计算的反应的路径应尽量是和实际的反应路径相一致。因此，在所有构建反应自由能图的文章中，计算的反应路径都应该和实验的结果进行半定性比较，并且计算的吸附能和产生过渡态的反应能之间应该建立一种称为 Brøsted–Evans–Polanyi 的线性关系。

当我们要估计电化学反应中的过电位时，通常通过查看相关的反应自由能图。图 2 显示的是 Pt(111) 表面分别在 U＝0，U＝0.78Ｖ，U=1.23V 三种电化学势条件下的氧还原反应路径。

第一个对应的是短循环的燃料电池反应，第二个对应的是用于 DFT 计算的 U_ORR，最后一个对应的是理论上的最大电势值。对于 Pt(111) 表面在 U=1.23V 电化学势条件下的氧还原反应路径，自由能图中存在两个带有明显上升趋势质子-电子转移步骤。这种描述反应中间步骤的电势-决定步骤（PDS）展示了最大的自由能差。

至于 Pt 的氧还原反应，PDS 是 H₂O 中 OH 基的氧化和 △G_max(1.23V)=0.45eV。因此，OH 基被认为是一种毒化反应的种类，占据了表面的活性位置。电势 U_ORR=1.23V-0.45V=0.78V 是反应自由能没有上升的步骤中的最高值，且与实验的出来的结果相一致。在 U_ORR 的电势范围内，闭合的伏安特曲线中有一个最高点，该最高点描述的是在酸性条件下 OH 基在电极表面的饱和状态（图 4）。当然，过电位的评价不总是能直接和实验数据相对比的，因为实验中存在的表面毒化，扩散能垒以及动力学效应都能影响电极电流。

图 4. Pt 电极在不同电解液中的伏安特曲线

目前，反应活性趋势图已经被应用于计算催化活性的研究中。例如，Abid-Pedersen 等通过 DFT 计算建立了含氢化合物分子在不同过渡金属表面的吸附能之间的近似线性关系，该近似线性关系能够展示反应活性趋势变化，结果如图 5 所示。图 5 显示了分子 AH_x的吸附能 △E_AHx与原子 A 的吸附能之间的线性关系。

图 5. 中间产物 CH_x、NH_x、OH、SH 在台阶式表面（红色）和密堆积表面（黑色）的吸附能与 C、N、O、S 在对应表面的吸附能的关系

当斜率 y 越大时反映的是分子 AH_x 的吸附能 △E_AHx 与原子 A 的吸附能之间越密切。这也说明了原子 A 能够在很大程度上影响 AH_x 分子在不同过渡元素表面的反应活性趋势。

d-band 模型及等效介质理论参数为分析反应活性趋势提供了依据。d-band 在这些关系中扮演了桥梁作用，吸附能的变化反映了不同金属的 d-band 特征的变化规律。当吸附在过渡金属表面的 AH_x 的吸附能越大时，过渡金属的 d-band 中心越靠近费米能级；反之，当吸附在过渡金属表面的 AH_x 的吸附能越小时，过渡金属的 d-band 中心越远离费米能级。

催化活性的趋势也可以通过图形来描绘，而在图形中催化活性可以依靠合适的描述符来描述。最好的催化活性与描述符的最佳值相对应，所以火山形状图能够直观地洞察反应并且揭示不同催化剂的活性趋势。图 6 显示了析氢反应中实验测量的交换电流密度与氢的吸附能之间的关系。

图 6. 过渡金属的交换电流密度与氢的吸附自由能的关系

从图中可以发现，不同过渡金属表面的氢的吸附能所对应的交换电流密度在图中呈现出火山状分布。仔细观察后可以会发现，析氢反应过程中具有较好催化活性的过渡金属大部分集中在火山状图形的中间，而析氢反应过程中催化活性较差的过渡金属主要集中在火山状图形的两边，如 Au，Ag，Nb 等。这主要是由于以氢的吸附自由能作为析氢反应过程中催化活性的描述符，当氢的吸附自由能趋近于 0 时，析氢反应的催化效果往往较好。

可以用以下方式理解火山形状图的两部分：

（1）对于火山图左边部分来说，反应速率随着氢的吸附能降低而减小，这主要是由于电极表面缺乏可利用的 H 与 H 结合的位置，导致氢键太强；

（2）对于火山图右边部分来说，反应速率随着氢的吸附能 △E_H 增加而减小，这主要是由于随着氢在电极表面变得越来越不稳定从而使得氢离子吸附越来越困难导致氢键太弱。而 Pt 正好处在最优位置附近，因此在析氢反应过程中，每一步的反应速率都是接近最优的。

图 6 中显示的火山形状图以氢的吸附能 △E_H 或者氢的吸附自由能 △G_H*为析氢反应描述符，能够有效地区分了那些使用系统性的数据库来研究以本征的 H-金属之间反应能作为反应描述符的析氢反应，并且通过该图可以清晰地观察不同过渡金属表面的析氢反应趋势。

线性关系与反应描述符的综合利用为快速寻找高效催化剂提供了可能，因为火山形状图的顶部只能容纳少量的相关的结构，而计算少量的吸附能就意味着可以快速鉴别和去除低效的催化剂。因此，在筛选高效催化剂之前应该考虑下是否在火山顶部存在着电化学催化剂。

一种高效的催化剂应该处在火山状图的顶部，因为吸附能正好在这个位置上处在中等强度，吸附键既不是太强也不是太弱，这就是所谓的 Sabatier 原则。如果一种催化剂能够强烈地吸附反应中间物，这样虽然有利于其活化反应物，但是它的脱附能力将会变得非常困难，不能活化的气体快速脱离催化剂表面，造成催化性能的下降。如果在相反的情况下，当一种催化剂不能有效地吸附反应中间物，这也就意味着其不能有效活化反应物，催化性能自然就下降了。

Ni 元素在图 6 所示的火山图中处于顶部的左边，这也就意味着 Ｈ 在 Ni 表面的吸附能力很强，脱附能力很弱。为了提高 Ni 表面的析氢性能，应该通过合金化等方式降低 H 的吸附能，使其吸附自由能更靠近 ０。由于 Ni 在火山图中处于远离顶部的位置，因此它的催化性能在理论上还有较大的提升空间。

当然，在筛选高效催化剂之前也应该考虑下是否有零过电位的电化学催化剂。图 7 显示了在过渡金属表面发生的氧还原反应，它的催化活性存在着基本上限，即：OH 与 OOH 吸附能之差被限定在 0.4eV 以内，换句话说，OH 与 OOH 吸附能不能完全独立地变化。理想的氧还原反应催化剂要求 OH 与 OOH 吸附能之差为两倍的平衡势能，即：2.46eV。在这种情况下，△G_i(1.23V)=0，即：理想催化剂的活性应该位于图 7 中纵坐标的 0 位置，但是这种情况并没有发生在任何一种过渡金属表面。

图 7. 氧还原活性趋势与氧结合能的关系

密度泛函理论（DFT）被广泛应用于评价催化剂的稳定性，但是以往大部分的催化剂稳定性的理论设计研究都依赖较为单一的性质描述符（如偏析能），难以对高效稳定的催化剂进行有效设计，结果是理论计算预期的催化剂的稳定性并不能满足要求，这是迄今许多高效催化剂未取得实际应用的原因之一。目前对氧还原反应（ORR）催化材料的研究中，欧利辉等人提出了以合金形成能和过渡金属表面偏析能作为筛选 ORR 电催化材料稳定性的描述符。通过这些描述符能够有效地筛选出具有较高催化活性和稳定性的 ORR 电催化材料。

图 8. 用于表面偏析计算的薄板模型 

用 DFT 计算二元合金表面偏析能的方法有很多，其中一种计算方法如图 8 所示。在原始 Pt₃M(111) 超晶胞薄板中，每一层都含 4 个原子，其中 3 个为 Pt 原子，一个为 M 原子。如果通过将第一层的 M 原子与第二层的 Pt 原子进行位置交换就可以得到表面偏析的纯 Pt 表面结构（红色箭头换位方向），交换所引起的薄板总能量改变即为 Pt 的表面偏析能。如果通过将第一层的 Pt 原子与第二层的 M 原子进行位置交换可以得到表面偏析的富 M 表面结构（蓝色箭头换位方向），交换所引起的薄板总能量改变即为 M 的表面偏析能。

图 9. 用于表面偏析计算的薄板模型 

如果合金不是 LI2 结构，那么也可以选择另外一种方法来计算 M 的表面偏析能（Esegr）。以 PdM(111) 合金体系为例，M 的表面偏析能可以定义为 M 位于不同表面层的体系总能与 M 位于体相的体系总能之差，其计算模型如图 9 所示，在图中，位于表面第一层和第二层的表面偏析能可以通过以下公式来计算:

式中，E_{PdM(M,x-layer)} 代表 M 位于第一或第二层的 PdM 合金体系总能，E_{PdM(M,3rd-layer)} 代表 M 位于第三表面层（体相）的 PdM 合金体系总能。