一种对反应分子动力学催化重组机制发生次数的统计方法

技术领域

本发明涉及化学分子反应动力学技术领域，尤其涉及一种对反应分子动力学催化重组机制发生次数的统计方法。

背景技术

分子动力学是一门结合物理，数学和化学的一套广泛应用的分子模拟方法，其基本原理是由牛顿经典运动方程计算物理系统中各个原子的运动轨迹，然后采用统计的方法计算出体系的各种性质。分子动力学在纳米尺度、纳秒级别有着无可比拟的优势，是科研领域里不可或缺的研究方法。分子动力学模拟软件LAMMPS(Large-scale Atomic/MolecularMassively Parallel Simulator)支持多种物质形态下、多种系综下大规模原子分子体系模拟计算，且具有良好的扩展性，被广泛应用于分子动力学模拟工作。

反应分子动力学综合了传统分子动力学可以进行大规模体系运算和密度泛函理论可以进行化学计算的优点，AdrivanDuin和William A.Goddard III设计的Reax FF反应力场形式已经成功的应用于大量的反应分子动力学模拟的研究中，应用体系包括碳氢有机小分子体系、高分子体系、高能材料体系、金属氧化物体系以及过渡金属催化剂体系等。

高速飞行情况下，激波层和边界层内的气体温度急剧升高。使得气体分子发生能量激发、离解、电离等一系列复杂物理化学变化，即“真实气体效应”，离解后的O、N原子在表面热防护材料的催化效应作用下发生复合反应并释放热量，不同的壁面催化特性会极大地影响飞行器表面热流。因而有必要发展相应的高超声速下近壁面高温气体效应模拟的精确、高效模型，以满足未来高超声速飞行器气动力/热环境准确、快速模拟的需求。

基于ReaxFF反应力场的反应分子动力学可以有效地从原子尺度模拟气固界面催化反应过程。反应分子动力学可以对催化反应过程的碰撞、吸附、重组、解吸附等过程可以进行准确模拟，但是在研究气固界面分子重组Langmuir-Hinshelwood(L-H)机制或Eley-Rideal(E-R)机制的研究中，由于无法直接通过lammps软件获得催化反应机制发生次数的有关参数，对不同反应机制的统计和研究造成了困难。

发明内容

有鉴于此，本发明提供了一种对反应分子动力学催化重组机制发生次数的统计方法，用以研究E-R机制和L-H机制并统计两种反应机制所占比例。

本发明提供的一种对反应分子动力学催化重组机制发生次数的统计方法，包括如下步骤：

S1：在整个模拟期间，每间隔预设时长在通量边界模型中导出一次所有双原子气体原子的轨迹文件，合并保存为dump.txt文件；

S2：对dump.txt文件进行预处理，仅保留保存数据时的计算步数、原子序号、原子种类和原子坐标，保存为dump.xlsx文件，将dump.xlsx文件以数值矩阵格式导入MATLAB程序中，命名为dump；

S3：在通量边界模型中导出所有分子的所有原子的轨迹文件，对导出的所有分子的所有原子的轨迹文件进行预处理，仅保留原子种类所在行，保存为atom.txt文件，将atom.txt文件以数值矩阵格式导入MATLAB程序中，命名为atom；

S4：通过文件dump内的原子坐标判断文件dump内的原子之间是否形成分子；若是，则标记形成分子的原子的序号、坐标，并保存在lh.txt文件中，将lh.txt文件以数值矩阵格式导入MATLAB程序中，命名为lh；若否，则不进行标记和保存；

S5：遍历文件lh中形成分子的每个原子序号，通过文件atom内与文件lh中的原子序号对应的原子在不同时刻的z轴坐标，判断每个时刻下文件atom内的原子与壁面的距离是否小于距离阈值；若是，则将文件atom内原子序号和该时刻下原子坐标保存在lh2.txt文件中；若否，则不保存；待遍历文件lh中的所有原子序号后，将lh2.txt文件以数值矩阵格式导入MATLAB程序中，命名为lh2；

S6：对于文件lh2中形成分子的所有原子均进行如下判断：

在文件lh2中获取形成某一分子的原子a与壁面的距离第一次小于距离阈值时的坐标，判断获取的原子a的坐标与文件lh2中保存的形成同一分子的原子b的所有坐标间的距离是否小于断裂阈值；若是，则原子a与原子b形成分子的机制为E-R机制；若否，则在文件lh2中获取原子b与壁面的距离第一次小于距离阈值时的坐标，判断获取的原子b的坐标与文件lh2中保存的原子a的所有坐标间的距离是否小于断裂阈值，若是，则原子a与原子b形成分子的机制为E-R机制，若否，则原子a与原子b形成分子的机制为L-H机制；

