CN109192250A - 一种多相催化中克服表面物种快速迁移的加速模拟方法 - Google Patents

Abstract

本发明涉及一种多相催化中克服表面物种快速迁移的加速模拟方法，包括以下步骤：1)开始模拟并初始化催化剂晶体网格模型；2)将多相催化过程中的事件分为快事件和慢事件两类；3)遍历表面位点进行事件匹配，对慢事件列表创建相应的慢事件概率列表，并采用轮盘赌算法选择当前轮迭代的事件；4)对慢事件概率列表执行VSSM标准算法来进行体系推进模拟，包括选取事件以及反应位点、体系时间推进以及局部环境的重新匹配并更新事件列表；5)对快事件进行快速迁移物种的重分布，并更新表面构型；6)重复步骤4)直到达到最大迭代步数或者最大体系推进时间后进行数据后处理并终止。与现有技术相比，本发明具有实用广泛、快速推进、更具物理意义等优点。

Description

一种多相催化中克服表面物种快速迁移的加速模拟方法

技术领域

本发明涉及统计学的数值模拟计算领域，尤其是涉及一种多相催化中克服表面物种快速迁移的加速模拟方法。

背景技术

动力学蒙特卡洛(kinetic Monte Carlo)模拟在连接微观尺度第一性原理计算与宏观动力学现象中扮演了至关重要的角色，该方法目前已经被广泛而且成功的应用在了复杂表面多相催化的各种研究中。目前研究人员可以通过商业的量子化学软件计算出催化体系中各个可能的基元反应的热力学和动力学能量数据，基于此数据和动力学以及蒙特卡洛的方法可以使用计算机快速的模拟出该催化体系的各种宏观动力学性质，例如转化频率/反应速率，催化剂表面的物种中覆盖度以及表面构型等等。

动力学蒙特卡洛是一种基于统计学的数值模拟计算方法，它本身是对分子动力学模拟的一种粗粒化近似，该方法把催化剂表面抽象成“态”从而构建马尔科夫链来模拟出催化体系随时间不断演化的过程。整个动力学蒙特卡洛模拟遵从Master Equation公式，程序的具体实现有很多不同的算法其中在多相催化体系的模拟中 VSSM算法是最常用的算法之一。在模拟中我们需要对表面发生的事件进行选取并执行，选取时我们通过轮盘赌等类似的选择算法根据时间发生的反应速率进行选取，若事件发生的速率越快，被选中的概率就越大，反之亦然。在执行之后我们根据指数分布抽取随机数计算并获取体系实际时间的推进。具体计算的公式如下所示：

t_draw＝-(1/k_tot)ln(r)

若所有的事件中存在一个或者多个非常快的事件，则动力学蒙特卡洛模拟则会陷入一种效率及其底下的状态。若某一个或者多个的反应速率非常大(下文中称为快事件)，根据轮盘赌算法，算法总会以极高的概率在事件列表中选择出快事件。随后算法会在表面执行此快事件并更新表面构型，然后根据上面的公式获取体系推进的时间并对时间进行更新，由于快事件的反应速率值非常大就会导致公式中k_tot的值非常的大，进而导致时间推进步长t_draw会非常的小，通常小于10^-10s。这时， kMC模拟想要达到s级别的时间推进程度则需要十亿甚至更多的迭代次数，使得 kMC模拟几乎在有限的时间内达不到想要的宏观结果。在多相催化中，物种在催化剂表面的迁移过程通常具有很小的能垒，也就是具有非常快速的迁移速率。而迁移过程又是表面反应中几乎不可避免的过程，这就给多相催化的动力学蒙特卡洛模拟带来了极大的困难

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种多相催化中克服表面物种快速迁移的加速模拟方法。

本发明的目的可以通过以下技术方案来实现：

一种多相催化中克服表面物种快速迁移的加速模拟方法，包括以下步骤：

1)开始模拟并初始化催化剂晶体网格模型；

2)将多相催化过程中的事件分为快事件和慢事件两类；

