A method and system for predicting chemical reaction pathways

Abstract

The application relates to a chemical reaction path prediction method and system, the method comprising the following steps: defining a reaction pair comprising reactants and potential products; performing joint optimization processing on the reactants and the potential products by using a molecular dynamics method and a quantum chemistry method to generate initial reactant configurations and initial product configurations; determining preliminary reaction path information connecting the initial reactant configurations and the initial product configurations by using a quantum chemistry method based on the alignment; and locating and verifying accurate transition states of the reactants and the potential products based on the preliminary reaction path information by using a quantum chemistry method. The application provides a chemical reaction path prediction method and a chemical reaction path prediction system using the method, which can improve the initial guess quality for high-precision calculation, thereby greatly improving the convergence success rate and accuracy of transition state searching, and further reducing the overall calculation cost under the premise of ensuring the prediction accuracy.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of computational chemistry, and in particular to a chemical reaction path prediction method and system. BACKGROUND

[0002] In modern chemical research and industrial applications, such as new drug design, catalyst development, new material creation, and combustion mechanism analysis, a deep understanding of complex chemical reaction networks is crucial. Theoretical computational chemistry, especially quantum mechanics-based computational methods, provides core tools for revealing reaction mechanisms, finding reaction paths, and locating transition states at the atomic level. High-precision computational methods such as density functional theory (DFT) can provide reliable energy and structural information, but their computational cost is extremely high, often increasing exponentially or higher with the number of atoms in the system.

[0003] Therefore, when faced with complex systems containing numerous atoms or involving multiple-step reactions, relying solely on high-precision methods for global reaction path exploration is computationally infeasible. This "brute force" strategy not only consumes a huge amount of time, but also often fails to converge within a reasonable time, resulting in a severe lack of coverage and greatly limiting the practical application range and efficiency of automated reaction prediction technology. SUMMARY

[0004] The present application aims to overcome the above-mentioned defects of the prior art and provide a chemical reaction path prediction method and system, which significantly improves the quality of initial guesses for high-precision calculations through a systematic joint optimization process, thereby greatly improving the convergence success rate and accuracy of transition state searches, further reducing overall computational cost while ensuring prediction accuracy.

[0005] To achieve the above-mentioned purpose, the present application adopts the following technical solutions:

[0006] A chemical reaction path prediction method, comprising:

[0007] (A) Defining a reaction pair containing reactants and potential products;

[0008] (B) Generating initial configurations using molecular mechanics methods and performing joint optimization processing of the three-dimensional geometric configurations of the reactants and potential products using a quantum chemistry method with a first computational accuracy level to generate a pair of spatially aligned initial reactant configurations and initial product configurations;

[0009] (C) Based on the aligned initial reactant configurations and initial product configurations, determining a preliminary reaction path and a preliminary transition state on the path using the quantum chemistry method with the first computational accuracy level;

[0010] (D) using a quantum chemical method of a second computational accuracy level, optimizing and verifying the preliminary transition state to obtain an accurate transition state connecting the reactant and the potential product;

[0011] wherein the second computational accuracy level is higher than the first computational accuracy level.

[0012] Preferably, the quantum chemical method of the first computational accuracy level is a semi-empirical quantum chemical method.

[0013] Preferably, the semi-empirical quantum chemical method is a GFN2-xTB method.

[0014] Preferably, the quantum chemical method of the second computational accuracy level is a density functional theory method or a post-Hartree-Fock method.

[0015] Preferably, the joint optimization process in the step (C) specifically comprises:

[0016] generating an initial three-dimensional geometry of the potential product using a force field method based on the molecular graph structure of the potential product;

[0017] using the initial three-dimensional geometry of the potential product as a template, optimizing on a force field potential energy surface corresponding to the bond-electron matrix of the reactant to generate an initial three-dimensional geometry of the reactant;

[0018] using the quantum chemical method of the first computational accuracy level, optimizing the initial three-dimensional geometries of the reactant and the potential product (20) respectively;

[0019] aligning the centers of mass of the optimized reactant and potential product by rotation and translation operations, and minimizing the root mean square deviation (RMSD) of their atomic coordinates.

[0020] Preferably, the step (C) determines the preliminary reaction path using a growing string method (GSM).

[0021] Preferably, the step (D) optimizes the preliminary transition state using a Berny optimization algorithm.