S7：统计通过E-R机制形成分子的个数和所有原子序号、以及通过L-H机制形成分子的个数和所有原子序号储存在lh3.txt文件中，得到所有双原子气体分子通过E-R机制和L-H机制重组的比例。

在一种可能的实现方式中，在本发明提供的上述对反应分子动力学催化重组机制发生次数的统计方法中，步骤S4，通过文件dump内的原子坐标判断文件dump内的原子之间是否形成分子，具体包括：

通过文件dump内的原子坐标，计算文件dump内的原子之间的距离，判断文件dump内的原子之间的距离是否小于断裂阈值。

在一种可能的实现方式中，在本发明提供的上述对反应分子动力学催化重组机制发生次数的统计方法中，步骤S1中，预设时长的范围为0.1ps～5ps。

在一种可能的实现方式中，在本发明提供的上述对反应分子动力学催化重组机制发生次数的统计方法中，步骤S5中，阈值的范围为

本发明提供的上述对反应分子动力学催化重组机制发生次数的统计方法，根据通量边界模型输出的轨迹文件，使用MATLAB程序进行后处理，研究E-R机制和L-H机制并统计两种反应机制所占比例，为高超声速下使用反应分子动力学对气固界面反应机制的研究提供了可靠方法。从反应分子动力学模拟过程输出的dump文件中提取数据，通过E-R重组机制原理对不同的原子催化重组过程进行机理判定，探究不同产物随条件变化重组机制的变化规律，设计了一种判定双原子气体分子是否形成的方法，弥补了通过lammps自带的species.out命令无法对通量边界模型进行准确的种类统计的缺陷，为反应分子动力学模拟探索各种工况下气固界面反应机制、明确产物变化的规律提供了理论基础。并且，由于在模拟计算过程中分别对不同步骤的计算结果进行独立输出，因此，输出结果亦可用于特定时刻原子/分子种类和数量的统计中。

附图说明

图1为本发明提供的一种对反应分子动力学催化重组机制发生次数的统计方法的流程图；

图2为本发明实施例1中气体内的原子信息文件dump.xlsx文件的数据格式；

图3为本发明实施例1中预处理后的原子轨迹文件atom.txt文件的数据格式；

图4为本发明实施例1中形成O

图5为本发明实施例1中形成O

图6为本发明实施例1中结果文件lh3.txt文件的数据格式；

图7为本发明实施例1中两种催化重组机制E-R机制和L-H机制的比例图。

具体实施方式

下面将结合本发明实施方式中的附图，对本发明实施方式中的技术方案进行清楚、完整的描述，显然，所描述的实施方式仅仅是作为例示，并非用于限制本发明。

本发明提供的一种对反应分子动力学催化重组机制发生次数的统计方法，如图1所示，包括如下步骤：

S1：在整个模拟期间，每间隔预设时长在通量边界模型中导出一次所有双原子气体原子的轨迹文件，合并保存为dump.txt文件；

S2：对dump.txt文件进行预处理，仅保留保存数据时的计算步数、原子序号、原子种类和原子坐标，保存为dump.xlsx文件，将dump.xlsx文件以数值矩阵格式导入MATLAB程序中，命名为dump；

S3：在通量边界模型中导出所有分子的所有原子的轨迹文件，对导出的所有分子的所有原子的轨迹文件进行预处理，仅保留原子种类所在行，保存为atom.txt文件，将atom.txt文件以数值矩阵格式导入MATLAB程序中，命名为atom；

S4：通过文件dump内的原子坐标判断文件dump内的原子之间是否形成分子；若是，则标记形成分子的原子的序号、坐标，并保存在lh.txt文件中，将lh.txt文件以数值矩阵格式导入MATLAB程序中，命名为lh；若否，则不进行标记和保存；

S5：遍历文件lh中形成分子的每个原子序号，通过文件atom内与文件lh中的原子序号对应的原子在不同时刻的z轴坐标，判断每个时刻下文件atom内的原子与壁面的距离是否小于距离阈值；若是，则将文件atom内原子序号和该时刻下原子坐标保存在lh2.txt文件中；若否，则不保存；待遍历文件lh中的所有原子序号后，将lh2.txt文件以数值矩阵格式导入MATLAB程序中，命名为lh2；

S6：对于文件lh2中形成分子的所有原子均进行如下判断：

在文件lh2中获取形成某一分子的原子a与壁面的距离第一次小于距离阈值时的坐标，判断获取的原子a的坐标与文件lh2中保存的形成同一分子的原子b的所有坐标间的距离是否小于断裂阈值；若是，则原子a与原子b形成分子的机制为E-R机制；若否，则在文件lh2中获取原子b与壁面的距离第一次小于距离阈值时的坐标，判断获取的原子b的坐标与文件lh2中保存的原子a的所有坐标间的距离是否小于断裂阈值，若是，则原子a与原子b形成分子的机制为E-R机制，若否，则原子a与原子b形成分子的机制为L-H机制；

