基于晶格化的动力学Monte Carlo方法研究钛酞菁配合物氧化还原反应催化能力

周 苏, 胡 哲, 高 扬, 姜 缜

(同济大学 a. 汽车学院； b. 中德学院， 上海 201804)

一般情况下，比较催化剂的性能，可以通过实验获得以下两个参数：① 单位时间内单位活性中心所转化的反应物分子数(Turnover Frequency，TOF)，即单位时间内在催化剂单位面积每个活性位置的转化数；② 活性位置密度(Site Density，SD)，即催化剂单位面积上的活性位置.但是，对于未知或者具有新型结构的催化剂，很难通过实验得到TOF值或SD值.例如，在研究一个可能的催化反应机理时，整个反应过程持续时间很短，即使超高速摄像机也只能获得若干张宏观变化的图像.若想了解其微观层级的变化并进一步推测其催化性能，并与其他催化剂进行比较，需借助其他方法或工具.

随着计算机技术的发展，计算化学已经能够在没有实验情况下给出定性或定量结论，为寻找或设计潜在的新型催化剂材料提供方向性的信息[1].例如，以钛酞菁配合物(TiPc)为研究对象，周苏等[2]基于密度泛函理论构建了TiPc氧化还原反应机理模型，通过计算其Gibbs自由能等指标，分析预测了TiPc的催化能力.动力学Monte Carlo方法(KMC)也是一种应用广泛的计算化学工具[3].例如，汤星舟[4]结合动态Monte Carlo法与分子动力学研究了C面蓝宝石衬底的氮化现象，分析总结了表面氮化的相关机理；Franziska Hess等[5]利用KMC研究了催化剂RuO2对于CO和HCl氧化催化反应的作用.

本文以TiPc为研究对象，基于晶格化的KMC并借助Kmos软件，通过求解Markov链主方程，计算TiPc的TOF值并以其作为催化活性的参考评价指标，研究反应温度和氢氧气体分压比对TiPc催化活性的影响；通过比较未知TiPc和已知Pt两者的TOF值，预测未知催化剂TiPc的催化活性；另一方面，从计算化学工具的角度，以TiPc作为对象，旨在验证KMC用于研究未知催化剂催化活性的可行性.

1 模型与方法

1.1 计算软件

图1 Kmos软件GUI界面(氢氧反应图)Fig.1 GUI interface of Kmos software (hydrogen and oxygen reaction)

Kmos是基于Python的一款开源软件，其模型建立是由用户通过Python编写py文件作为输入，然后由Kmos模拟给出结果.为了展示催化剂表面的微观变化，该软件也支持UI界面，图1所示为氢气和氧气在催化剂表面某一反应时刻的截图，图中深色点代表氧原子，浅色点代表氢原子，白色点代表空位，基底圆形代表催化剂.

1.2 基于晶格化的KMC及应用的假设条件

基于晶格化的KMC将整个催化剂薄膜简化为10×10个晶胞所组成的平面，同时以单个晶胞作为研究对象，以基本反应的动力学参数作为输入，在晶胞表面采用KMC求解Markov链主方程[5-6,7,9]及微观动力学模型；根据给定的反应条件(如温度和所有反应物的分压)，可以计算或获取催化剂的TOF值、晶胞表面的反应物或/和生成物成分、反应链中的步骤、反应机理等信息[10-11].针对所研究的对象，应用晶格化的KMC时，通常给出如下4个方面的假设.

1.2.1特殊位置的吸附和晶格排布 假设所有反应物和中间物都吸附在特殊位置，即晶体表面的活性位置.在实际的KMC模拟时，只考虑这一个晶格，整个晶体表面可以看成由这个晶格重复排列组合而成.例如，低指数单晶表面中一个晶格是由多个相同的晶胞组成.图2所示为一个面心立方排列(100)的催化剂模型，它有2个活性位置[5].

图2 催化剂表面CO氧化反应(包含吸附、扩散、反应和解吸附过程)模型示意图Fig.2 Catalyst surface diagram of CO oxidation reaction on catalyst surface (including adsorption, diffusion, reaction and desorption process)

1.2.2Rare-event动力学和Markov链主方程 假设时间演化是基于Rare-event动力学的.利用时间尺度分离，根据微观动力学理论，假设事件的发生是相互独立的，可以基于Markov近似来推演反应过程.因此，事件(从一种组态到另一种组态)的发生可以通过Markov链主方程描述：