3)遍历表面位点进行事件匹配，对慢事件列表创建相应的慢事件概率列表，并采用轮盘赌算法选择当前轮迭代的事件；

4)对慢事件概率列表执行VSSM标准算法来进行体系推进模拟，包括选取事件以及反应位点、体系时间推进以及局部环境的重新匹配并更新事件列表；

5)对快事件进行快速迁移物种的重分布，并更新表面构型；

6)重复步骤4)直到达到最大迭代步数或者最大体系推进时间后进行数据后处理并终止。

所述的步骤1)中，初始化催化剂晶体网格模型包括晶格抽象、表面物种构型初始化以及位点类型初始化。

所述的步骤2)具体包括以下步骤：

21)计算所有事件的反应速率；

22)根据计算得到的反应速率将所有的事件分为快事件和慢事件两类；

23)对应创建慢事件列表和快事件列表，所述的慢事件列表用以执行常规算法推进体系时间，所述的快事件列表用以进行快速物种的表面重分布处理；

所述的步骤21)中，反应速率r的计算式为：

其中，G_a为过程发生所需的吉布斯自由能垒，k_B为玻尔兹曼常数，h为普朗克常数，T为温度。

所述的步骤22)具体包括以下步骤：

221)手动分类：根据经验设定反应速率的阈值，超过该阈值，则判定为快事件，否则判定为慢事件；

222)自动分类：当前事件没有设定反应速率的阈值，则设定事件在模拟过程中正逆过程连续发生次数的阈值，若超过该阈值，一般设定为10，则判定为快事件，否则判定为慢事件。

所述的步骤5)中的重分布具体包括以下步骤：

51)对当前构型进行匹配，在当前构型中动态的识别出快速迁移物种与慢物种；

52)从当前表面构型中抽离所有快速迁移物种，并将其放入构建的重分布物种栈中；

53)更新抽离快速迁移物种后的局部环境对应的事件列表；

54)确定快速迁移物种重分布的空间；

55)重新撒点内循环，完成重分布。

所述的步骤55)具体包括以下步骤：

551)随机打乱所有重分布空间中的点以及重分布物种栈中物种的顺序；

552)进行重分布尝试，具体为：

5521)从重分布物种栈中取出物种；

5522)从重分布空间中随机选取一个位点；

5523)尝试重分布操作，将物种放入该位点，若被禁止，则重复执行步骤5522) 直到此物种被放入重分布空间中；

所述的禁止是对随机选中的位点附近的局部环境进行判别，若此局部环境满足该物种预先定义的吸附事件所发生的条件，则允许该物种重新放入；反之，若局部环境不满足，则被禁止；

5524)由于物种重分布导致局部环境变化，对局部环境相关的事件列表进行更新；

553)重复步骤552)直到所有重分布物种被重新放入重分布空间中。

与现有技术相比，本发明具有以下优点：

一、实用广泛：本发明是在广泛使用的VSSM算法上的改进，因此具有广泛的实用性，可以应用在二维表面甚至是三维晶格的动力学蒙特卡洛模拟中

二、快速推进：本发明采用的时间尺度解耦的思想，能够使得体系会随着慢事件的发生进行快速的推进，不会受到快事件的影响而使体系陷入极慢的推进状态。

三、更具物理意义：在时间解耦的基础上，增加了空间表面物种的重分布处理，更加真实的模拟了快速迁移物种真实的行为，使得模拟结果更具有物理意义。

四、在重分布内循环中，本发明使用了基于吸附过程的重分布方法，此方法可以将表面物种的相互作用因素考虑进来，进而能够使得重分布处理包含物种相互作用信息，同样使模拟更具有物理意义，本发明的重分布操作，本质上是一种基于吸附过程的重分布处理，将快物种去除并保存然后重新在可行空间的撒点，可以看作物种重新吸附在表面的过程，而吸附的过程，则可以考虑表面物种相互作用的影响，因此我们的重分布过程可以包含物种相互作用的影响，从而得到比完全随机分布更具有物理意义的分布结果进而可以更好的继续动力学蒙特卡洛迭代。