[0022] Preferably, the step (A) defines the reaction pair by applying a preset chemical bond recombination rule to the molecular graph structure of the reactant to automatically enumerate the molecular graph structure of the potential product.

[0023] Preferably, the above method further comprises:

[0024] (E) integrating the accurate transition state and its corresponding reaction path information into a reaction network (40);

[0025] (F) When there are multiple intermediates that can be used as new reactants in the reaction network, determine the exploration priority of each intermediate based on the highest activation energy barrier in the path connecting the initial reactants to each intermediate, and select the intermediate with the highest priority as the new reactant, and iteratively perform steps (A) to (E).

[0026] Further, a chemical reaction path prediction system is also provided, which comprises at least one processor;

[0027] and a memory connected in communication with the processor, the memory having stored therein instructions executable by the processor, the instructions, when executed, causing the system to implement the above-described chemical reaction path prediction method.

[0028] Compared with the prior art, the present application has significant and multiple benefits:

[0029] 1. Significantly improve the success rate of convergence and "hit rate": The present application ensures that the reactant and product configurations used for low-precision path search are accurately aligned in space through systematic joint optimization processing; the joint optimization processing adopts a unique template-induced optimization strategy: that is, taking the three-dimensional configuration of the product as the template, but optimizing on the force field potential surface of the bonding relationship (i.e. connection matrix) of the reactant to generate the initial configuration of the reactant. This seemingly "contradictory" optimization method cleverly ensures that the reactant configuration is highly similar to the product configuration in topology from the beginning, making the preliminary transition state determined in subsequent steps closer to the real accurate transition state in geometry.

[0030] 2. Further reduce the overall computing cost: Due to the significant improvement in the success rate of convergence and hit rate, invalid and time-consuming second-precision-level (such as DFT) calculation tasks are significantly reduced. This means that when exploring the same number of potential reaction paths, the total computing resources required by the present application are much less than those of the prior art.

[0031] 3. Expand the exploration range and improve network integrity: Extremely low unit path calculation cost makes it possible to explore larger and more complex chemical reaction systems. In addition, the iterative exploration strategy based on the highest activation energy barrier proposed by the present application can intelligently guide the computing resources to the most favorable reaction area in dynamics, avoiding wasting time on ineffective paths with too high energy. The combination of these two advantages enables the present application to efficiently construct a more complete and reliable chemical reaction network than the prior art, thereby being able to handle complex systems that traditional methods cannot handle.

[0032] In order to make the purpose, technical scheme and advantages of the present application more clear and explicit, the following will further detail the present application through specific embodiments and related drawings. BRIEF DESCRIPTION OF DRAWINGS

[0033] Figure 1 is a functional module schematic diagram of an automated reaction pathway prediction system according to an embodiment of the present application.

[0034] Figure 2 is a general flowchart of an automated reaction pathway prediction method according to an embodiment of the present application.

[0035] Figure 3 is Figure 2 is a detailed flowchart of a preferred implementation of step (B) "joint optimization process" in

[0036] Figure 4 is an iterative strategy schematic diagram for exploring complex reaction networks according to an embodiment of the present application. DETAILED DESCRIPTION

[0037] In order to make the objectives, technical solutions and advantages of the present application clearer, the present application will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments herein are only used to explain the present application, but not to limit the present application.

[0038] The technical solutions of the present application can be deployed on computing devices including but not limited to personal computers, high-performance computing workstations, server clusters or cloud computing platforms, etc. These devices usually include at least one processor, memory, input / output interface, etc. hardware. EMBODIMENTS

[0039] The present embodiment provides a chemical reaction pathway prediction method. Please refer to Figure 2 The method is executed by a processor, and the specific process is as follows:

[0040] Step (A) 201: Define a reaction pair containing a reactant 10 and a potential product 20.

[0041] This step is the starting point of reaction pathway prediction. First, the system receives the molecular structure information of the initial reactant 10 defined by the user through its input interface. (The molecular structure information is represented as a molecular graph data structure in the memory, which contains the types, numbers of atoms and bonding relationships between atoms (i.e. connectivity).)

[0042] In order to automatically explore chemical possibilities, the system needs to automatically generate a series of potential products 20, thereby forming a plurality of reaction pairs to be studied. In this embodiment, the reaction pair is defined by applying a preset chemical bond recombination rule to the molecular graph structure of the reactant 10 to automatically enumerate the molecular graph structure of the potential product 20.

