A method for automatic calculation of catalytic materials and reaction paths

Abstract

The present invention relates to a method for automated calculation of catalytic materials and reaction pathways, comprising the following steps: automated sectioning; surface active site identification; enumeration of intermediates and elementary reactions; construction of intermediate adsorption configurations; placement of transition state structures; preparation of DFT calculation files; construction of reaction pathways; and microscopic dynamics. Compared with the prior art, the method of the present invention uses saturation to predict the stability of different surface terminals, thereby improving sectioning efficiency. The introduction of "genetic bonds" enables automated construction of reasonable adsorption models, improving computational efficiency and saving computational costs. The catalytic reaction pathways are comprehensively considered, and the method can be linked to different computing software with a high degree of integration. Furthermore, the related method has a wider range of application scenarios, which is expected to effectively improve the computational efficiency of catalytic materials and catalytic reactions, and promote theoretical calculations from a "workshop-style" approach to an automated / standardized level.

Description

Technical Field

[0001] The present invention relates to the technical field of catalytic materials, and in particular to a method for automatic calculation of catalytic materials and reaction paths. Background Art

[0002] The development of high-performance catalytic materials is of great research significance. The rational design of catalysts is the dream of catalytic scientists. Compared with the traditional experimental trial and error method, the rational design strategy based on first-principles calculations has been proven to accelerate the development of catalytic materials. However, due to the diverse composition / structure of heterogeneous catalysts and the complex catalytic reaction network, the traditional theoretical design method has many shortcomings, resulting in the rational design of catalysts still relying on expert experience and knowledge; for example, the calculation of relevant material properties and catalytic elementary reactions is still mainly completed manually, which is inefficient. At the same time, under human interference, there are often inconsistent calculation accuracy and error standards in the theoretical calculation process, which brings unavoidable problems to the construction of subsequent standardized databases and scientific research. Therefore, it is necessary to develop an automated calculation method for the physical properties of catalytic materials and reaction pathways, realize high-throughput calculation of material / reaction data with less human participation, and greatly accelerate the construction of large-scale catalytic databases.

[0003] Although automated cross-sectioning of catalytic materials has been implemented in some programs (e.g., ASE and CatKit), its selection of surface terminations for complex materials (e.g., metal oxides and metal nitrides) remains insufficiently intelligent, requiring manual intervention to select stable terminations, resulting in significant inefficiencies. Furthermore, the identification of surface active sites is prone to bias or errors (e.g., the surface / subsurface distinction is ambiguous, making it difficult to accurately locate reactive sites), and the matching compatibility between adsorbed molecules and surface sites is also imperfect. These factors lead to significant errors in the placement of intermediate adsorption structures and transition state structures, requiring significant additional resources for theoretical calculations. Furthermore, no program currently offers a one-stop automated calculation method for outputting the optimal reaction pathway from a primitive unit cell cross-section. This results in a complex and complex catalytic data collection process from various sources, hindering data cleaning and accumulation and hindering the development of data-driven intelligent materials design.

[0004] CN111177915B provides a high-throughput calculation method and system for catalytic materials based on big data and machine learning algorithms. The method comprises the following steps: screening a catalytic material to be confirmed that meets the target catalytic performance based on a catalyst structure-activity relationship model constructed based on big data and machine learning algorithms; and determining that the catalytic material to be confirmed is a catalytic material that achieves the target catalytic performance if the deviation between the predicted result and the calculated result of the catalytic performance of the catalytic material to be confirmed is within a predetermined deviation range. However, the implementation of this method requires the accumulation of a large amount of standardized quantitative calculation and simulation data in the early stage, and it is unable to achieve complex catalytic reaction mechanism evaluation and prediction.

[0005] Therefore, in the context of the rapid development of high-performance computers and big data, there is an urgent need for a method for automated calculation of catalytic materials and reaction pathways to improve the computational efficiency of catalytic materials and catalytic reactions and to promote theoretical calculations from a "workshop-style" to an automated / standardized level. Summary of the Invention

[0006] The purpose of the present invention is to provide a method for automatic calculation of catalytic materials and reaction paths in order to overcome the defects of the above-mentioned prior art. Against the background of the rapid development of high-performance computers and big data, the invented method for automatic calculation of catalytic materials and reaction paths is expected to effectively improve the calculation efficiency of catalytic materials and catalytic reactions, and promote theoretical calculations from "workshop-style" to automation / standardization level.

[0007] The purpose of the present invention can be achieved by the following technical solutions:

[0008] The object of the present invention is to provide a method for automatically calculating catalytic materials and reaction paths, the method comprising the following steps:

[0009] Step 1: Based on the unit cell structure and physical properties of the material database, pre-cutting is performed to obtain the surface structure, and then weighted statistics are performed on the total saturation of the material and the surface atomic unsaturation;

[0010] Step 2: Based on the surface structure obtained in step 1, identify all possible unsaturated sites on the surface as active sites and calculate the degree of unsaturation of each active site. Select the active site with the largest degree of unsaturation as the candidate adsorption site, and exclude sites with the same chemical environment;

[0011] Step 3: Enumerate intermediates and elementary reactions, generate all possible intermediates based on the maximum values ​​of elements in reactants and products, and construct elementary reactions using the intermediate relationship matrix;

[0012] Step 4: Construct the intermediate adsorption configuration;