(1)

1.2.3基本反应的确定 假设前文提到的一个晶格中组态变化只考虑1～10个位置.因此，可以根据局部产物Ea的晶格状态和变化后的局部产物Pa的晶格状态以及相应的伴随速率Ka来确定一个基本反应，即

(2)

1.2.4状态转移矩阵的大小和结构 式(1)中w可分解成一系列基本反应：

(3)

(4)

式中：ka为基本反应a的反应速率.

以上假设用于简化Markov链主方程的求解.

1.3 模型

1.3.1催化反应速率计算 KMC模拟的主要特点之一就是可以计算各种催化剂的反应速率，包括TOF值.如果基本反应只有一个最终生成物，那么TOF的计算只考虑这一种生成物.

对于基本反应a，在任意时刻t，

TOFa=

(5)

1.3.2参数计算 基于密度泛函理论(DFT)，利用Amsterdan Density Functional(ADF)等软件，根据构造的TiPc在氧化还原过程中形成的各粒子能量变化，可计算TiPc的氧吸附能：

EbO2=EML-O2-(EO2+EML)

(6)

式中：EML-O2为氧分子和吸附粒子的总能量；EO2为孤立的氧分子的总能量；EML为吸附粒子的能量.此式也可用于计算Pt的氧吸附能.

同理，计算TiPc的氢吸附能：

EbH2=EML-H2-(EH2+EML)

(7)

式中：EML-H2为氢分子和吸附粒子的总能量；EH2为孤立的氢分子的总能量；EML为吸附粒子的能量.此式也用于计算Pt的氢吸附能.

2 结果与讨论

2.1 TiPc的TOF值

在模拟TiPc表面氢气和氧气反应时，反应温度T和TiPc表面的氢气分压与氧气分压的比值pH2/pO2是两个重要的参数.在Kmos软件运行环境下，取T在470～570 K(间隔5 K)之间的21个值和pH2/pO2在 0.1～10之间的40个值，可以计算得到相应的TOF值.以温度为470、520和570 K的情况为例，分析反应温度T和氢氧气体分压比值pH2/pO2对TiPc催化剂TOF值的影响.

2.1.1TOF值随反应温度T的变化 由化学反应基本原理可知，反应温度T是影响催化剂活性(TOF值)的重要因素.温度较低时，催化剂活性较小，对于反应促进效果较差；随着温度上升，催化效果逐渐增强.但是，大多数催化剂都有其催化活性适宜温度范围.温度过高，会使催化剂发生烧结，从而破坏其催化活性.为了分析反应温度对TOF值的影响，图3所示为3种温度条件下不同氢氧气体分压比值对应的TOF值(实线部分).与温度条件520 K及570 K相比，温度470 K下的TOF值过低，所以图3的虚线部分也给出了lg(TOF)值，以便于对比分析.

图3 不同温度下pH2/pO2对应的TiPc催化剂TOF及lg(TOF)值Fig.3 TOF and lg(TOF) values of TiPc under different pH2/pO2 at 470, 520, 570 K

由图3可明显看出：① 随着反应温度上升，同一气体分压比值对应的TOF值也依次增高，这符合催化剂基本特性.但反应温度超过催化剂适宜反应温度区间后，催化剂催化性能下降.Choi等[14]及陈维民等[15]通过各自的研究发现，催化剂会随着温度上升而发生聚集，超过一定程度会降低催化剂催化活性；② 与温度区间520～570 K相比，在温度区间470～520 K内TOF值的增长幅度较小，这说明TiPc的活性虽然与温度正相关，但并不是线性相关；③ lg(TOF)值变化随温度变化近似呈线性关系.

2.1.2TOF值随氢氧气体分压比值pH2/pO2的变化 除反应温度T外，氢氧气体分压比值pH2/pO2同样会影响催化剂活性(即TOF值).为了研究相关规律，将3种温度条件下TOF值(图3实线部分)在横坐标 0.1～0.33 区间内的曲线放大得到图4.

图4 不同温度下pH2/pO2在 0.1～0.33 区间时对应的TiPc催化剂TOF值Fig.4 TOF values of TiPc within pH2/pO2 0.1-0.33 at 470, 520, 570 K