附图说明

图1为本发明的方法流程图。

图2为使用量子化学软件计算出的CO在Pt100表面氧化的可能基元反应路径以及能量数据。

图3为表面物种覆盖度演化曲线，其中，图(3a)为使用标准VSSM对事件的模拟结果，图(3b)为使用本发明的改进时间尺度解耦重分布方法对事件的模拟结果。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

实施例

如图1所示，本发明提供一种多相催化中克服表面物种快速迁移的加速模拟方法，包括以下步骤：

1)开始模拟并初始化催化剂晶体网格模型；

2)将多相催化过程中的事件分为快事件和慢事件两类，分类方式包含手动分类和程序自动判别:1.手动分类是在定义事件的配置文件中添加fast域，如果为 True则为快事件若为False则为慢事件；2.自动判别是通过在程序运行中判断如果一个事件的正逆过程连续发生n次则自动判别为快事件，其中n为预先设定的参数；

3)遍历表面位点进行事件匹配，对慢事件列表创建相应的慢事件概率列表，并采用轮盘赌算法选择当前轮迭代的事件；

4)对慢事件概率列表执行VSSM标准算法来进行体系推进模拟，包括选取事件以及反应位点、体系时间推进以及局部环境的重新匹配并更新事件列表；

5)对快事件进行快速迁移物种的重分布，并更新表面构型；

6)重复步骤4)直到达到最大迭代步数或者最大体系推进时间后进行数据后处理并终止。

本例中，以CO在Pt100进行表面催化氧化进行模拟实例，图2所示为CO在 Pt100可能发生的所有事件路径，由量子化学软件VASP计算获得，具体基元反应表达式以及能量数据如表1所示。以此我们可以将整个过程的事件代入到动力学蒙特卡洛流程中进行动力学模拟。

表1所有可能发生的基元反应表达式以及能量数据

下面为具体时间尺度分离重分布针对CO在Pt100表面模拟过程的详细实施过程:

1、对于Pt100表面我们可以抽象建立一个大小为60*60的网格来模拟。

2、快慢事件的分类分为手动和自动，从表1中可以轻易看出，CO的迁移反应的迁移能垒为0.17eV，能垒很低，我们可以手动将其相关的事件定义中添加fast 标签将其标记为一个快速迁移反应；若没有手动标记，程序在执行过程中会判断此迁移反应事件是否连续重复发生多次，若是，则自动标记为快事件，反之亦然。

3、程序扫描网格的所有格点，并对所有分类的慢事件进行匹配创建事件概率列表，并使用标准动力学蒙特卡洛算法进行事件选取，执行事件和体系时间推进。例如执行CO吸附事件，则是在发生反应的格点上放置一个CO物种，并更新表面构型和概率列表，根据公式计算出推进时间t_draw对体系事件进行累加。

4、在扫描格点的同时，针对每个格点的局部环境和快事件的定义进行匹配，若匹配则将该位点参与迁移的物种标记为快物种，在本例中就是网格中能够发生迁移的所有CO吸附物。程序将其从网格中取出并存储在额外的数据结构中，以便接下来的重分布处理。

5、抽取了网格中所有能够发生快速迁移的CO物种后在网格可行空间中实施基于吸附的CO物种重分布，一次性将所有CO快物种的迁移进行了处理，并可以将CO和其他物种的相互作用影响。具体实施过程：从栈中取出栈顶的CO物种，尝试对网格可行空间中的格点进行随机选取并执行CO吸附事件，若事件与格点匹配，则执行成功，更新表面构型，进而对栈中的下一个CO进行重复操作；若匹配未成功，则再次重复选区网格点进行尝试直至成功。

6、重复上述快慢事件的处理，直到体系模拟满足终止条件即可。

为了对比进行对比，我们分别使用标准VSSM和改进的Time-scale decouplingmethod对同一组事件进行模拟，得到的表面物种覆盖度演化曲线如下图3所示。