[0013] Step 5: Determine whether it is necessary to examine the transition state reaction. If so, place the transition state structure.

[0014] Step 6: Based on the intermediate adsorption configuration constructed in step 4, match the specific calculation parameters of different material elements and complete the modification of the input files of the relevant calculation tasks;

[0015] Step 7: Process the elementary reaction list provided in the step into point form, i.e., nodes. Based on the provided main reactants and products, starting from the reactant nodes, obtain all possible reaction paths involved in the reaction, and calculate the energy barriers and enthalpy changes of all reaction paths;

[0016] Step 8: The calculated reaction path energy is transferred to CATKINAS, and the key kinetic factors are calculated using the catalytic material structure performance and reaction path automatic calculation program and the microscopic dynamics general program CATKINAS software.

[0017] Furthermore, the method for automatically calculating the catalytic material and the reaction path specifically comprises the following steps:

[0018] Step 1: Automatic sectioning. Based on the unit cell structure and physical properties of the material database, different materials are pre-sectioned, and then the total saturation (S all ) and surface atomic unsaturation (S slab ) to perform weight statistics.

[0019] Step 2: Identification of surface active sites. Based on the surface structure obtained in step 1, all possible unsaturated sites on the surface are identified as active sites (top, bridge, hollow) with the maximum bond length between atoms in the original unit cell of the material as the standard, and the degree of unsaturation (S) of each active site is calculated. active-site ), select S active-site The largest site was selected as a candidate adsorption site, and sites with the same chemical environment were excluded to prevent duplicate counting.

[0020] Step 3: Enumerate intermediates and elementary reactions. Read the overall equation for the catalytic reaction, generate all possible intermediates based on the maximum values ​​of the elements in the reactants and products, and then clean these intermediates based on basic chemistry knowledge (element valence information) (intermediate cleaning). Then, use the intermediate relationship matrix to construct elementary reactions, and again clean out unreasonable elementary reactions (elementary reaction cleaning).

[0021] Step 4: Construction of intermediate adsorption configuration. By introducing a "gene bond", the surface site adsorption direction (host) is calibrated and the adsorption molecule configuration (guest) is initialized, and the host and guest are matched to achieve a more precise placement of the adsorption structure. For the surface, the direction of the unsaturated bonds ("gene bonds") of all atoms on the surface is obtained by expanding the model by 1 times in the z-axis direction; for the adsorbed intermediate, all unsaturated bonds of the intermediate molecule are confirmed as "gene bonds" based on the pre-constructed small molecule configuration database. Finally, combined with the mechanism of step-by-step matching, different adsorption configurations are fine-tuned according to the number and direction of the "gene bonds" of the surface adsorption sites and adsorbed molecules.