图5 温度470至570 K条件下不同氧氢气体分压比值对应的TiPc催化剂lg(TOF)值Fig.5 lg(TOF) values of TiPc under different O2 and H2 partial pressure ratios between 470 K and 570 K

结合图3和4的信息可知，3种不同温度条件下TOF值的最大值出现在氢氧气体分压比值0.289～0.326区间内.在氢氧气体分压比值0.1～0.2之间，TOF值单调缓慢上升；在氢氧气体分压比值0.2～2区间内，TOF值先急速上升至最高点，然后单调急速下降；从氢氧气体分压比2开始，TOF值单调平稳下降至接近0值.催化反应主要经历以下3个步骤[12-13]：① 反应物被催化剂吸附在表面；② 反应物在催化剂表面进行反应；③ 反应生成物从催化剂表面脱附.其中，步骤①对应前文所计算吸附能，步骤③保证反应生成物能够及时富集及保证催化剂再生.根据 1.2 中的假设，所有的反应物和中间物都吸附在表层催化活性位置.因此，根据化学平衡移动原理，当氢氧分压比从零开始升高时，催化剂表面反应物氢气的浓度增加，反应倾向于向正反应方向移动，对于催化剂表面的活化中心具有很好的活化作用.但是，氢氧分压比进一步增大，使得氢气的覆盖度过大，会出现“氢中毒”现象，反而使催化剂活性下降，抑制反应的进行.本文所用方法不能从定量上解释TOF值急速变化的原因，但是，从上述的分析中可以推测，TiPc催化剂的活性区域在氢氧分压比值 0.2 ～2之间，此区间内氢气和氧气在催化剂表面以较为合适正向反应的比例共存.

2.1.3采用2D Map图分析TOF 为了能够同时研究反应氢氧气体分压比及反应温度对TiPc催化活性或反应速率的影响，利用线性插值，将上述21组数据绘制成2D Map图；为了分析方便，取氢氧气体分压比值的倒数即氧氢气体分压比值为横坐标，纵坐标为反应温度(见图5)；二维图中区域的颜色深浅则对应TOF值(反应速率)，颜色越浅代表TOF值越大.从TiPc的计算数据及上文的结论可知，在氢氧气体分压比值 0.1～0.2 区间内，TiPc的TOF值非常小，接近于零.为了增强这种区别，取TOF值的常用对数.图5所示为反应温度 470～570 K 下不同氧氢气体气体分压比值对应的 lg(TOF) 值.

图5中暗部区域对应的横坐标范围为氧氢气体分压比值 3.8～10，对应于氢氧气体分压比值 0.1～0.26，基本符合 2.1.1 和 2.1.2 得出的结论，该区域内TiPc反应活性远低于左侧亮部区域，且亮部区域和暗部区域突变位置即为TOF值陡降陡升位置.

2.2 Pt的TOF值

采用上述的相同条件，在Kmos软件运行环境下计算得到相应的Pt催化剂TOF值，用于分析反应温度T及氢氧气体分压比值pH2/pO2对TOF值的影响，从而进一步比较TiPc与Pt两者的催化活性.

2.2.1TOF值随反应温度T的变化 由Pt催化剂的计算数据可知，Pt催化剂TOF值在氢氧气体分压比值 0.8～10 区间内几乎不变，为了便于分析，只选用氢氧气体分压比值 0～0.8 区间内的计算数据.图6所示为3种温度条件下不同氢氧气体分压比值对应的TOF值(实线部分)，且与TiPc催化剂情况相同，为突出470 K与520 K温度下TOF值曲线间的差别，图6的虚线部分也给出了 lg(TOF)值.

图6 不同温度下pH2/pO2对应的Pt催化剂TOF及lg(TOF)值Fig.6 TOF and lg(TOF) values of Pt under pH2/pO2 at 470, 520, 570 K

由图6可看出，不同温度条件下，Pt催化剂表现出与TiPc催化剂相似的特性：① 随着反应温度上升，同一气体分压比值对应的TOF值也依次增高；② 虽然TOF值与温度不是线性正相关，但是 lg(TOF) 值随温度变化近似呈线性关系.