[0043] These rules are pre-set based on basic chemical principles, for example, a common rule is the "b2f2" rule, which selectively breaks two chemical bonds in the reactant molecule and forms two new chemical bonds at the same time. The processor will traverse all valid chemical bond combinations in the molecule graph of the reactant 10, apply this rule, and generate a series of new, chemically reasonable molecule graph structures, which are the potential products 20. Compared with random search, this graph-based enumeration method can systematically and exhaustively generate all products that meet the specific reaction mode, laying a foundation for subsequent comprehensive exploration.

[0044] Step (B) 202: Generate an initial configuration using molecular mechanics method, and jointly optimize the three-dimensional geometric configuration of the reactant 10 and the potential product 20 using a quantum chemistry method of the first calculation precision level to generate a pair of spatially aligned initial reactant configuration and initial product configuration.

[0045] This step is one of the core innovations of the present application, and its purpose is to provide a pair of high-quality, geometrically matched end-point structures for subsequent path search.

[0046] In this embodiment, the quantum chemistry method of the first calculation precision level is a semi-empirical quantum chemistry method combined with molecular dynamics method, specifically, the GFN2-xTB method can be used as a semi-empirical quantum chemistry component. (The reason for choosing GFN2-xTB method is that it has achieved an excellent balance between calculation speed and description accuracy of molecular geometric configuration. Its calculation cost is much lower than that of DFT method, but it is significantly better than traditional force field method, and it is very suitable for fast and reliable geometric optimization in this stage.)

[0047] Please refer to Figure 3 The joint optimization processing step (B) 202 specifically includes the following sub-steps:

[0048] Step 301: Based on the molecule graph structure of the potential product 20, use the force field method to generate its initial three-dimensional geometric configuration.

[0049] (For each molecule graph (two-dimensional connection information) of the potential product 20 generated in step (A) 201, the processor first calls a program module based on the force field (for example, the universal force field UFF). This module can quickly convert the two-dimensional molecule graph into an initial coordinate set with a roughly reasonable three-dimensional geometric configuration. The calculation speed of the force field method is extremely fast, and it is suitable for processing the preliminary configuration generation of a large number of products.)

[0050] Step 302: Use the initial three-dimensional geometry of the potential product 20 as a template to perform optimization on the force field potential energy surface corresponding to the bond-electron matrix of the reactant 10, to generate the initial three-dimensional geometry of the reactant 10.

[0051] (This is a key technique in the present application scheme to ensure the initial alignment. The processor uses the initial three-dimensional configuration of the product 20 as a template, but switches the bonding relationship to that of the reactant 10, and then performs force field optimization on this "contradictory" structure. Due to the presence of the template, the optimized reactant 10 configuration will naturally maintain a high degree of similarity in spatial orientation with the product 20 configuration, laying a solid foundation for subsequent accurate alignment.)

[0052] Step 303: Use the quantum chemistry method of the first calculation precision level combined with molecular dynamics method to optimize the initial three-dimensional geometry of the reactant 10 and potential product 20, respectively.

[0053] The processor calls the GFN2-xTB calculation module to perform geometric optimization on the rough three-dimensional configuration generated in sub-steps 301 and 302. (This optimization process will more accurately adjust the bond length, bond angle, and dihedral angle, so that the configurations of the reactant 10 and the product 20 respectively reach their energy minimum points on the GFN2-xTB potential energy surface, obtaining more reliable energy and more refined structure of the geometric configuration.)

[0054] Step 304: Align the centers of mass of the optimized reactant 10 and potential product 20 by rotation and translation, and minimize the root mean square deviation (RMSD) of their atomic coordinates. This alignment step ensures the spatial matching of the reactant and product configurations to supplement the optimization results based on the connection matrix.

[0055] To achieve perfect spatial matching, the processor performs a rigid body transformation algorithm. It first calculates and aligns the centers of mass of the two molecules, and then finds a rotation matrix that minimizes the sum of the squares of the distances between corresponding atom pairs in the two molecules (i.e., RMSD) by rotating one of the molecules. After this step, the system finally outputs a pair of initial reactant and product configurations that are accurately aligned in space and optimized by low-precision methods.

[0056] Step (C) 203: Based on the aligned initial reactant and initial product configurations, use the quantum chemistry method of the first calculation precision level to determine the preliminary reaction path and the preliminary transition state on the path connecting the two.

[0057] In this embodiment, the growth string method (GSM) is used to determine the preliminary reaction path.