[0022] Step 5: Place the transition state structure (for reactions that do not require transition state examination, skip this step and jump to step 6. Use the NEB "interpolation point" method to interpolate 11 points in the "initial → final state" reaction path, select the 6th structure as the pre-guessed transition state configuration, and use the constrained search method to optimize the search to obtain the final transition state structure.

[0023] Step 6: Prepare the DFT calculation file. Based on the structural model constructed in steps 3 and 4 (based on the intermediate adsorption configuration constructed in step 4), automatically match the specific calculation parameters of different material elements, and complete the input file modification of related calculation tasks with one click. Based on the real-time supervision of the server cluster computing nodes, realize the automatic delivery and real-time collection of calculation jobs, and organize the relevant energy data and structural file data into the constructed catalytic reaction database. Complete the energy barrier and enthalpy change calculations for all elementary reactions (initial state, final state, transition state).

[0024] Step 7: Reaction Path Construction. The list of elementary reactions provided in Step 3 is converted into point form ("nodes"). Based on the provided main reactants and products, starting from the reactant node, the recursive search for the next node is continued using the newly accessible node as the new starting point until the product node is reached. "Backtracking" is not allowed during the recursive process. After path cleaning, all possible reaction paths involved in the reaction (composed of different elementary reactions) are finally obtained. Finally, the energy of all reaction paths (including energy barriers and enthalpy changes) is calculated.

[0025] Step 8: Microdynamics. Link the automated calculation program for catalytic material structure, properties, and reaction pathways with CATKINAS, a self-developed general-purpose microdynamics program. Automatically transfer the calculated reaction pathway energy (including reaction energy barriers) to CATKINAS, enabling interaction from DFT calculations to kinetic simulations. This quantitatively evaluates the rates of different reaction pathways under actual operating conditions, deriving the optimal reaction pathway. Key kinetic factors, such as surface intercalant coverage and rate-determining steps, are then verified and output to guide catalyst material modification design.

[0026] In step 1, the pre-cut surface will be re-expanded according to the size of the cut surface lattice (a, b) so that there are sufficient sites on the surface to adsorb small molecules. The weighted statistical calculation formulas for the total saturation of the material and the unsaturation of the surface atoms are S all =∑ allatoms CN and S slab =∑ slabatoms (CN max -CN) / CN max , where CN and CN max They represent the coordination number of the corresponding atom and its maximum coordination number in the bulk phase, respectively.

[0027] Furthermore, in step 2, identifying all possible unsaturated sites (top, bridge, hollow) on the surface is specifically performed by searching for all atoms that may jointly constitute bridge / hollow sites with the longest bond length as the radius for each atom on the surface. The calculation formula for active site unsaturation is: Where m is the number of atoms corresponding to each type of active site (i.e., top = 1, bridge = 2, hollow = 3).

[0028] Furthermore, in step 3, intermediate cleaning includes but is not limited to element valence mismatch, peroxides and other unreasonable intermediates; elementary reaction cleaning includes but is not limited to replacement reactions, reactions in which the bonding relationship between reactants and products is a non-subset relationship (i.e., the "elementary reaction" involves the breaking of a bond to generate a non-one-step direct conversion).

[0029] Furthermore, the small molecule configuration database in step 4 mainly contains basic hydrogen-saturated small molecules such as H2O, NH3, CH4, etc., and records their basic structural information. The rotation initialization is mainly performed with the help of bond angles and arrangements.

[0030] Furthermore, the step-by-step matching mechanism in step 4 is specifically manifested as giving priority to sites on the surface with the same degree of unsaturation as the bonding atoms in the adsorbed molecules (such as hollow sites corresponding to unsaturation 3). If the surface cannot provide this type of site, the low-saturation sites are matched downward in sequence (i.e., hollow→bridge→top).

[0031] Furthermore, the fine-tuning of different adsorption sites is specifically manifested as follows: (i) top site, rotating the adsorbed molecule so that its “gene bond” coincides with the “gene bond” direction of the adsorption site, and the distance is the sum of the atomic radii of the two bonding parties; (ii) bridge site, projecting the adsorbed molecule and the surface bridge site onto the plane to obtain the corresponding plane vector, calculating the angle between the two vectors and rotating the adsorbed molecule so that the adsorption molecule plane is perpendicular to the bridge site vector, and then matching the two “gene bonds” to coincide, with an initial distance of (iii) Hollow site: rotate the adsorbed molecule so that its “gene bond” coincides with the normal vector of the adsorption site plane. The initial distance is

[0032] Furthermore, the initial state and the final state in step 5 are the reactants and the products in the elementary reaction as the initial state and the final state of the reaction respectively. For those containing multiple species, a co-adsorption method is adopted to build the structure.

[0033] Furthermore, the specific calculation parameters in step 6 include, but are not limited to, U value, spin, and magnetism. Related calculation tasks include, but are not limited to, the required input files for software such as VASP, Gaussian, and CATKINAS. Real-time monitoring is based on the Slurm job scheduling system, with a monitoring interval of 60-300 seconds and a maximum of 20 concurrent jobs.

[0034] Furthermore, the path cleaning implementation method in step 7 is to add up all elementary reactions in each reaction path, compare them with the overall equation of the catalytic reaction, and eliminate inconsistent reaction paths.

[0035] Furthermore, in step 8, the reaction path automated calculation program refers to the calculation program of steps 1 to 7. Compared with the prior art, the present invention has the following beneficial effects:

[0036] 1) The method for automatic calculation of catalytic materials and reaction paths provided by this technical solution combines the total saturation of the material and the unsaturation of surface atoms to predict the stability of different surface terminals, which greatly improves the cutting efficiency.

[0037] 2) The method for automatic calculation of catalytic materials and reaction pathways provided by this technical solution and the introduction of "gene keys" can automatically build a more reasonable adsorption model, which can improve calculation efficiency and save calculation costs.

[0038] 3) The method for automatic calculation of catalytic materials and reaction pathways provided by this technical solution takes the catalytic reaction pathways into comprehensive consideration.

[0039] 4) The method for automated calculation of catalytic materials and reaction pathways provided by this technical solution can be linked with different computing software and has a high degree of integration, greatly reducing the workload of scientific research staff. At the same time, the relevant methods have more application scenarios (such as precise sectioning, reaction network construction, etc.). BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 Flow chart of the main steps of the method for automatic calculation of catalytic materials and reaction pathways of the present invention.

[0041] Figure 2 Schematic diagram of the crystal structure of traversing different terminal surfaces in Example 2 of the present invention.

[0042] Figure 3 Schematic diagram of surface site identification in Example 2 of the present invention, wherein (a) the RuO2(110) surface active site identified by the Catkit (ASE) program, and (b) the RuO2(110) surface active site identified by the automatic site identification method of the present invention.

[0043] Figure 4 Schematic diagram of the structure of the complex surface (Ni2Mo3N) adsorption intermediate in Example 3 of the present invention, including (a) the problem of the complex surface adsorption intermediate structure of Catkit / ASE, and (b) the adsorption structure arranged based on the "gene bond" method of the present invention.

[0044] Figure 5 This is an example diagram of the process of constructing the transition state structure of the elementary reaction (2) in Example 4 of the present invention. DETAILED DESCRIPTION

[0045] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments. This embodiment is implemented based on the technical solution of the present invention, and provides a detailed implementation method and specific operation process, but the protection scope of the present invention is not limited to the following embodiments.

[0046] Any features not explicitly stated in the present invention, such as component models, material names, connection structures, control methods, etc., shall be deemed as common technical features disclosed in the prior art.

[0047] Example 1

[0048] like Figure 1 As shown, this embodiment provides a method for automatically calculating catalytic materials and reaction paths, including the following steps:

[0049] Step 1: Automated sectioning. Based on the unit cell structure of the material database and the crystal plane index corresponding to the strongest peak in the corresponding XRD, the Python-based ASE software (Atomic Simulation Environment) was used to pre-section the different materials, and then the total saturation (S all ) and surface atomic unsaturation (S slab ) to perform weighted statistics. The section will be re-expanded p(x, y) according to the size of the section lattice (a, b) so that there are sufficient sites on the surface to adsorb small molecules, for example: and is the lattice size of the actual model. The weighted statistical calculation formulas for the total saturation of the material and the unsaturation of the surface atoms are S all =∑ allatoms CN and S slab =∑ slabatoms (CN max -CN) / CN max , where CN and CN max They represent the coordination number of the corresponding atom and its maximum coordination number in the bulk phase, respectively.

[0050] Step 2: Identification of surface active sites. Based on the surface structure obtained in step 1, all possible unsaturated sites on the surface are identified as active sites (top, bridge, hollow) with the maximum bond length between atoms in the original unit cell of the material as the standard, and the degree of unsaturation (S) of each active site is calculated. active-site ), select S active-site The largest site is selected as a candidate adsorption site, and sites with the same chemical environment are excluded to prevent repeated calculations. Identify all possible unsaturated sites (top, bridge, hollow) on the surface. Specifically, for each atom on the surface, take it as the center and the longest bond length as the radius to search for all atoms that may jointly constitute bridge / hollow sites. The calculation formula for active site unsaturation is: Where m is the number of atoms corresponding to each type of active site (i.e., top = 1, bridge = 2, hollow = 3).

[0051] Step 3: Enumeration of intermediates and elementary reactions. Read the overall equation of the catalytic reaction, generate all possible intermediates based on the maximum values ​​of the elements in the reactants and products, and clean the intermediates based on the element valence information; then, use the intermediate relationship matrix to construct the elementary reactions, and clean up the unreasonable elementary reactions again. Intermediate cleaning includes but is not limited to removing unreasonable intermediates such as element valence mismatch; elementary reaction cleaning includes but is not limited to removing replacement reactions and reactions where the bonding relationship between reactants and products is a non-subset relationship (i.e., the "elementary reaction" involves the breaking of a bond to generate a non-one-step direct conversion).

[0052] Step 4: Construction of the intermediate adsorption configuration. By introducing "genetic bonds," the surface adsorption direction (host) is calibrated and the configuration of the adsorbed molecule (guest) is initialized. Host-guest matching is used to achieve more precise placement of the adsorbed structure. For the surface, the z-axis of the model is expanded by a factor of 1 to determine the directions of the unsaturated bonds ("genetic bonds") of all surface atoms. For the adsorbed intermediate, all unsaturated bonds of the intermediate molecule are identified as "genetic bonds" based on a pre-built small molecule configuration database. Finally, a hierarchical matching mechanism is employed to fine-tune different adsorption configurations based on the number and orientation of "genetic bonds" between the surface adsorption site and the adsorbed molecule. The small molecule configuration database primarily contains basic hydrogen-saturated small molecules such as H2O, NH3, and CH4, recording their basic structural information. Rotational initialization is primarily performed using bond angles and configurations. The hierarchical matching mechanism prioritizes sites on the surface where the unsaturation of the bonding atoms in the adsorbed molecule is the same (e.g., hollow sites corresponding to an unsaturation of 3). If the surface does not provide such sites, the matching proceeds downwards to less saturated sites (i.e., hollow → bridge → top). The fine-tuning of different adsorption sites is specifically manifested as follows: (i) top site, rotating the adsorbed molecule so that its "gene bond" coincides with the "gene bond" direction of the adsorption site, and the distance is the sum of the atomic radii of the two bonding parties; (ii) bridge site, projecting the adsorbed molecule and the surface bridge site onto the plane to obtain the corresponding plane vector, calculating the angle between the two vectors and rotating the adsorbed molecule so that the adsorbed molecule plane is perpendicular to the bridge site vector, and then matching the two "gene bonds" to coincide, with an initial distance of (iii) Hollow site: rotate the adsorbed molecule so that its “gene bond” coincides with the normal vector of the adsorption site plane. The initial distance is