S7：统计通过E-R机制形成分子的个数和所有原子序号、以及通过L-H机制形成分子的个数和所有原子序号储存在lh3.txt文件中，得到所有双原子气体分子通过E-R机制和L-H机制重组的比例。

在具体实施时，在执行本发明提供的上述对反应分子动力学催化重组机制发生次数的统计方法中的步骤S4，通过文件dump内的原子坐标判断文件dump内的原子之间是否形成分子时，具体可以通过以下方式来实现：通过文件dump内的原子坐标，计算文件dump内的原子之间的距离，判断文件dump内的原子之间的距离是否小于断裂阈值(即原子间键长的1.2倍)。

在具体实施时，在本发明提供的上述对反应分子动力学催化重组机制发生次数的统计方法中的步骤S1中，预设时长的范围可以控制在0.1ps～5ps范围，具体可以根据不同使用者对于通量边界模型每次删除气柱内分子的时间间隔来进行设定。

在具体实施时，在本发明提供的上述对反应分子动力学催化重组机制发生次数的统计方法中的步骤S5中，阈值的范围随着不同原子间键长而改变，大小为原子间键长，阈值的范围可以控制在

下面通过一个具体的实施例对本发明提供的上述对反应分子动力学催化重组机制发生次数的统计方法的具体实施进行详细说明。

实施例1：

以O原子垂直壁面入射为例，分别在七种工况下对通量边界模型计算过程气固界面反应机制进行研究，七种工况如下：入射动能分别为0.05eV、0.1eV、0.5eV、1eV、2eV、3eV和4eV，壁面温度为1000K。

第一步，对数据预处理，并获取形成分子的对应原子的序号。

(1)针对入射动能0.05eV，壁面温度1000K的模拟过程，每隔1ps在通量边界模型中导出一次所有气体原子的轨迹文件，合并保存为dump.txt文件。

(2)对dump.txt文件进行预处理，仅保留保存数据时的计算步数、原子序号、原子种类和原子坐标，保存为dump.xlsx文件，文件格式如图2所示，将dump.xlsx文件以数值矩阵格式导入MATLAB程序中，命名为dump，文件dump用以判断原子是否结合成为分子。

(3)在通量边界模型中导出所有分子的所有原子的轨迹文件，对导出的所有分子的所有原子的轨迹文件进行预处理，仅保留原子种类所在行，保存为atom.txt文件，文件格式如图3所示，将atom.txt文件以数值矩阵格式导入MATLAB程序中，命名为atom，文件atom用于获取原子在模拟期间的轨迹。

(4)设置断裂阈值

第二步，截取原子a与原子b在壁面附近的轨迹。

遍历文件lh中形成O

第三步，判断所有气体分子中原子重组成为分子时的反应机制。

(1)在文件lh2中获取原子a第一次靠近壁面时的坐标，即在文件lh2中获取原子a与壁面的距离第一次小于距离阈值

(2)在文件lh2中获取原子b第一次靠近壁面时的坐标，即在文件lh2中获取原子b与壁面的距离第一次小于距离阈值

(3)记录该O

(4)依据lh3.txt文件中的记录信息，统计本模拟过程(入射动能0.05eV，壁面温度1000K)的E-R机制和L-H机制的比例，然后依次进行其他六种工况(入射动能分别为0.1eV、0.5eV、1eV、2eV、3eV和4eV，壁面温度为1000K)的反应机制的统计，直至全部工况统计完成，统计后的催化重组比例图如图7所示。

综上，通过对所有双原子气体分子在某一工况下通过E-R机制和L-H机制重组的数目进行统计，并获取不同机制在同一工况下所占的比例，进而分析在不同工况下E-R机制和L-H机制对壁面催化的作用，从而可以根据飞行器运行的不同环境条件对所使用的壁面材料进行选择。

本发明提供的上述对反应分子动力学催化重组机制发生次数的统计方法，根据通量边界模型输出的轨迹文件，使用MATLAB程序进行后处理，研究E-R机制和L-H机制并统计两种反应机制所占比例，为高超声速下使用反应分子动力学对气固界面反应机制的研究提供了可靠方法。从反应分子动力学模拟过程输出的dump文件中提取数据，通过E-R重组机制原理对不同的原子催化重组过程进行机理判定，探究不同产物随条件变化重组机制的变化规律，设计了一种判定双原子气体分子是否形成的方法，弥补了通过lammps自带的species.out命令无法对通量边界模型进行准确的种类统计的缺陷，为反应分子动力学模拟探索各种工况下气固界面反应机制、明确产物变化的规律提供了理论基础。并且，由于在模拟计算过程中分别对不同步骤的计算结果进行独立输出，因此，输出结果亦可用于特定时刻原子/分子种类和数量的统计中。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。