2.2.2TOF值随氢氧气体分压比值pH2/pO2的变化 为了研究氢氧气体分压比值pH2/pO2对催化剂活性(即TOF值)的影响关系，将3种温度条件下TOF值(图6实线部分)在横坐标 0.1～0.26 区间内的曲线放大得到图7.

由图7可以看出，3种不同温度条件下TOF值的最大值均出现在氢氧气体分压比值0.229处，此时Pt催化剂的催化活性最高；与TiPc催化剂相比，Pt催化剂TOF值的最大值明显小几个数量级，所以也可以预测TiPc催化剂的活性高于Pt催化剂.

图7 不同温度下pH2/pO2 在 0.1～0.26 区间时对应的Pt催化剂TOF值Fig.7 TOF values of Pt within pH2/pO2 0.1-0.26 at 470, 520, 570 K

2.2.3采用2D Map图分析TOF 将Pt催化剂21组原始数据利用线性插值绘制成2D Map图.由图6实线部分可知，在氢氧气体分压比值 0.1～0.3 区间内，Pt催化剂TOF值远大于其他区域.为了增强这种区别，取TOF值的常用对数.图8所示为反应温度470～570 K下不同氧氢气体分压比值对应的lg(TOF)值.

图8中亮部区域对应的横坐标范围为氧氢气体分压比值 3.9～10，对应氢氧气体分压比值 0.1～0.256，基本符合本节 2.2.1 和 2.2.2 所得结论，该区域内Pt催化剂反应活性远高于左侧暗部区域，且亮部区域和暗部区域突变位置即为TOF值陡降陡升位置.

图8 470～570 K条件下不同氧氢气体分压比值对应的Pt催化剂lg(TOF)值Fig.8 lg(TOF) values of Pt under different O2 and H2 partial pressure ratios between 470 K and 570 K

2.3 分析与讨论

TiPc催化剂TOF值与反应温度T及氢氧气体分压比值pH2/pO2存在一定复杂关系.TiPc催化剂TOF值会随温度变化而发生改变，随着反应温度由470 K增长到570 K，TOF值同样呈现出递增的趋势；氢氧气体分压比值同样会影响催化剂活性(TOF值)，TiPc催化剂TOF值的最大值出现在氢氧气体分压比值 0.289～0.326 区间内，氢氧气体分压比值在 0.1 至该区间，TOF值首先接近于零，之后单调递增；氢氧气体分压比值在该区间至10，TOF值单调递减直至接近0值.此外，TiPc催化剂TOF值的最大值出现在反应温度570 K和氢氧气体分压比值 0.326 处，其值为 2.034×107.

与TiPc催化剂不同，Pt催化剂TOF值在氢氧气体分压比值 0.1～0.229 之间单调递增，在氢氧气体分压比值 0.229～0.26 之间单调递减，可以推测氢氧气体分压比值 0.1～0.26 区间是Pt催化剂的活性区域；当氢氧气体分压比值处在 0.26～10区间内，TOF值接近于零且单调递减.Pt催化剂TOF值的最大值出现在反应温度570 K和氢氧气体分压比值 0.229 处，其值约为82.

与Pt催化剂相比，TiPc催化剂TOF值的最大值以及TiPc催化剂的2D Map图中亮部区域对应氢氧气体分压比值区间都明显更大，可以推测TiPc催化剂的活性区域更广且其催化活性更高.

3 结语

基于晶格化的动力学Monte Carlo方法，采用TOF值作为参考评价指标，评价TiPc催化剂在反应中的催化活性，模拟气-固界面氢、氧分子在催化剂表面的反应.讨论了TiPc在不同反应温度及不同氢氧气体分压比值的反应条件下的催化活性.TiPc催化剂TOF值随反应温度上升而递增，同时也随氢氧气体分压比值上升而首先递增之后递减.TiPc催化剂TOF值的最大值出现在反应温度为570 K和氢氧气体分压比值 0.326 处.采用与TiPc催化剂相同的反应条件，在Kmos软件运行环境下计算得到相应的Pt催化剂TOF值，并与TiPc催化剂TOF值进行比较，得到了与文献[2]相同的预测结果，即TiPc催化剂催化活性强于Pt催化剂.

因计算方法和计算能力的局限性，本文仅研究了TiPc在相对较简单的气-固界面上反应的催化性能，没有涉及质子交换膜燃料电池中催化剂在固-液界面的真实反应环境.