[0053] Step 5: Place the transition state structure (for reactions that do not require transition state examination, skip this step and jump to step 6. Use the Nudged Elastic Band (NEB) "interpolation" method in the ASE library to interpolate 11 points in the "initial→final state" reaction path, select the 6th structure as the pre-guessed transition state configuration, and use the constrained search method to optimize the search to obtain the final transition state structure. The initial and final states are the reactants and products in the elementary reaction as the initial and final states of the reaction, respectively. For reactions containing multiple species (such as OH+H), the structure is constructed by co-adsorbing the intermediates on the surface.

[0054] Step 6: DFT calculation file preparation. Based on the structural model built in steps 3 and 4, the specific calculation parameters of different material elements are automatically matched, and the input file modification of related calculation tasks is completed. Based on the real-time supervision of the server cluster computing nodes, the automatic delivery and real-time collection of calculation jobs are realized, and the optimized energy of the initial state, transition state, and final state structures and the corresponding structural coordinate files are stored in the constructed catalytic reaction database. Complete the reaction energy calculation, that is, complete the energy barrier and enthalpy change calculation of all elementary reactions (initial state, final state, transition state). Specific calculation parameters include but are not limited to U value, spin, magnetism, etc., and related calculation tasks include but are not limited to the necessary input files of software such as VASP, Gaussian, CATKINAS, etc. Real-time supervision is based on the Slurm job scheduling system, the supervision time interval is 60-300s, and the upper limit of the number of jobs allowed to exist at the same time is 20.

[0055] Step 7: Reaction path construction. The list of elementary reactions provided in step 3 is processed into point form ("nodes"). Based on the main reactants and products provided, starting from the reactant node, the new node that can be reached is used as the new starting point to continue recursively searching for the next node until the product node is reached. "Backtracking" is not allowed during the recursive process. After path cleaning, all possible reaction paths involved in the reaction (composed of different elementary reactions) are finally obtained, and finally the energy barriers and enthalpy changes of all reaction paths are calculated. The path cleaning implementation method is to add up all the elementary reactions in each reaction path, compare them with the total equation of the catalytic reaction, and eliminate the reaction paths that do not match.

[0056] Step 8: Microdynamics. Link with CATKINAS, a general microdynamics program, and automatically transfer the calculated reaction path energies (including reaction barriers) to CATKINAS, enabling interaction from DFT calculations to kinetic simulations. This quantitatively evaluates the rates of different reaction paths under actual operating conditions, verifies and outputs key kinetic factors such as surface intercalant coverage and rate-determining steps of the catalytic material, and guides the design of catalyst material modifications.

[0057] The above method has the following advantages:

[0058] 1) Using saturation to predict the stability of different surface terminals significantly improves the efficiency of the cut. Simultaneously, calculating the saturation of active sites and adsorbed molecules improves the rationality of model construction and saves computational costs for subsequent structural optimization.

[0059] 2) The introduction of the "gene bond" greatly improves the matching compatibility between the adsorbate and the catalytic site. Combined with the step-by-step adsorption matching mechanism and the bond angle analysis of the small molecule database, it realizes the fine-tuning of the configuration of the adsorbed molecules on different adsorption sites, greatly improving the success rate of constructing reasonable adsorption structures and avoiding the generation of bad structures and the waste of computing resources.

[0060] 3) Cleaning of intermediates, elementary reactions and reaction paths improves computational efficiency and saves computational costs.

[0061] 4) The search for transition state structures combines the path sampling method with the independently developed constrained search method, which greatly reduces the computing resources required by the traditional NEB method and improves the efficiency of transition state search.

[0062] 5) Comprehensive consideration of catalytic reaction pathways, especially the investigation of complex reaction networks, has unique advantages, avoiding omissions that may be caused by artificially assumed pathways.

[0063] Example 2

[0064] This embodiment uses step 1 and step 2 as an example for calculation, and the specific process is as follows:

[0065] Automatic sectioning (110) of RuO2 unit cell and identification of surface active sites.

[0066] For a given RuO2 unit cell file (accept cif or POSCAR format), set the slice Miller index (110), the atomic layer defaults to 8 (layer=8), the fixed layer defaults to half an atomic layer (fixed=4), and the vacuum layer defaults to (All parameters above are provided as input parameters in the interface.) Pre-cutting was performed. The initial cut surface model was obtained, with unit cell parameters of (a: 3.14, b: 6.43, c: 18.92), totaling 18 atoms. Considering the small initial surface area, automatic cell expansion p(3, 1) was performed. After cell expansion, the model had unit cell parameters of (a: 9.42, b: 6.43, c: 18.92), totaling 54 atoms.

[0067] Then the different terminal surfaces are traversed, such as Figure 2 As shown in FIG, there are three types of terminal surfaces (iterm=1, 2, 3). The total saturation (S all ) and surface atomic unsaturation (S slab ) to perform weighted statistics (for Ru, CN = 5 / 6, CN max =6; for O, CN=1 / 2 / 3, CN max =3):

[0068] iterm=1:S all =∑ allatoms CN=234,S slab =∑ slab atoms (CN max -CN) / CN max =4.125

[0069] iterm=2:S all =∑ allatoms CN=234,S slab =∑ slab atoms (CN max -CN) / CN max =4.125

[0070] iterm=3:S all =∑ allatoms CN=240,S slab =∑ slab atoms (CN max -CN) / CN max =2.75