表2不同算法平均每蒙特卡洛步迭代体系时间推进程度

	Standard VSSM	Time-scale Decoupling
			Tps(s)	7.9×10<sup>-13</sup>	1.8×10<sup>-4</sup>

通过图3和表2可以看到，改进的算法可以在更少的迭代次数下获取与标准算法相同的稳态结果，同时改进的算法可以使得体系时间推进的更快(快8～9数量级)。

Claims (7)

1.一种多相催化中克服表面物种快速迁移的加速模拟方法，其特征在于，包括以下步骤：

1)开始模拟并初始化催化剂晶体网格模型；

2)将多相催化过程中的事件分为快事件和慢事件两类；

3)遍历表面位点进行事件匹配，对慢事件列表创建相应的慢事件概率列表，并采用轮盘赌算法选择当前轮迭代的事件；

4)对慢事件概率列表执行VSSM标准算法来进行体系推进模拟，包括选取事件以及反应位点、体系时间推进以及局部环境的重新匹配并更新事件列表；

5)对快事件进行快速迁移物种的重分布，并更新表面构型；

6)重复步骤4)直到达到最大迭代步数或者最大体系推进时间后进行数据后处理并终止。

2.根据权利要求1所述的一种多相催化中克服表面物种快速迁移的加速模拟方法，其特征在于，所述的步骤1)中，初始化催化剂晶体网格模型包括晶格抽象、表面物种构型初始化以及位点类型初始化。

3.根据权利要求1所述的一种多相催化中克服表面物种快速迁移的加速模拟方法，其特征在于，所述的步骤2)具体包括以下步骤：

21)计算所有事件的反应速率；

22)根据计算得到的反应速率将所有的事件分为快事件和慢事件两类；

23)对应创建慢事件列表和快事件列表，所述的慢事件列表用以执行常规算法推进体系时间，所述的快事件列表用以进行快速物种的表面重分布处理。

4.根据权利要求3所述的一种多相催化中克服表面物种快速迁移的加速模拟方法，其特征在于，所述的步骤21)中，反应速率r的计算式为：

其中，G_a为过程发生所需的吉布斯自由能垒，k_B为玻尔兹曼常数，h为普朗克常数，T为温度。

5.根据权利要求3所述的一种多相催化中克服表面物种快速迁移的加速模拟方法，其特征在于，所述的步骤22)具体包括以下步骤：

221)手动分类：根据经验设定反应速率的阈值，超过该阈值，则判定为快事件，否则判定为慢事件；

222)自动分类：当前事件没有设定反应速率的阈值，则设定事件在模拟过程中正逆过程连续发生次数的阈值，若超过该阈值，则判定为快事件，否则判定为慢事件。

6.根据权利要求1所述的一种多相催化中克服表面物种快速迁移的加速模拟方法，其特征在于，所述的步骤5)中的重分布具体包括以下步骤：

51)对当前构型进行匹配，在当前构型中动态的识别出快速迁移物种与慢物种；

52)从当前表面构型中抽离所有快速迁移物种，并将其放入构建的重分布物种栈中；

53)更新抽离快速迁移物种后的局部环境对应的事件列表；

54)确定快速迁移物种重分布的空间；

55)重新撒点内循环，完成重分布。

7.根据权利要求6所述的一种多相催化中克服表面物种快速迁移的加速模拟方法，其特征在于，所述的步骤55)具体包括以下步骤：

551)随机打乱所有重分布空间中的点以及重分布物种栈中物种的顺序；

552)进行重分布尝试，具体为：

5521)从重分布物种栈中取出物种；

5522)从重分布空间中随机选取一个位点；

5523)尝试重分布操作，将物种放入该位点，若被禁止，则重复执行步骤5522)直到此物种被放入重分布空间中；

5524)由于物种重分布导致局部环境变化，对局部环境相关的事件列表进行更新；

553)重复步骤552)直到所有重分布物种被重新放入重分布空间中。