[0058] (GSM is an efficient double-ended transition state search algorithm. The processor takes the aligned reactant and product configurations output by step (B) 202 as the two endpoints of the GSM algorithm. The GSM algorithm constructs a minimum energy path (MEP) by iteratively "growing" and optimizing a series of intermediate nodes on the multidimensional potential energy surface connecting the two endpoints.) This path represents the most probable transition pathway from reactants to products on the GFN2-xTB potential energy surface. The node with the highest energy on this path is the preliminary transition state at the first level of computational accuracy. The advantage of the GSM algorithm is that it directly utilizes high-quality endpoint information, is highly efficient, and is well suited for finding preliminary reaction pathways.

[0059] Step (D) 204: The preliminary transition state is optimized and verified using a quantum chemical method at a second level of computational accuracy to obtain an accurate transition state 30 connecting the reactants 10 and the potential product 20.

[0060] This step is a high-accuracy confirmation of the low-accuracy result. In the present embodiment, the quantum chemical method at the second level of computational accuracy is a density functional theory method (DFT) or a post-Hartree-Fock method (post-HF), with DFT being a further preferred method.

[0061] The Berny optimization algorithm is used to optimize the preliminary transition state.

[0062] (The processor takes the atomic coordinates of the preliminary transition state obtained in step (C) 203 as input and initiates a transition state optimization task based on the DFT potential energy surface. The Berny optimization algorithm is a highly efficient second-order optimization method that can efficiently find a saddle point on the potential energy surface based on first- and second-order derivative (Hessian matrix) information of the energy with respect to the coordinates.) After optimization converges, a candidate transition state structure at the DFT level of accuracy is obtained.

[0063] Subsequently, the system must verify the candidate transition state. (The verification process includes two aspects: first, a vibration frequency analysis is performed, and a true transition state must have and only have one imaginary frequency in its Hessian matrix, which corresponds to the vibration mode in the reaction coordinate direction. Second, an intrinsic reaction coordinate (IRC) calculation is performed. The IRC calculation traces the minimum energy path from the candidate transition state in both the positive and negative directions until it reaches the energy minimum point. If the two endpoints of the IRC path converge to the initial reactants 10 and potential product 20, respectively, it is proved that the transition state is the accurate transition state 30 connecting the two.)

[0064] Steps (E) 205 and (F) 206: Reaction network construction and iterative exploration.

[0065] After finding one or more precise transition states 30, the method further comprises:

[0066] (E) 205: integrating the precise transition states 30 and their corresponding reaction pathway information into a reaction network 40.

[0067] (In memory, the system maintains a graph data structure to represent the reaction network 40. Molecules (reactants, products, intermediates) are nodes of the graph, while each validated reaction pathway (defined by reactants, precise transition states, products, and activation energy barriers) is a directed edge connecting the nodes.)

[0068] (F) 206: when there are multiple intermediates in the reaction network 40 that can serve as new reactants, determine the exploration priority of each intermediate based on the highest activation energy barrier in the path connecting the initial reactant to each intermediate, and select the intermediate with the highest priority as the new reactant 10, and iterate steps (A) to (E).

[0069] To efficiently explore the complex reaction space, please refer to Figure 4 The present application also provides an intelligent network exploration strategy. As shown in Figure 4 The system calculates a "exploration cost" for each intermediate node (such as intermediate 1, intermediate 2) in the network, which is defined as the highest activation energy barrier that needs to be crossed from the initial reactant (initial reactant A) to reach the intermediate. The system maintains a priority queue, always selects the intermediate with the lowest "exploration cost" as the new reactant 10, and returns to step (A) 201 to start a new round of pathway exploration.

[0070] The advantage of this strategy is that it can ensure that the computational resources are preferentially used to explore the most kinetically favorable, i.e. the lowest energy barrier, reaction region, thereby revealing the most core and most important reaction channels in the entire reaction network with the highest efficiency. Embodiment

[0071] The present embodiment provides a chemical reaction pathway prediction system. Please refer to Figure 1 The system can be a software system deployed on a server, whose hardware foundation includes at least one processor 101 and a memory 102 in communication with the processor 101. The memory 102 stores computer program instructions executable by the processor 101.

[0072] When the instructions are executed, the system implements the method as described in Embodiment One. Specifically, the system can include the following functional modules (these modules are functional divisions of program instructions):

[0073] Reaction pair defining module 11: responsible for performing step (A), receiving initial reactant information, and generating potential products according to chemical bond recombination rules.