[0071] Specifically, S all The larger (S slab The smaller the value, the fewer unsaturated bonds in the material (upper and lower surfaces), which theoretically makes the surface more stable. This surface is therefore more easily exposed and therefore preferred as the substrate for subsequent catalytic reactions. This method ensures the stability of the selected surface from two perspectives. It rapidly selects the endpoint through a simple saturation calculation, greatly improving the efficiency of the sectioning process. It also provides an interface for surface energy calculations for verification purposes.

[0072] Then the surface sites are identified. Although the Catkit (ASE) program can identify most of the active sites on the RuO2 (110) surface, it lacks intelligent judgment, which makes the display and subsequent calculation efficiency low. In addition, some sites are misidentified, such as surface saturated tri-coordinate oxygen, Ru-O mixed sites, etc., which bring great problems to the construction of subsequent intermediate adsorption structure and transition state structure, that is, the structure is unreasonable (such as Figure 3 a).

[0073] The automatic site identification method of the present invention first uses the longest bond length between atoms in the original unit cell of the material as the standard ('Ru_O': 'Ru_Ru': ), and additionally provide As a tolerance factor, all atoms on the surface whose distance is less than the bond length standard can be considered as active sites. For each atom on the surface, with it as the center and the longest bond length as the radius, all atoms that may jointly constitute bridge / hollow sites are searched. If the result is 0, it becomes a top site independently. In addition, the present invention removes the same sites by comparing the chemical environment of each site to prevent repeated calculations and waste of resources. Finally, the automatically obtained RuO2(110) surface has a total of one bridge site and two top sites (such as Figure 3b). Then, the unsaturation S of the active site site The most unsaturated sites (i.e., more empty orbitals, easier to adsorb / activate small molecules) are selected as candidate adsorption sites. Take the unique bridge site as an example, and mark its site name as "Ru_Ru", and the site unsaturation Mark the site coordinates Ru1: [3.21, 3.14, 17.67], Ru2: [3.21, 0.00, 17.67], take the average of the two plane coordinates ((Ru1+Ru2) / 2) and provide the coordinates in the z-axis direction. The initial distance is used as the coordinate of the basic adsorption point of the small molecule: [3.21, 1.57, 17.67+1.5], which provides a basis for the subsequent placement of small molecules.

[0074] By comparison, it can be found that the method of the present invention can locate potential active sites more accurately and has higher efficiency than other methods without missing any.

[0075] Example 3

[0076] This embodiment uses step 4 as an example for calculation, and the specific process is as follows:

[0077] Complex surface (Ni2Mo3N(100) as an example) adsorption intermediate (CH x )’s structure is shown in comparison.

[0078] For other existing automated programs (Catkit / ASE), in addition to the aforementioned possible misidentification of adsorption sites, there are also problems such as incorrect small molecule adsorption height causing molecules to enter the subsurface, and the small molecule adsorption direction / angle is not intelligent enough (such as Figure 4 a), none of them can reasonably and automatically place the adsorbed molecules to the correct position, resulting in low efficiency of subsequent structure optimization and large errors in energy calculation.

[0079] For the method of the present invention, the surface active site identification is first performed on the Ni2Mo3N unit cell structure file in the same manner as in Example 1. Before matching the small molecule adsorption, the surface site adsorption direction (host) is calibrated by introducing a "gene key" and the adsorbed molecule configuration (guest) is initialized. The host-guest matching is used to achieve a more accurate placement of the adsorption structure (such as Figure 4 b).

[0080] First, the entire model is replicated in the z-axis direction (original layer, replica layer). In the replica layer, atoms with the same xy coordinates as the surface atoms of the original layer and a constant z-axis distance are found. These atoms are called surface atoms in the replica layer. By comparing the bonding information between the surface atoms of the original layer and the replica layer, the difference in coordinates between the two is all the unbonded directions of the surface atoms of the original layer (i.e., "genetic bonds"). This "genetic bond" is more consistent with the empty orbital direction of the atom's current chemical environment and will have a higher degree of hybridization, i.e., binding energy, when bonding with small molecules. Using the above method, the "genetic bond" on the surface of Ni2Mo3N(100) is:

[0081] 40-Ni: [0.11428785, -1.27503548, 2.4357831]

[0082] 41-Mo: [-1.16074763, 1.16074763, 2.32149525]

[0083] 41-Mo: [1.27503548, -0.11428785, 2.4357831]

[0084] 42-Mo: [2.16931792, 0.50428515, 1.66503277]

[0085] …

[0086] 43-Ni: [-0.77999459, 1.66503277, 1.66503277]

[0087] 43-Ni: [-0.3899973, -0.89428244, 2.55931521]

[0088] 44-N: [0.50428515, 1.66503277, 1.16074763]

[0089] …

[0090] Then, we start to initialize the small molecule configuration. Through the pre-built small molecule database, we know the bond angle information of most small molecules (such as CH4: 109.28° in this example) and their spatial configuration rules (such as tetrahedral configuration in this example). From this, we can obtain the unbonded directions of all C-containing intermediates, that is, the "genetic bond" of the small molecule. Figure 4As shown in b, the original molecule (i.e., the reaction intermediate) is first read, and the most unsaturated atom is first confirmed as the bonding atom through the unsaturation calculation, and one of the "gene bonds" is selected to form a bond. The molecule is initialized to its "gene bond" direction of (0, 0, -1) by spatial rotation. Subsequently, according to the unsaturation of the adsorbed small molecule (CH: 3, CH2: 2, CH3: 1), the best adsorption sites on the surface (CH: hollow, CH2: bridge, CH3: top) are preferentially matched. For the top site and the hollow site, it is only necessary to overlap the surface and the molecule "gene bond". For the bridge site, additional special processing is required: the surface bridge site diatomic and small molecule coordinates are projected onto the xy plane, and the "site vector" and "molecule original vector" can be obtained respectively, and the "molecule target vector" is the vertical vector of the "site vector", from which it can be known that the angle θ that the adsorbed molecule needs to rotate is made, so that the adsorption configuration of CH2 is more reasonable. The automated construction method of the adsorption structure of the present invention has "expert experience knowledge" to a certain extent, and the placement structure is more reasonable, saving costs and improving efficiency for subsequent calculations.

[0091] Example 4

[0092] This embodiment provides a method for automatically calculating catalytic materials and reaction pathways, which realizes the full automatic calculation of the H2O electrolysis (HER) reaction pathway of Pt(111), specifically comprising the following steps:

[0093] For a given Pt unit cell file, a (111) direction pre-cut is performed. Through the previous steps, the model is first automatically expanded to p(4, 3) so that the model lattice parameters are (a: 11.25, b: 7.31, c: 16.89). Then, through the unsaturation investigation, the surface terminal that is most easily exposed is selected as the surface model for subsequent reactions: iterm1: S all =∑ all atoms CN=504,S slab =∑ slab atoms (CN max -CN) / CN max =6.

[0094] Subsequently, the surface active sites were identified, and the maximum bond length between atoms in the primitive unit cell ('Pt_Pt': ), and additionally provide As a tolerance factor, the distance between atoms on the surface that is less than the bond length standard can be considered as an active site. By comparing the chemical environment, the top / bridge / hollow sites on the Pt(111) surface are identified, and the corresponding site unsaturation is The corresponding adsorption point coordinates are:

[0095] Top: [0.00, 0.00, 11.50]

[0096] Bridge: [1.43, 0.00, 11.00]

[0097] Hollow: [1.43, 0.83, 11.00]

[0098] At the same time, the surface "gene key" is:

[0099] 1-Pt: [0.00, 1.62350955, 2.29598923]

[0100] 1-Pt: [-1.40600052, -0.81175478, 2.29598923]

[0101] 1-Pt: [1.40600052, -0.81175478, 2.29598923]

[0102] …

[0103] 11-Pt: [-1.40600052, -0.81175478, 2.29598923]

[0104] 11-Pt: [1.40600052, -0.81175478, 2.29598923]

[0105] 11-Pt: [0.00, 1.62350955, 2.29598923]

[0106] The overall reaction equation for H2O dissociation is: H2O→1 / 2H2+OH - , the maximum elemental composition of all intermediates is determined to be (H:2, O:1), and a total of five intermediates are generated (H, O, H2, OH, H2O). These intermediates all meet basic chemical knowledge and do not require cleaning. By additionally considering the instructions of adsorption and desorption reactions, a total of the following six elementary reactions are generated:

[0107] (1) *H+*H→*H2

[0108] (2) *H+*OH→*H2O

[0109] (3)*H2+*OH→*H+*H2O

[0110] (4)*+H2O(aq)→*H2O

[0111] (5)*H2→*+H2(aq)

[0112] (6)*OH→*+OH(aq)

[0113] Where (aq) represents that the small molecule does not interact with the surface, that is, it is in the liquid phase (that is, steps 4-6 are adsorption-desorption reactions). By reading each elementary reaction, the reactant is the initial state (IS) and the product is the final state (FS), and the transition state structure (TS) of the corresponding elementary reaction is generated by path sampling and interpolation. Figure 5 The figure shows an example of the transition state structure construction process for the elementary reaction (2). The TS-6 structure was ultimately selected as the initial transition state structure, and an optimized search was performed using a self-developed constrained search method. During the automated adsorption structure construction phase, all intermediates were initialized using the small molecule database (e.g., H2O: 104.5°, with a planar "V" spatial configuration in this example), and then constructed by matching the "genetic bonds" between the surface and the small molecule.

[0114] The slurm job scheduling system monitors the number of remaining cores in the server cluster in real time, and the VASP calculation input files (including INCAR, KPOINTS, POSCAR, POTCAR) automatically generated according to the above structure are delivered to the server for calculation through the job submission system. At the same time, the completion status of the delivered jobs is monitored every 60 seconds and status information is returned in a timely manner. For the optimized jobs, the optimized structure is copied and the energy information is saved. After all the elementary reaction jobs are optimized, the reaction energy barrier (E a ) and enthalpy change (ΔH), as shown in Table 1:

[0115] Table 1 Reaction energy barriers and enthalpy changes of each elementary step.

[0116]

[0117] Note: (m) indicates that the intermediate is not bonded to the catalyst surface and exists in molecular form.

[0118] Based on the above elementary reactions, all possible complete reaction paths (essentially the addition process of different elementary reactions) are enumerated with the starting material "H2O" as the starting point and the product "H2" as the end point, as shown in Table 2:

[0119] Table 2 All possible complete reaction pathways.

[0120]

[0121]

[0122] Note: (m) indicates that the intermediate is not bonded to the catalyst surface and exists in molecular form.

[0123] From this, we can know that the optimal reaction path of HER on Pt(111) is described in path 1, with a total reaction energy barrier of 0.581 and the rate-determining step being the dissociation of H2O on the surface. In order to further verify whether the reaction path determined by energy is correct, the E of the HER reaction path is a The system automatically fills in the input file of the independently developed microdynamics CATKINAS software with ΔH and ΔH, and automatically submits the job to analyze the microdynamics of the reaction, quantitatively evaluate the rates of different reaction paths under actual working conditions, verify and output key kinetic factors such as the surface coverage of the catalytic material and the rate-determining step, and guide the modification design of the catalyst material.

[0124] The above description of the embodiments is intended to facilitate understanding and use of the invention by those skilled in the art. It will be apparent that those skilled in the art can readily make various modifications to these embodiments and apply the general principles described herein to other embodiments without requiring inventive effort. Therefore, the present invention is not limited to the above-described embodiments. Improvements and modifications made by those skilled in the art based on the disclosure of the present invention, without departing from the scope of the present invention, should be within the scope of protection of the present invention.

Claims

1. A method for automatic calculation of catalytic materials and reaction paths, characterized in that: The method comprises the following steps: Step 1: Based on the unit cell structure and physical properties of the material database, pre-cutting is performed to obtain the surface structure, and then weighted statistics are performed on the total saturation of the material and the surface atomic unsaturation; Step 2: Based on the surface structure obtained in step 1, identify all possible unsaturated sites on the surface as active sites and calculate the degree of unsaturation of each active site. Select the active site with the largest degree of unsaturation as the candidate adsorption site, and exclude sites with the same chemical environment; Step 3: Enumerate intermediates and elementary reactions, generate all possible intermediates based on the maximum values ​​of elements in reactants and products, and construct elementary reactions using the intermediate relationship matrix; Step 4: Construct the intermediate adsorption configuration; Step 5: Determine whether it is necessary to examine the transition state reaction. If so, place the transition state structure. Step 6: Based on the intermediate adsorption configuration constructed in step 4, match the specific calculation parameters of different material elements and complete the modification of the input files of the relevant calculation tasks; Step 7: Process the elementary reaction list provided in the step into point form, i.e., nodes. Based on the provided main reactants and products, all possible reaction paths involved in the reaction are obtained starting from the reactant nodes, and the energy of all reaction paths is calculated; Step 8: The calculated reaction path energy is transferred to CATKINAS, and the key kinetic factors are calculated using the catalytic material structure performance and reaction path automatic calculation program and the microscopic dynamics general program CATKINAS software.

2. The method for automatic calculation of catalytic materials and reaction paths according to claim 1, characterized in that: In step 1, the pre-cut surface will be re-expanded according to the size of the cut surface lattice (a, b) p (x, y) allows for sufficient sites on the surface to adsorb small molecules.

3. The method for automatic calculation of catalytic materials and reaction paths according to claim 1, characterized in that: In step 2, the identification of all possible unsaturated sites on the surface is to search for all atoms that may jointly constitute bridge sites or hollow sites with each atom on the surface as the center and the longest bond length as the radius.

4. The method for automatic calculation of catalytic materials and reaction paths according to claim 1, characterized in that: In step 3, the specific process of enumerating intermediates and elementary reactions is as follows: Read the overall equation of the catalytic reaction, generate all possible intermediates based on the maximum values ​​of the elements in the reactants and products, and clean the intermediates based on basic chemical knowledge; then, use the intermediate relationship matrix to construct the elementary reactions, and clean up the unreasonable elementary reactions again.

5. The method for automatic calculation of catalytic materials and reaction paths according to claim 4, characterized in that: The method of cleaning the intermediates is to remove unreasonable intermediates with mismatched element valences; The method of cleaning unreasonable elementary reactions includes removing displacement reactions and reactions in which the bonding relationship between reactants and products is a non-subset relationship.

6. The method for automatic calculation of catalytic materials and reaction paths according to claim 1, characterized in that: In step 4, the intermediate adsorption configuration is constructed as follows: By introducing "gene bonds", the adsorption direction of surface sites is calibrated and the configuration of adsorbed molecules is initialized to achieve more precise placement of the adsorption structure; for the surface, the model is expanded by 1 times in the z-axis direction to obtain the unsaturated bonds of all surface atoms as the "gene bonds" of the adsorption sites; for the adsorbed intermediates, based on the pre-constructed small molecule configuration database, all unsaturated bonds of the intermediate molecules are confirmed as the "gene bonds" of the adsorbed molecules; finally, combined with the step-by-step matching mechanism, different adsorption configurations are fine-tuned according to the number and direction of the "gene bonds" of the surface adsorption sites and adsorbed molecules.

7. The method for automatic calculation of catalytic materials and reaction paths according to claim 1, characterized in that: In step 5, the method of placing the transition state structure is as follows: Using the NEB "point interpolation" method, 11 points were interpolated in the "initial→final state" reaction path. The sixth structure was selected as the pre-guessed transition state configuration, and the "constrained search method" was used to optimize the search to obtain the final transition state structure.

8. The method for automatic calculation of catalytic materials and reaction paths according to claim 1, characterized in that: The specific process of step 6 is as follows: Based on the intermediate adsorption configuration constructed in step 4, the specific calculation parameters of different material elements are matched, the input file modification of the relevant calculation tasks is completed, and the real-time supervision of the server cluster computing nodes is based on the realization of automatic delivery and real-time collection of calculation jobs. The relevant energy data and structure file data are sorted and stored in the constructed catalytic reaction database to complete the energy barrier and enthalpy change calculations for all elementary reactions.

9. The method for automatic calculation of catalytic materials and reaction paths according to claim 1, characterized in that: The specific process of step 7 is as follows: The elementary reaction list provided in step 3 is processed into point form, i.e., "nodes". Based on the main reactants and products provided, starting from the reactant node, the reachable new node is used as the new starting point to continue recursively searching for the next node until the product node is reached. "Backtracking" is not allowed during the recursive process. After path cleaning, all possible reaction paths involved in the reaction are finally obtained, and the energy of all reaction paths is calculated.

10. The method for automatic calculation of catalytic materials and reaction paths according to claim 1, characterized in that: The specific process of step 8 is as follows: The automated calculation program for the structural properties and reaction paths of catalytic materials is linked to the general microdynamics program CATKINAS software, and the calculated reaction path energy is transferred to CATKINAS, enabling interaction from DFT calculation to kinetic simulation. This allows for quantitative evaluation of the rates of different reaction paths under actual operating conditions, verification and output of key kinetic factors, and guidance for the modification design of catalyst materials.

Images