当今85%的化学品生产和全部的现代化燃油精炼都依赖于催化过程.尽管如此, 催化的意义远不止化学工业和石油精炼, 它更被视为解决社会挑战、为可持续未来开辟道路的核心学科.催化的挑战[1]主要来源于两方面, 一是实现能源和环境相关的特定应用目标, 二是催化的方法论.从分子尺度到材料尺度理解催化, 是催化方法论所面临的挑战.它既是基础研究的主要内容, 也是克服应用挑战的关键问题; 它既是理性设计催化剂的基础, 又是新催化剂从基础研究走向工业生产的促进剂.工业催化过程的复杂性, 使得从分子尺度到材料尺度理解催化十分困难.
在1922年, Langmuir[2]提出一个可行的简化途径, “大多数高分散的催化剂的结构都非常复杂. 为了简化我们从理论上研究表面催化反应, 让我们先将注意力集中于发生在平整表面上的反应. 如果这些例子上的规律能够被很好地理解, 那么继续将这些理论拓展到多孔的例子是可能的. 总之, 我们需要将表面看成一个棋盘…”然而, Langmuir脑中的“表面科学”在那个年代的实验上并不能实现.直到上世纪60年代, 随着超高真空技术的发展和各种表面敏感的物理方法的发展, 利用表面科学手段, 研究催化问题才成为可能[3, 4].在低温超高真空下, 通过一系列的表面表征技术可以获得理想状态下的催化剂在分子层次的详细信息[5~7].理论计算可以方便地获得表面吸附的平衡结构、吸附能、反应路径以及活化能等[8].表面表征技术和理论计算的结合, 使得人们能够在原子、分子层次认识理解催化过程, 从而也使得催化剂的理性设计成为可能[9].
近年来, 原位实验技术的巨大发展, 提供了大量原位条件下催化剂状态的信息[10, 11].这些实验表明, 超高真空、单晶体系和实际催化体系之间存在显著的差异[12~16].不同于低温超真空单晶表面, 原位条件下真实催化剂存在更多更复杂的影响因素[13], 如(1)催化剂表面的组成和结构可能发生改变; (2)催化剂表面多种反应位点共存; (3)催化剂表面覆盖度高, 物种间相互作用影响显著; (4)对于液-固催化体系(如, 电催化体系通常在水溶液中进行)还存在显著的溶剂化效应.因此, 原位实验的进展, 迫切需要相匹配的理论模拟的协同发展[12].
为了实现原位条件下理论模拟催化过程, 相关的理论方法近年来得到较快的发展.本文将简要综述这些理论方法的一些概念, 以及应用于多相催化研究的新进展和挑战, 主要包括: (1)限制性的第一性原理热力学(AITD, ab initio thermodynamic)分析方法; (2)偏倚的第一性原理分子动力学(AIMD, ab initio molecular dynamic); (3)微观反应动力学(microkinetics).其中, 限制性的AITD用于近似地获得原位条件下催化剂的结构和组成; 偏倚的AIMD主要应用于准确地描述势能面平缓(高温或者考虑液/固界面溶剂化效应)的催化反应过程; 最终, 需要微观反应动力学将(各种方法得到的)微观基元反应组成的复杂反应网络与催化剂的宏观性能有效地关联起来.通过对催化反应的动力学过程的研究, 指认控制催化性能的关键因素, 再针对关键因素进行定向改进.这是目前改进和设计催化剂的一条基本思路.
然而, 当前常用的微观动力学方法往往不能兼顾精度和效率.本文着重介绍一种微观反应动力学新方法——拓展的唯象动力学(XPK, extended phenomenological kinetics). XPK方法将不同时间尺度的事件分开处理, 极大地提高了效率.同时, 建立不同时间尺度之间的精确映射, 保证了精度.通过模型体系和复杂的真实催化体系, XPK的精度和效率已得到了验证. XPK方法在动力学方法上的进步, 将促进复杂多相催化体系的精确动力学研究.
最后值得注意的是, 本文主要关注于低温超高真空单晶表面的理论模拟和原位条件下真实催化剂理论模拟的差异, 因此并未讨论第一性原理的电子结构计算, 特别是密度泛函理论(DFT).事实上, 上述的三个方法, 目前在具体实现时均主要依赖于DFT理论.发展更精确、更高效、适用于更复杂体系的电子结构计算方法, 同样是多相催化理论研究的重要研究方向.此外, 鉴于原位条件下真实催化体系的复杂性, 全局优化算法和机器学习也将成为催化研究中非常重要的工具.
有限压强(p)和温度(T)对催化剂组成和结构的影响, 可以通过限制性的AITD分析方法引入[17-20].在AITD中, 特定表面构型的吉布斯自由能, 被表示为一些单独组分化学势的函数.而这些单独组分化学势通常是由反应条件确定的.由此, 即可比较不同表面构型在给定条件下的表面能.其中表面能最低的构型被认为是实验上应该观测到的.更多关于限制性AITD分析方法的细节, 可以参考综述[19]和[20].限制性的AITD分析方法可以通过很小的计算量, 提供原位条件下催化剂表面近似的组成和结构.因此, 尽管这是一个近似的方法, 目前也越来越流行[21~23].
然而, 限制性的AITD分析方法, 并不考虑动力学对表面结构和组分的影响.这是方法内禀的近似, 可能带来显著的误差.此外, 在实际应用过程中, 限制性的AITD分析方法的准确性还依赖于DFT方法的精度以及挑选的表面构型.
对于处于高温的体系或者需要考虑(液/固界面)溶剂化效应的体系, 通过少数几个构型描述体系是不够的, 而需要通过系综统计描述体系.这是由于此类体系通常具有比较平缓的势能面.在研究反应过程的自由能变化时, AIMD是目前最常用的方法.然而, 化学反应的时间尺度通常远大于AIMD中基元分子运动的时间步长, 导致直接的AIMD模拟难以实现.因此, 人们发展了一系列的增强采样方法[24~26].在多相催化的理论研究应用最广泛的是添加偏倚势方法.如伞形采样[27, 28], 温度加速动力学[29~31]和多元动力学[32~34]等.值得一提的是, 近年来人们开始发展组合增强采样方法以期获得更理想的采样效率[35~37].
目前AIMD主要应用于模拟一些简单的反应, 或者考虑溶剂化效应对基元步骤的影响[38].受限于效率, 通过AIMD直接模拟整个催化过程, 目前并不是一个很好的选择.与此相反, MD可以自然地与微观动力学模拟结合.以AIMD得到的准确速率系数作为输入, 通过微观动力学可以有效地关联分子层次的信息和整个催化反应网络的宏观表现.
表面催化的反应机理通常包含大量的基元步骤, 这些基元步骤或串行或并行组成复杂的反应网络.基于DFT的电子结构计算, 可以对一个催化体系提供诸如结构、基元反应焓变和反应能垒等极其有用的信息.然而, 如何将这些分子尺度上的微观性质, 和催化体系的宏观变量(例如给定温度、压力下生成特定产物的转化率)进行关联, 还需要进行进一步的微观反应动力学模拟.通过对催化反应的动力学过程的研究, 阐明控制催化性能的关键因素.这对于催化剂的改进和理性设计十分重要.因此, 如何有效地获得完整精准的动力学过程, 是原位条件下理论模拟多相催化过程的首要目标.微观反应动力学模型结合线性比例关系[39, 40], 可以描绘出(重要吸附物种)吸附能和活性的火山型曲线, 目前被广泛应用于催化剂理性筛选和理性设计[41].除了微观反应动力学方法本身的精度, 速率系数的精度同样会显著影响最终的结果.目前, 通过理论方法获得的速率系数的精度, 既依赖于DFT的精度, 又取决于是否准确地考虑了原位条件下复杂的环境对速率系数的影响.
目前, 国内多个理论催化的研究小组, 针对原位条件下多相催化的理论模拟, 在相关方法的发展和具体体系的研究, 均已取得显著的成绩, 具有明显的特色.例如, 北京大学高毅勤发展的温度加速动力学增强采样方法[29], 已是AIMD模拟中重要的采样方案之一.清华大学李隽课题组通过构建基于化学势的热力学模型, 并结合DFT计算、AIMD及微观反应动力学模拟, 发现Au/CeO2体系在催化CO氧化时活性中心为催化过程中动态生成的单原子Au, 并提出了“动态单原子催化”的概念[42].中国科学技术大学侯中怀和罗毅课题组, 通过反应-迁移的动力学模型, 成功地揭示了尖端增强CO2电还原的机理[43].中国科学技术大学李微雪课题组, 发展了反应条件下负载纳米催化剂Ostwald熟化和分解的一般性理论, 提出了抑制烧催化剂烧结、加速催化剂再生的一般性策略和途径[44].复旦大学刘智攀课题组近期致力于随机势能面行走(SSW)全局搜索的结构预测方法的发展与应用, 该方法在催化材料的结构预测方面已显现出了很好的应用前景[45].大连化学物理所邓伟桥课题组通过实验和理论的结合, 报道了在常温常压下, Co配位的共轭微孔高聚物可以催化环氧乙烷和CO2之间的反应, 且性能优于均相的Salen-cobalt催化剂[46].基于DFT和AIMD, 东南大学王金兰课题组揭示了光诱导单层或薄层黑磷在环境条件下降解的机理[47].南开大学周震课题组, 从理论上预言Ti2CO2单层上锚定的Ti可以作为催化CO氧化的单原子催化剂, 并通过AIMD验证催化剂在给定温度下的稳定性[48].华东理工大学龚学庆课题组基于DFT计算, 通过引入氢氧根的影响, 成功地解释了表面表征实验在CeO2表面观测到的氧空位簇结构[49].厦门大学傅钢课题组, 在一系列催化体系中, 从理论上阐明了真实催化剂表面存在的配位化合物对催化剂催化性能的影响[50].华东理工大学的胡培君课题组, 在微观动力学方法的发展和基于计算的催化剂理性设计取得了十分丰富的成果[51]等等.
常见的动力学模型大体上可以分为两类(图 1): (1)基于平均场近似的唯象动力学(PK, phenomenological kinetics), 如Langmuir-Hinshelwood型模型[52], Sabatier分析[53]以及平均场的微观动力学模型[51, 54].由于实现简便、并且高效的缘故, 唯象动力学是目前应用最广泛的动力学模型; (2)基于第一性原理的严格理论[55]——KMC(kinetic Monte Carlo).在表面催化的KMC模型中, 表面的催化位点被映射成一个网格, 而表面吸附物种占据网格的格点.网格和相应位点占据的吸附物种共同构成一个网格化的表面构象.基于网格构象的KMC模拟也称为显格子KMC.显格子KMC无需引入平均场近似、热力学平衡态近似或动力学稳态近似, 考虑体系在某一个构象下的所有可能的基元步骤(如吸附、脱附、扩散与反应等), 并由此决定它的下一步的演化去向.如图 1所示, 精度的提高, 往往伴随着效率的降低, 反之亦然.精度和效率无法同时兼顾.
原位条件下的多相催化过程, 存在内禀的不均性, 无法通过平均场近似精准地描述.首先, 工业催化剂通常为负载的纳米颗粒, 同时存在多种表面活性位; 其次, 不同于超真空条件, 原位条件下催化剂表面通常都具有一定覆盖度的表面物种.这些表面物种之间的相互作用会显著地影响基元反应的热力学和动力学性质[56, 57]、表面物种的排序[58, 59]、优势路径[60]、产物选择性[61]以及催化剂设计中“火山型曲线”的形状[53].因此, 需要更复杂的统计方法, 特别是显格子KMC[62~66]模拟.
显格子KMC模型可以在原子尺度精细地描述表面催化过程, 严格地考虑催化剂表面的不均匀性[57, 64, 67].然而, 在许多情况下, 表面的基元事件之间通常存在巨大的时间尺度分离.例如, 迁移事件通常远快于反应事件.据文献[68]报道, 在金属表面, 原子或者小分子(H, C, N, O, CO, NO等)迁移能垒仅为吸附能的0.12倍.此外, 反应网络中的快平衡步骤又远快于决速步骤.这使得直接的显格子KMC在描述许多真实催化体系的动力学时效率极低.因此, 发展能兼顾精度和效率、原位模拟催化过程的新型微观动力学方法, 是目前催化研究方法论发展所面临的挑战性关键问题, 近年来也受到越来越广泛的关注[51, 57, 64].
为了克服显格子KMC效率方面的难题, 同时保持其精度, 最近我们发展了XPK方法[69].
通过表面(通常为网格化)构象x的主方程, 可以得到表面催化详细的动力学信息.该完备态的主方程可以从第一性原理推导而来[55], 是显格子KMC的基础, 同时也是推导发展XPK方法的起点.另一方面, 化学生产通常更关心表面物种个数(或覆盖度)n的演化和各产物生成的速率.相比于x, n的主方程的维度和复杂性大大降低.通过对比x和n的主方程, 我们首先准确地确定了n的主方程中的反应趋势的具体内容, 为精确的动力学演化能够基于n的主方程提供了理论基础.
对于表面物种高速迁移的体系, 通常可以认为在任一化学反应发生之前, 体系总是通过迁移充分混合达到迁移准平衡.这意味着给定n下x的分布P(x|n)可视为准定态分布, 且该分布以及相关的平均性质均仅为n的函数.例如, 对于某一单位点反应(Single-site, SS)A→B或者双位点反应(Double-site, DS)A+C→D+E, 其反应趋势RSS(n)和RDS(n)分别为:
其中, NA是n的一个元素, 表示表面物种A的个数. NA, C(x)表示构象x上(A, C)对的个数. ki|SS(x)表示在构象x上的第i个A→B的速率系数. ki|DS(x)与ki|SS(x)具有类似的含义.据此, 通过一个仅含迁移的显格子KMC, 即可得到给定n下x的精确分布及精确的反应趋势.这为精确的动力学演化能够基于n的主方程提供了操作上的可行性.
我们将精确统计得到的反应趋势改写成PK的形式, 如下式所示
其中, Ns为表面的总位点数, Nba为格子上位点的配位数.由此得到的表观速率系数kapp, 包含了基于格子构象的精确反应趋势.通过推导发现, kapp尽管原理上应该是n的函数, 但实际上只对覆盖度较大物种的变化敏感[69].即使反应网络再复杂, 表面覆盖度较大物种的种类也很少[70].据此, 我们可以通过仅在少数的几个方向上, 对kapp进行泰勒展开即可获取一定范围内精确的反应趋势, 并由此基于n的主方程高效地演化体系的动力学.尽管XPK在原理上是严格精确的(基于迁移准平衡), 在通过泰勒展开获取kapp的时候将产生一定的误差.泰勒展开的级数越高、范围越小, 则误差越小, 但是相应的计算量会增大.考虑到该部分容易实现高效并行, 精度的提高对模拟的总时间不会有显著的影响.
同显格子KMC相比, XPK方法克服了迁移事件和反应事件之间的时间尺度分离, 并且基于维度大幅度降低的n的主方程演化动力学体系.因此, 在原理上XPK的效率将远高于显格子KMC方法.值得注意的是, 目前的XPK方法着重关注于克服迁移事件和反应事件之间的时间尺度分离, 而对于一些反应网络十分复杂的多相催化剂体系, 反应事件之间的时间尺度分离同样会大幅度降低KMC的效率.因此, 我们已进一步验证了[69], 通过引入准平衡假设, 同样可以克服快反应事件和决速步之间的时间尺度分离.未来将进一步完善该部分内容.
在实现上, XPK混合使用了仅含迁移的显格子KMC和隐格子KMC, 前者精确地统计反应趋势, 而后者高效地演化覆盖度和计算速率.最终实现精度和效率的兼顾.在一个含相互作用能的两步加氢模型体系中, 通过与显格子KMC的对比, 我们验证了XPK的精度和效率[69, 71].进一步, XPK被应用于原位下真实的催化体系——催化氨分解的火山型曲线描绘.在0.1 MPa, 850 K下, 基于线性比例关系, 我们通过XPK方法分别模拟了不考虑表面吸附物种间相互作用能(Lateral interactions, LI)和包含相互作用能的火山型曲线.如图 2所示, 在无相互作用能体系XPK可以给出和显格子KMC几乎一致的速率和覆盖度.但因为不考虑表面吸附物种间相互作用, 催化剂随覆盖度过快上升而迅速失活.而在考虑相互作用能之后, 原始的显格子KMC效率很低, 很难运行; 而XPK依然可以高效地模拟真实催化剂表面、考虑覆盖度影响的动力学过程.由图 2可见, 相互作用能会显著地影响火山型曲线的形状, 特别是最重要的火山顶位置从N*吸附能为-5 eV移动到了-6 eV.如果不考虑相互作用能, 将会大大的低估火山型曲线左半边催化剂的活性, 最终可能导致在筛选和理性设计催化剂的时候忽略掉一些能更好地平衡其他因素(如价格、稳定性等)的催化剂.这些结果显示了XPK方法的精度和效率, 同时也充分说明了原位条件下的动力学模拟对催化剂理性设计的重要性.因此, 在原位条件下多相催化动力学模拟和基于计算的催化剂理性研究中, XPK有望成为强有力的工具.
实验与理论的紧密结合, 已成为催化研究的范式.通常, 基于一个活性中心模型, 利用DFT的电子结构计算, 可以得出一套催化反应的微观机理, 从而构成了催化体系的一幅静态图像.事实上, 随着催化反应的进行, 表面的微观图像不断地在改变.不同的表面吸附物种, 可以引起表面不同方式的变化.特别地, 表面吸附物种因相互作用而动态地改变着表面的微环境, 直接影响着催化反应机理.近年来, 原位及工况条件下的表面表征新技术不停涌现、并日益发展成熟, 这迫切需要相匹配的理论模拟的协同发展.这需要催化理论的研究范式, 从传统的静态图像到动态图像的转变.
不同于低温超真空单晶表面, 原位条件下真实催化剂存在更多更复杂的影响因素.为了实现原位条件下的理论模拟多相催化过程, 多种关键的方法被发展并应用于考虑这些复杂的因素对基元反应的影响.其中, 微观反应动力学可以将微观基元反应组成的复杂反应网络和催化剂的宏观性能有效地关联, 指认控制动力学过程的决定性因素, 是原位条件下理论模拟多相催化的关键科学问题.
直接基于主方程的KMC方法, 虽然提供了一条原理上精确的模拟路径, 但受限于效率, 在实际应用上举步维艰.而另一方面, 实际应用最多的是基于平均场近似的PK.然而, PK模拟的结果即使与动力学实验数据表面上吻合, 也不能保证PK模拟所基于的反应机理、决速步骤、和表面最丰物种等假定, 在定性上是正确的.针对微观反应动力学方法往往不能兼顾效率和精度的难题, 我们课题组最近发展了XPK方法. XPK方法混合使用显格子KMC和隐格子KMC, 前者精确地统计反应趋势, 而后者高效地演化覆盖度和计算速率, 最终实现精度和效率的兼顾.因此, 在原位条件下多相催化动力学模拟和基于计算的催化剂理性研究中, XPK有望成为强有力的工具.
动态图像的准确描述往往需要更复杂、更精细的统计方法, 如分子动力学、动力学蒙特卡洛等.因此, 发展相关的可以兼顾效率和精度的新方法, 使之适用于常见的重要多相催化体系, 是理论研究的一个重要方向.另一方面, 在新方法的基础上, 准确地描述重要体系的催化过程, 并从复杂的模拟结果中抽提出简洁直观的催化概念, 将十分有利于实验上理解催化过程、改进或设计催化剂.