[0074] Joint optimization and alignment module 12: responsible for performing step (B), performing systematic joint optimization processing on reactants and products, and generating aligned initial configurations.

[0075] Preliminary path search module 13: responsible for performing step (C), using algorithms such as GSM to determine preliminary reaction paths and preliminary transition states on low-precision potential energy surfaces.

[0076] Precise verification module 14: responsible for performing step (D), using methods such as DFT and Berny optimization, IRC verification, etc. to obtain precise transition states.

[0077] Network management and exploration module 15: responsible for performing steps (E) and (F), constructing and maintaining reaction networks, and intelligently performing iterative exploration according to priority strategies.

[0078] Those of ordinary skill in the art can understand that all or part of the steps in the above embodiments can be instructed by a computer program to complete the relevant hardware, and the program can be stored in a computer readable storage medium such as a hard disk, an optical disk, or a flash memory.

[0079] The above is only a preferred embodiment of the present application, and is not intended to limit the protection scope of the present application. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present application shall be included in the protection scope of the present application.

Claims

1. A method for predicting chemical reaction pathways, characterized in that, include: (A) Define a reaction pair comprising reactant (10) and potential product (20); (B) The initial configuration is generated using molecular mechanics methods, and the three-dimensional geometric configurations of the reactant (10) and the potential product (20) are jointly optimized using quantum chemical methods with the first level of computational precision to generate a pair of spatially aligned initial reactant configurations and initial product configurations. (C) Based on the aligned initial reactant and initial product configurations, a quantum chemical method of the first computational accuracy level is used to determine the initial reaction pathway connecting the two and the initial transition state along the pathway. (D) The preliminary transition state is optimized and verified using a quantum chemical method with a second level of computational precision to obtain a precise transition state (30) connecting the reactant (10) and the potential product (20); Wherein, the second calculation accuracy level is higher than the first calculation accuracy level; The joint optimization process in step (B) specifically includes: (B1) Based on the molecular diagram structure of the potential product (20), its initial three-dimensional geometry is generated using a force field method; (B2) Using the initial three-dimensional geometry of the potential product (20) as a template, optimization is performed on the force field potential energy surface corresponding to the bond-electron matrix of the reactant (10) to generate the initial three-dimensional geometry of the reactant (10). (B3) The initial three-dimensional geometry of the reactant (10) and the potential product (20) is optimized using a quantum chemical method with the first level of computational precision. (B4) Align the centroids of the optimized reactant (10) and potential product (20) by rotation and translation operations, and minimize the root mean square deviation of their atomic coordinates.

2. The chemical reaction pathway prediction method according to claim 1, characterized in that, The quantum chemical method with the first level of computational accuracy is a semi-empirical quantum chemical method.

3. The chemical reaction pathway prediction method according to claim 2, characterized in that, The semi-empirical quantum chemical method is the GFN2-xTB method.

4. The chemical reaction pathway prediction method according to claim 1, characterized in that, The quantum chemical method for the second level of computational accuracy is either the density functional theory method or the post-Hartley-Fock method.

5. The chemical reaction pathway prediction method according to claim 1, characterized in that, The initial reaction path is determined in step (C) using the growth string method.

6. The chemical reaction pathway prediction method according to claim 1, characterized in that, The Bernie optimization algorithm is used to optimize the initial transition state in step (D).

7. The chemical reaction pathway prediction method according to claim 1, characterized in that, The determination of reaction pairs in step (A) is achieved by applying a preset chemical bond recombination rule to the molecular graph structure of the reactant (10) to automatically enumerate and generate the molecular graph structure of the potential product (20).

8. The chemical reaction pathway prediction method according to claim 1, characterized in that, Also includes: (E) Integrate the precise transition state (30) and its corresponding reaction path information into a reaction network (40); (F) When there are multiple intermediates in the reaction network (40) that can be used as new reactants, the exploration priority of each intermediate is determined based on the highest activation energy barrier in the path connecting the initial reactant to each intermediate, and the intermediate with the highest priority is selected as the new reactant (10), and steps (A) to (E) are executed iteratively.

9. A chemical reaction pathway prediction system, characterized in that, include: At least one processor; The system also includes a memory communicatively connected to the processor, the memory storing instructions executable by the processor, which, when executed, cause the system to implement a chemical reaction path prediction method as described in any one of claims 1-8.

Images