Transition state and catalytic network construction method

Abstract

The invention discloses a transition state and catalysis network construction method, which comprises the following steps of: realizing automatic enumeration and screening of reaction intermediates by integrating an element reaction rule base and a graph theory algorithm; intelligent matching and weight optimization of the DFT functional are realized by fusing quantum calculation and machine learning features; by establishing a differentiated transition state search strategy of a multi-complexity scene, high-precision energy barrier calculation from a simple system to a complex system is realized; and dynamic alignment of the catalytic network and real chemical behaviors is realized through path calibration driven by experimental data and an iterative convergence mechanism, and finally, a comprehensive, self-consistent and verifiable catalytic reaction network is constructed.

Description

Technical Field

[0001] This invention relates to the field of metal catalysis research technology, and more specifically, to a method for constructing transition states and catalytic networks. Background Technology

[0002] In organic synthesis and drug development, palladium-catalyzed Heck reactions are a key method for constructing carbon-carbon bonds and are widely used in the large-scale production of fine chemicals and drug molecules. Palladium-catalyzed Heck reactions achieve cross-coupling of aryl halides with alkenes via palladium catalysts, and the reaction pathway typically involves elementary steps such as oxidative addition, migration insertion, and reductive elimination. Current palladium-catalyzed Heck reaction systems often rely on manual modeling and calculations by researchers. For example, when attempting to predict the Heck reaction pathway of 4-bromostyrene and ethyl acrylate catalyzed by a Pd / dppf binuclear catalyst, researchers need to manually construct possible intermediate structures. This process heavily relies on personal experience and often overlooks crucial intermediates such as ligand bridging and metal-metal bond formation, resulting in an incomplete reaction network. Furthermore, the stable coordination number of palladium centers varies between 2 and 6 under different ligand environments, and current methods struggle to systematically enumerate all reasonable structures that meet the coordination number requirements.

[0003] Furthermore, in the transition state search phase, researchers need to select and perform calculations using density functional theory methods individually for each possible reaction step. For example, common functionals such as B3LYP and M06-2X perform well in handling mononuclear palladium complexes, but when dealing with metal-metal interactions in binuclear palladium systems, the calculation results of different functionals may produce errors exceeding 10 kJ / mol. This uncertainty makes the energy barrier prediction unreliable, directly affecting the determination of the rate-determining step of the reaction.

[0004] In addition, when the calculated energy barrier deviates from the gas-phase experimental data to a certain extent, researchers often need to readjust the calculation parameters or change the theoretical method for calibration. This trial-and-error process is time-consuming and inefficient. This deviation is particularly noticeable in reaction systems with significant solvation effects (such as the Heck reaction in DMF solvent). Summary of the Invention

[0005] To address the challenges in constructing catalytic reaction networks, this invention provides a method for constructing transition states and catalytic networks. By integrating a primitive reaction rule base with graph theory algorithms, it achieves automated enumeration and screening of reaction intermediates. By fusing quantum computing and machine learning features, it enables intelligent matching and weight optimization of DFT functionals. By establishing differentiated transition state search strategies for multi-complexity scenarios, it achieves high-precision energy barrier calculations for systems ranging from simple to complex. Through experimental data-driven path calibration and iterative convergence mechanisms, it achieves dynamic alignment between the catalytic network and real chemical behavior, ultimately constructing a comprehensive and verifiable catalytic reaction network.

[0006] The technical solution of this invention is as follows: A method for constructing transition states and catalytic networks includes the following steps: Step S1. Obtain the basic data of the target catalytic system uploaded by the user, load the built-in database of elementary reaction rules containing multiple bond breaking and formation modes, parse the data and extract microstructure information, preprocess the data and unify the data format, establish a dedicated task ID and parameter index table, generate logs that record initial parameters, start time and system resource allocation and back them up to local and cloud simultaneously. Step S2. Use graph theory algorithm to traverse the bond interaction combinations between catalyst and substrate, combine elementary reaction rules to screen breakable bond pairs to enumerate reaction intermediates, calculate the relative energy and bond change cost of intermediates, and screen and remove duplicates of enumerated intermediates based on preset energy threshold, metal stable coordination number range, atomic spatial steric hindrance conditions and electron cloud overlap principle to obtain a set of screened intermediates.

[0007] Step S3. Based on quantum chemical calculation and machine learning models, extract the local metal features of the intermediate and the global features of the substrate, quantize the feature data to the 0-1 range to obtain a standardized vector, remove abnormal samples, and obtain the feature dataset in which the samples and intermediates have a one-to-one mapping relationship. Step S4. Call the classification prediction model based on the random forest algorithm, divide the intermediates into 9 basic categories according to the relevant information of the number of transition metal nuclei and the tooth size of ligands, evaluate the error and efficiency of several DFT functionals in the energy calculation of various intermediates, select the optimal functional, calculate and label the weight coefficient of each functional according to the reciprocal of the error, and construct a density functional approximate consensus set. Step S5. Based on the DFT approximation error value and the number of intermediate atoms, intermediates are divided into three complexity scenarios: low, medium, and high. Differentiated transition state generation methods are adopted so that intermediates in the three different complexity scenarios obtain transition states using different methods. For the low complexity scenario, the initial transition state is generated by interpolation, the structure is optimized by calling the functional with the highest weight, and the imaginary frequency characteristics are verified by frequency analysis to determine the qualified transition state. For the medium complexity scenario, candidate structures are generated by molecular dynamics simulation, and high-quality structures are obtained through energy calculation and structure screening and iterative optimization. Finally, the low-energy structure that meets the imaginary frequency condition is selected as the optimal transition state. For the high complexity scenario, functional cross-validation is added on the basis of the medium complexity method. Abnormal structures are eliminated through multifunctional parallel computation and variance analysis, and the transition state is selected after iteration. Step S6. First, generate a complete path based on the structural matching between the intermediate and the transition state in step S5. Use the atomic coordinate RMSD value as the matching criterion to form a continuous sequence from reactants to products and record energy change data. Then, compare the original paths pairwise. First, filter and eliminate redundant paths according to the energy barrier difference threshold of key steps. Then, retain different types of paths by judging the mechanism difference through the bond change mode sequence to obtain the core path. Finally, call the experimental database to match the data at a set frequency, classify the deviation level between the path and the experimental data, and take differentiated iterative optimization to form a closed loop. Step S7. Iteratively converge the reaction paths that have undergone multiple calibration processes in Step S6. Integrate the reaction paths with smaller deviations in multiple consecutive path calibrations that meet the rated threshold into a catalytic reaction network and label them. Calculate the accuracy of the DFT approximate dynamic selection of the catalytic reaction network and the median of the overall reaction barrier calculation error. If both the accuracy of the DFT approximate dynamic selection of the catalytic reaction network and the median of the overall reaction barrier calculation error are satisfied, output the catalytic reaction network and send a notification report. If neither of these thresholds is satisfied, return to the path calibration in Step S6 for optimization and adjustment.

[0008] In the aforementioned transition state and catalytic network construction method, the core data received in step S1 includes a CIF format crystal structure file uploaded by the user that records the space group and three-dimensional coordinates of the catalyst crystal, an XYZ format file that records the molecular structure, and a substrate SMILES format string that expresses the molecular topology. Material property data must include at least the substrate purity parameter to identify and remove impurities whose content is below a set threshold. The reaction condition parameters cover the temperature range, pressure threshold and solvent type, and are transformed into calculation boundary parameters through density functional theory. Data preprocessing includes using a built-in parser to extract atomic coordinates, atomic connections, and electronic states from the structure file, correcting minor deviations in atomic coordinates caused by measurement errors, and converting all data into a unified format that the system can recognize. Establishing the associated index involves building a parameter index table to associate rule IDs, intermediate pre-numbers, and computing node information, and generating logs that record initial parameters, task start times, and system resource allocations, which are then synchronously backed up to local and cloud storage.

[0009] Furthermore, the general rule set includes state transition patterns and constraints verified in the field, a database of elementary reaction rules containing more than 120 typical catalytic reaction bond patterns, and can automatically select a subset of rules specific to the research object.

[0010] In the above-mentioned transition state and catalytic network construction method, in step S2, when constructing the molecular graph model, atoms are labeled with atom-related information as nodes with attributes, and chemical bonds are labeled with chemical bond-related information as edges with parameters. When screening for breakable bond pairs, first screen the bonds that are directly connected to the metal-ligand of the catalyst and the functional group characteristic bonds of the substrate from the molecular diagram, and then pair them according to the combination forms in the rules, retaining only the combinations whose bond energies are within the allowable range of the rules. The intermediate's unique ID format is reaction type-metal type-serial number, with the serial number increasing in the order of generation. The intermediate output is an XYZ format molecular structure file. The first line of the file indicates the total number of atoms, the second line indicates the unique ID, and each subsequent line records the atom symbol and the X, Y, and Z coordinates in sequence. The energy threshold is set according to the system type, and the energy cost of bond change is calculated by subtracting the total bond energy of newly formed chemical bonds from the total bond energy of breaking all chemical bonds.

[0011] In the above-mentioned transition state and catalytic network construction method, during the metal stable coordination number check in step S2, each metal is judged according to its inherent stable coordination number range, and those exceeding the range are judged as structurally abnormal and eliminated. The shortest distance between atoms is calculated as the distance between the ligand atom and the substrate atom. If this distance is less than the sum of the van der Waals radii of the two atoms, the steric hindrance is deemed too great and the atom is rejected. The bond formation rationality check is based on the principle of electron cloud overlap, and bond formation structures that violate this principle are directly eliminated; Molecular fingerprint similarity calculation is used for deduplication. When the molecular fingerprint similarity between two intermediates is greater than or equal to a set threshold, the individual with lower energy is retained and duplicate individuals are removed.

[0012] In the above-mentioned transition state and catalytic network construction method, in step S3, the local metal features are obtained by density functional theory calculations and include at least a number of indicators, including the electron cloud distribution of the metal nucleus, the electronegativity contribution of the coordinating atoms, the bond order and bond energy of the metal-ligand bond, the number and distribution of empty orbitals of the metal atoms, the bond angle of the metal-ligand bond, and the charge distribution of the coordinating atoms. The global characteristics of the substrate are extracted using molecular topology analysis tools and include at least several indicators such as the topological index of the molecule, the type and number of functional groups, the uniformity of charge distribution, the molecular dipole moment, the bond length and bond angle distribution, the molar mass and volume, and the length of the conjugated system. If the proportion of missing features in a certain intermediate exceeds a preset threshold, the intermediate is identified as an abnormal sample and removed. For the few missing features in the retained samples, the average feature value of the same type of intermediate is used to fill them in, and finally a high-dimensional feature dataset is formed that is associated with the intermediate ID one by one.

[0013] Furthermore, there are a total of 18 key indicators for local metal features and 24 key indicators for global substrate features, which together constitute a high-dimensional feature dataset with 42 feature indicators; the preset threshold for missing features is 5%, that is, features with ≥3 missing features are removed.

[0014] In the aforementioned transition state and catalytic network construction method, in step S4, the classification prediction model is constructed based on the random forest algorithm, and its training data comes from a large amount of intermediate feature data of the catalytic system; Metal core counts are classified into three categories: single-core, dual-core, and multi-core. Ligand dentation is classified into three types: monodentate, dipaldentate, and multidentate. Nine basic categories are formed by combining two dimensions: the number of metal nuclei and the dentation of ligands. For intermediates containing mixed dentation ligands, they are classified into subcategories of the corresponding basic categories based on their dominant structural features. For each type of intermediate, we evaluated a variety of functionals, including B3LYP, B3LYP-D3, PBE, PBE0, M06-2X, TPSS, B97X-D, CAM-B3LYP, B97D3, and M11. The evaluation metrics included the average error of energy calculations with reference to high-precision quantum chemical calculation methods and the average computation time under the same hardware configuration. The optimal functional is selected by comprehensively considering error and efficiency, ranking and scoring the functionals, and then choosing the one with the highest score. Constructing a density functional approximate consensus set involves selecting 3-5 optimal functionals for each type of intermediate and assigning a weight coefficient to each functional based on the reciprocal of its error. The calculation formula is: weight coefficient of a functional = (1 / error of the functional) / (1 / sum of errors of all selected functionals).

[0015] Furthermore, taking the Pd-catalyzed Heck reaction as an example, the classification prediction model divides the intermediates into mononuclear-monodentate, binuclear-bidentate, and binuclear-bidentate-mixed ligands; For mononuclear-monodental intermediates, the optimal functionals selected include B3LYP-D3, M06-2X, and PBE, with weighting coefficients of approximately 0.35, 0.40, and 0.25, respectively.

[0016] In the aforementioned method for constructing a transition state and catalytic network, the criteria for scene division in step S5 are as follows: The low-complexity scenario is one where the DFT approximation error is ≤5kJ / mol and the number of intermediate atoms is ≤50; Medium-complexity scenarios are characterized by DFT approximation errors between 5 and 15 kJ / mol or intermediate atom counts between 50 and 100. High-complexity scenarios are those where the DFT approximation error is >15 kJ / mol or the number of intermediate atoms is >100; For low-complexity scenarios, an initial transition state is generated by interpolating the coordinates of reactant and product structures. The density functional approximates the functional with the highest weight in the consensus set for structural optimization. Frequency analysis is used to verify whether there is a unique imaginary frequency with the same direction as the reaction coordinates to determine a qualified transition state. For medium-complexity scenarios, candidate structures are generated through molecular dynamics simulations, their combined energies are calculated, and high-quality structures are selected for iterative optimization. Finally, low-energy structures that meet the imaginary frequency condition are selected as the optimal transition state through frequency analysis. For high-complexity scenarios, a functional cross-validation mechanism is added to the medium-complexity method. That is, the energy of candidate structures is calculated in parallel using all functionals in the consensus set, and abnormal structures with unstable energy data are eliminated through variance analysis. After iteration and frequency analysis, a highly reliable transition state is output.

[0017] Furthermore, in the transition state generation method for low-complexity scenarios, the corresponding atoms in the reactant and product structures are identified, a weight ratio from 0 to 1 is set, the average value of the corresponding atom coordinates is calculated to generate a series of intermediate structures, and the structure with energy between reactant and product is selected as the initial transition state. Structural optimization involves calculating the forces acting on atoms and iteratively adjusting atomic coordinates until the structural energy reaches its peak in the reaction coordinate direction and its minimum in other directions, and all bond lengths and bond angles fall within the chemically reasonable range determined based on authoritative structural databases. In frequency analysis, the criterion for a qualified transition state is that there is only one imaginary frequency, and the direction of the atomic vibration corresponding to the imaginary frequency is consistent with the stretching direction of the bond to be broken and the contraction direction of the bond to be formed in the reaction.

[0018] Furthermore, the reasonable chemical range for bond length refers to the range of ±3 standard deviations of the average value of the corresponding bond type in the authoritative structure database, and is relaxed to ±5 standard deviations for bonds to be broken / formed in the transition state.

[0019] Furthermore, in the transition state generation method for medium-complexity scenarios, molecular dynamics simulation needs to configure simulation temperature and time step parameters to generate candidate structure snapshots with significant structural differences. The criterion for judging structural differences is that the maximum value of the coordinate deviation of corresponding atoms in any two snapshots is greater than or equal to a set threshold. The calculation of the overall energy involves calling all functionals in the density functional approximation consensus set to calculate their energy separately, and then weighting and fusing them according to the weight coefficients of each functional. Screening for high-quality structures involves checking whether the metal coordination number is within the characteristic range, whether the non-bonding atomic spacing is greater than or equal to the sum of van der Waals radii, and retaining structures with the highest overall energy ranking. The iteration stops when the total energy difference between the retained structures in two consecutive rounds is less than or equal to a set threshold.

[0020] Furthermore, the parameter configuration standard for molecular dynamics simulation is to ensure that the simulation temperature can generate appropriate atomic vibrations without destroying the molecular framework, and the time step matches the atomic vibration period. In the screening of high-quality structures, the characteristic coordination number range for Pd intermediates is typically 4.

[0021] Furthermore, in the transition state generation method for high-complexity scenarios, the functional cross-validation mechanism is as follows: allocate independent computing resources to each functional in the density functional approximation consensus set, and simultaneously perform energy calculations on each high-quality candidate structure. After the calculation is completed, the energy data of the same structure under different functionals are extracted and the variance is calculated. If the variance result exceeds the set variance threshold, the energy data of the structure is determined to be unstable and is removed. The structure that meets the variance requirements and has the best overall energy is retained for subsequent iterations and frequency analysis.

[0022] In the above-mentioned transition state and catalytic network construction method, in step S6, a certain set value of the root mean square deviation of atomic coordinates is used as a judgment condition to splice the product of the preceding structure and the reactant of the following structure in the intermediate and the transition state, thereby forming a continuous sequence from reactant to product, and recording the relative energy and energy barrier data of each step. The selection of core paths includes: First, comparing the original paths pairwise. If the difference in energy barriers of the key steps of the two paths is less than a set threshold and the energy barriers of other steps are the same as the mechanism, they are considered redundant paths, and only the one with the lower energy is retained. Second, comparing the bond change pattern sequence of the retained paths. If the bond breaking type, bond formation type or bond position of any step is different, they are determined to be paths with different mechanisms and are all retained. In path calibration, the differences between the generated path and the experimental data in terms of the average deviation of the structural atomic coordinates and the energy barrier deviation are calculated, and the deviations are divided into three levels: small, medium and large. Differentiated iterative optimization strategies are adopted for different deviation levels, such as updating feature weights, increasing transition state sampling, replacing DFT functionals or retraining classification prediction models. The calibration frequency is set to be triggered after a certain number of core paths are generated.

[0023] Furthermore, the differentiated iterative optimization strategy is as follows: For paths with smaller deviations, update the cost function weights to prioritize their corresponding features; For paths with moderate deviation, if the deviation is structural, increase the number of transition state sampling rounds; if the deviation is energy, replace the functional in the density functional approximation consensus set. For paths with large deviations, the classification prediction model is retrained, and the relevant data is used as new samples to update the model parameters.

[0024] In the aforementioned transition state and catalytic network construction method, the criterion for determining iterative convergence in step S7 is: in three consecutive path calibrations, the proportion of paths with smaller deviations is greater than or equal to 80%. Using intermediates as nodes and reaction steps as edges, nodes are labeled with ID and relative energy and color-coded to distinguish energy levels, while edges are labeled with reaction type, energy barrier value and bond change mode, forming a complete network link from reactants to products. The performance evaluation includes at least the accuracy of the DFT approximate dynamic selection and the median of the overall reaction barrier calculation error. If either indicator fails to meet the standard, the process is backtracked to the path calibration step in step S6 for optimization.

[0025] Furthermore, the formula for calculating the accuracy of the DFT approximate dynamic selection is: (number of paths matching experimental data / total number of core paths) × 100%, where the matching criteria are structural RMSD ≤ 0.1 nm and energy barrier deviation ≤ 5 kJ / mol; The visualization of catalytic reaction networks supports rotating and zooming to view molecular configurations and can dynamically demonstrate structural changes in the order of reactions.

[0026] In the above-mentioned transition state and catalytic network construction method, in step S7, the multi-format output includes a visualization file that supports three-dimensional structure display and path roaming, a data table containing all core parameters, and a PDF report that records the calculation process in detail; The task completion notification must include the task ID, total execution time, core performance metrics, and file storage path information, and be sent to the user in a specified manner.

[0027] According to the above-described solution, the beneficial effects of this invention are as follows: 1. This invention, by constructing an automatic enumeration logic based on a chemical knowledge base, can effectively overcome the problem of intermediate omissions caused by the limitations of personal experience in manual modeling. It systematically generates intermediate candidate structures covering multiple metal coordination modes, oxidation states, and geometric configurations, thereby significantly improving the completeness of catalytic networks in chemical space exploration and providing a solid structural foundation for revealing unconventional but potentially existing reaction channels.

[0028] 2. By establishing an automatic matching mechanism of theoretical methods associated with the properties of the reaction center and the reaction type, this invention can select more suitable computational functionals and basis sets for transition state calculations with different electronic structures and bonding characteristics, thereby significantly improving the calculation accuracy of key energy barriers and relative energies. In particular, when dealing with systems involving metal-metal bonds, weak interactions and open-shell systems, the reliability of the calculation results and the consistency with chemical intuition are enhanced.

[0029] 3. This invention transforms experimentally observed byproducts, selectivity trends, and kinetic data into computable constraints and feeds them back into the path generation process. This enables continuous alignment and dynamic correction between the computational network and experimental evidence, allowing the theoretical model to not only explain known phenomena but also predict experimentally possible reaction paths that may not be detected, thereby strengthening the empirical basis of the computational conclusions.

[0030] 4. This invention develops an automated path connection algorithm based on structural fingerprints and energy gradients, which can autonomously identify reaction sequences that satisfy energy continuity and structural evolution rationality from massive amounts of intermediate and transition state data, and automatically eliminate false paths that do not satisfy the basic laws of thermodynamics or have structural conflicts, thereby constructing a logically consistent and mechanistic core reaction network map.

[0031] 5. This invention integrates a multi-criteria screening model that comprehensively considers the stability of intermediates, accessibility of transition states, and microscopic reversibility. While preserving the diversity of mechanisms, it can intelligently focus on the dominant reaction channels with the greatest kinetic advantages. This results in a catalytic network that avoids distortion due to oversimplification and clearly points to the path most likely to be followed in actual catalytic cycles. Detailed Implementation

[0032] To make the technical problems, technical solutions, and beneficial effects of this invention clearer, the invention will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0033] A method for constructing transition states and catalytic networks is shown below.

[0034] Step S1. Obtain the basic data of the target catalytic system uploaded by the user, load the built-in database of elementary reaction rules containing multiple bond breaking and formation modes, parse the data and extract microstructure information, preprocess the data and unify the data format, establish a dedicated task ID and parameter index table, generate logs that record initial parameters, start time and system resource allocation and back them up to local and cloud simultaneously.

[0035] The data acquired for the calculations includes user-uploaded structure files, including CIF format crystal structure files (recording information such as crystal space group and atomic three-dimensional coordinates) and XYZ format molecular structure files. These files record the relative positions of atoms in the molecule and also include SMILES format strings for the substrate, which use concise symbols to express the molecular topology. These files are the core basis for reconstructing the microstructure of matter and directly determine the accuracy of subsequent structure analysis.

[0036] The data required for the calculation also includes material property data, with the core being the purity parameter of the substrate. This is because impurities can interfere with the determination of reaction pathways and may lead to deviations in the calculation of side reaction pathways. Therefore, it is necessary to ask users to provide purity data and remove impurity data with contents below a certain value to ensure the purity of the substrate data.

[0037] The data required for the calculations also includes reaction condition parameters, covering temperature ranges, pressure thresholds, and solvent types. These parameters need to be transformed into computational boundaries using density functional theory (DFT).

[0038] The acquired data is preprocessed to remove data biases and standardize the format. A built-in CIF / XYZ parser extracts key information from the structure file, including the catalyst's atomic coordinates, atomic connections, and electronic states. The SMILES parsing tool extracts the functional group composition and topology from the substrate string, while purity parameters are used to identify impurities, preparing for subsequent removal operations.

[0039] Minor deviations in atomic coordinates caused by instrument measurement errors in the structural file were corrected. Although small, these deviations can affect the accuracy of bond length calculations and must be corrected. Impurity data with a concentration below 0.1% in the substrate were removed to avoid interference from impurities on the reaction pathway. Then, all data were converted to a unified format recognizable by the system to facilitate subsequent calculations and quantitative comparisons. (The corrected reference data comes from public databases such as the Cambridge Structural Database (CSD) and the Protein Database (PDB), or data from these public databases can be directly imported.) A data management and traceability mechanism was established, assigning a unique ID to this task and constructing a parameter index table. This table associates the parameter index table with information such as rule IDs, intermediate pre-numbers, and compute node IPs, ensuring that subsequent stages can quickly and accurately access the corresponding data. A startup log was generated for this task, detailing the initial parameters, task start time, and system resource allocation. This log was also backed up to a local SSD and a cloud NAS storage system, facilitating troubleshooting, ensuring no data loss, and providing a basis for result reproduction.

[0040] In this invention, the Pd-catalyzed Heck reaction is used as the object, and basic data uploaded by users is obtained. The catalysts are mononuclear Pd(0)-PPh3 complexes (CIF format, crystal space group P-1) and binuclear Pd2-dppf complexes (XYZ format, Pd-Pd bond length in the original data is 2.68 Å, including 0.02 Å coordinate bias); the substrates are 4-bromostyrene (SMILES: C1=CC=C(C1)Br, purity 99.8%, containing 0.05% 4-chlorostyrene impurity) and ethyl acrylate (SMILES: CCOC(=O)C=CH2, purity 99.5%, no obvious impurities); the reaction conditions are temperature 298-323 K, pressure 1 atm, and solvent DMF (dimethylformamide). The built-in elementary reaction rule database is loaded synchronously. This database contains more than 120 typical catalytic reaction bond modes and automatically filters out the Heck reaction-specific rule subset. The core of this Heck reaction-specific rule subset is the Pd-catalyzed double bond breaking and forming double bond mode, such as breaking Pd-P+C-Br to Pd-C+Pd-Br, and breaking Pd-C+C=C to C-C+Pd-C.

[0041] Key information was extracted from the uploaded data, including the coordinates of the single-core Pd atom (0.21, 0.35, 0.18), coordination number 4, and d¹ using a CIF parser. 0Hybridization state. The d orbital occupancy of the binuclear Pd atom (9.3) and the bidentate coordination mode of the dppf ligand with Pd were extracted using the XYZ resolver. The bromoaryl-vinyl conjugated structure of 4-bromostyrene, the ester group (-COOC2H5) of ethyl acrylate, and the oxygen atom charge (-0.32e) were identified using the SMILES resolution tool. The 4-chlorostyrene impurity in 4-bromostyrene was also labeled. The reaction conditions were converted into computational boundaries using the B3LYP-D3 functional of DFT. 298K corresponds to a vibrational entropy contribution of 0.0025 eV / K, DMF corresponds to a dielectric constant of 36.7, and 1 atm is set as a fixed unit cell volume.

[0042] The coordinate deviation of the Pd atom in the binuclear Pd2-dppf structure was corrected to 0.02 Å, ensuring the Pd-Pd bond length is precisely maintained at 2.68 Å. Data for the 0.05% 4-chlorostyrene impurity in 4-bromostyrene was removed. All structural data were uniformly converted to a three-dimensional array format of atom type-coordinate-charge, such as Pd atom recorded as Pd,0.21,0.35,0.18,0, and C atom recorded as C,0.35,0.42,0.25,-0.12, with parameter precision retained to 10. -4 Å.

[0043] The task is assigned a unique ID of Pd-Heck-20251117-001. Construct a parameter index table, which is associated with the rule with ID Heck-Pd-003, intermediates pre-numbered from Heck-Int-001 to Heck-Int-050, and compute nodes with IPs such as 192.AA.BB.102.

[0044] Generate a startup log, recording the initial parameters: Pd ligand is PPh3 / dppf, substrate purity is 99.5%+, reaction temperature is 298K, startup time is 2025-11-17 09:23:15, and system resource allocation information: 16-core CPU, 64GB memory, 2 NVIDIA A100 GPUs. Back up the log to the local SSD at D: / Catalysis / Logs / and the cloud NAS at / / NAS / TaskLogs / 20251117 / .

[0045] Step S2. Use graph theory algorithms to traverse the bond interactions between the catalyst and substrate, and use elementary reaction rules to screen for breakable bond pairs; enumerate reaction intermediates and label their attributes; calculate the relative energy and bond change cost, and remove individuals with excessive energy according to the energy threshold; check the stable coordination number of the metal and remove structures with abnormal coordination numbers; calculate the shortest distance between atoms and compare it with the sum of van der Waals radii, and remove those with excessive steric hindrance if the distance is less than the sum of van der Waals radii; check whether the bonding conforms to the principle of electron cloud overlap, and remove those that do not conform; finally, remove duplicates to obtain the set of screened intermediates.

[0046] The intermediate enumeration module is activated. This module uses graph theory algorithms as its core, combined with elementary reaction rules (obtained from an elementary reaction rule database), to traverse bond combination pairs. During this process, the structure files of the catalyst and substrate are first converted into molecular graph models. Atoms are treated as nodes with attributes, labeled with atomic number, charge, hybridization state, etc. Chemical bonds are treated as edges with parameters, labeled with bond order, bond energy, bond length, etc. Then, an adjacency matrix is ​​used to precisely record the connections between atoms, using 1 to represent bonding and 0 to represent no bonding, thus transforming the molecular structure into a computer-computable mathematical model.

[0047] During the process, breakable bond pairs are screened according to the rules of elementary reactions. Specifically, based on the bond change patterns corresponding to the reaction type, bonds directly linking the metal and ligand required for the catalyst are first screened from the molecular diagram, and then functional group characteristic bonds are screened from the substrate molecular diagram. The two types of bonds are then paired according to the combination forms specified in the rules. During pairing, the bond energy data of each bond is checked, retaining only combinations with bond energies within the allowed range of the rules, and eliminating invalid combinations with mismatched bond types or bond energies exceeding the range. This forms a candidate bond pair list, avoiding meaningless bond combination calculations and reducing subsequent workload.

[0048] After enumerating the intermediates, each intermediate is assigned a unique ID in the format of reaction type-metal type-serial number, with the numbers increasing sequentially according to the order of generation. The intermediates are output as a generic XYZ format molecular structure file. The first line of the file indicates the total number of atoms, the second line indicates the unique ID, and subsequent lines record the atom symbols and X, Y, and Z coordinates for easy viewing and calculation of the structure. Core attribute information is also recorded, including bonding mode (recorded in the format of central metal-connecting atoms / substituents, specifying which atoms the metal bonds to and the type of substituents), metal coordination number, and ligand type (classified as monodentate, bidentate, or multidentate based on the number of atoms bonded to the metal).

[0049] Calculate the energy-related parameters of the intermediate. Using the total energy of the catalyst and substrate existing alone as a baseline, calculate the relative energy of the intermediate (relative energy of intermediate = intermediate energy - baseline energy). A positive relative energy indicates that the intermediate's energy is higher than the baseline, suggesting thermodynamic instability; a negative relative energy indicates greater stability. Calculate the energy cost of bond changes, which is the sum of the bond energies of breaking all chemical bonds minus the sum of the bond energies of newly formed chemical bonds. The lower the energy cost, the easier it is for the intermediate to form.

[0050] The energy threshold is set according to the catalytic system. For mononuclear metal complex systems, the threshold is set to 20 kJ / mol because there is no intermetallic synergy. For polynuclear systems, the energy requirement is reduced due to metal synergy, so the threshold is reduced to 18 kJ / mol. If the relative energy or bond change cost of the intermediate exceeds the corresponding threshold, it is determined to be thermodynamically difficult to form and the intermediate is eliminated.

[0051] Check the coordination number of the metals. Each metal has a stable range of coordination numbers. If the coordination number of the intermediate is not within the stable range, it is considered a structural anomaly and should be discarded.

[0052] Calculate the shortest distance between the ligand atom and the substrate atom. If this distance is less than the sum of the van der Waals radii of the two atoms, it is determined that there is severe steric hindrance, the atoms repel each other, the structure is unstable, and it is discarded.

[0053] Check the bonding method. If bonding that violates the principle of electron cloud overlap occurs, it is judged as unreasonable and discarded.

[0054] By calculating molecular fingerprint similarity (a molecular fingerprint is an encoding that reflects the characteristics of molecular structure; the higher the similarity, the more similar the structure), with a threshold set at 95%, if the molecular fingerprint similarity between two intermediates is ≥95%, they are judged to be structurally highly repetitive, the one with lower energy is retained, and the duplicate individuals are removed.

[0055] By optimizing the structure of the retained intermediates using molecular mechanics methods, minor bond angle deviations are corrected, reducing the computational complexity of subsequent feature extraction steps, and ensuring that the final retained intermediates not only fully cover possible reaction pathways but are also streamlined and free of redundancy.

[0056] In this invention, for the Pd-catalyzed Heck reaction, the intermediate enumeration module converts the mononuclear Pd-PPh3 complex and 4-bromostyrene into a molecular diagram. In the mononuclear system, the Pd atom nodes of the mononuclear Pd-PPh3 complex are labeled with atomic number 46, charge 0, and sp³ hybridization; the P atom nodes are labeled with atomic number 15, coordination number 3, and charge -0.3; the Pd-P bond edge is labeled with bond order 1.1, bond energy 28.5 kJ / mol, and bond length 2.25 Å; while the C-Br bond edge of 4-bromostyrene is labeled with bond order 1.8, bond energy 31.2 kJ / mol, and bond length 1.91 Å. In the binuclear system, the Pd atom nodes of the binuclear Pd2-dppf complex are labeled with atomic number 46, charge +0.6, and sp² hybridization, while the P atom nodes of the dppf ligand are labeled with atomic number 15 and coordination number 2. The Pd-P bond edge is labeled with bond order 1.0 and bond energy 22.3 kJ / mol. The C-Br bond labeling information of 4-bromostyrene is consistent with that of the mononuclear system.

[0057] According to Heck's rules for oxidative addition, in the mononuclear system, Pd-P bonds with a bond energy of 28.5 kJ / mol are selected, falling within the 25-30 kJ / mol range. C-Br bonds with a bond energy of 31.2 kJ / mol are also selected, falling within the 30-35 kJ / mol range. These two bonds form a candidate pair. Combinations of Pd-P bonds with a bond energy of 35 kJ / mol (out of range) and C-Br bonds with a bond energy of 28 kJ / mol (out of range) are excluded, resulting in a valid candidate pair combination. For the binuclear system, two Pd-P bonds (each with a bond energy of 22.3 kJ / mol) and one C-Br bond (with a bond energy of 31.2 kJ / mol) need to be selected to form a candidate pair, because the reaction mechanism of the binuclear system requires two Pd-P bonds to participate.

[0058] In the mononuclear system, after the Pd-P bond breaks, the Pd atom generates one unpaired electron, and after the C-Br bond breaks, both the C and Br atoms generate one unpaired electron. The unpaired electron of the Pd atom pairs with the C and Br atoms, respectively, to form a Pd-C bond (bond order 1.7, bond energy 35.7 kJ / mol) and a Pd-Br bond (bond order 1.6, bond energy 22.3 kJ / mol). Subsequently, the adjacency matrix of the molecular diagram is updated, changing the connection between Pd and P from 1 to 0, and the connections between Pd and C and Br from 0 to 1.

[0059] In a binary system, after the two Pd-P bonds break, each of the two Pd atoms generates one unpaired electron. These electrons will pair with the C and Br atoms after the C-Br bond breaks. At the same time, a Pd-Pd metallic bond will be formed between the two Pd atoms (the bond order is 0.8 and the bond energy is 14.9 kJ / mol). The adjacency matrix of the molecular diagram will also be updated accordingly, changing the connection between the two Pd-P bonds to 0, and the connections between Pd and C, Br, and Pd and Pd to 1.

[0060] During the intermediate structure calibration stage, parameters are adjusted according to the system characteristics to ensure a reasonable structure. For the mononuclear system, atomic coordinates are calibrated with Pd-C bond lengths of 2.05±0.02 Å and Pd-Br bond lengths of 2.40±0.02 Å. Simultaneously, the benzene ring of the PPh3 ligand is rotated to maintain a reasonable distance from 4-bromostyrene, avoiding spatial overlap, ultimately generating a mononuclear intermediate. For the binuclear system, coordinates are calibrated with Pd-C bond lengths of 2.06±0.02 Å, Pd-Br bond lengths of 2.41±0.02 Å, and Pd-Pd bond lengths of 2.68±0.01 Å. The biphenyl structure of the dppf ligand is rotated to generate a binuclear intermediate. A total of 92 original intermediates are generated in both mononuclear and binuclear systems, with 43 in the mononuclear system and 49 in the binuclear system.

[0061] The ID of the mononuclear intermediate is Heck-Pd-018. The output XYZ file starts with 32 to indicate 32 atoms, the second line indicates the ID, and subsequent lines record the atom symbols and coordinates, such as Pd 0.22 0.36 0.19, C 0.35 0.42 0.25, etc. The structural properties recorded in the file are the bonding mode Pd-C(aryl)-Br-3PPh3, that is, Pd is bonded to aryl C, Br, and 3 PPh3 ligands, with a coordination number of 4 and a ligand type of monodentate. The energy calculation is based on the sum of the energies of the mononuclear catalyst Pd(PPh3)4 (E1), 4-bromostyrene (E2), and ethyl acrylate (E3) as the baseline energy E. 总 The relative energy is equal to the energy of the mononuclear intermediate (E). 中 Subtract E 总 The result is +1.7 kJ / mol, indicating that the intermediate is 1.7 kJ / mol higher than the initial state. The cost of bond change is the total bond energy of breaking old bonds minus the total bond energy of forming new bonds, that is, (28.5+31.2)-(35.7+22.3)=59.7-58.0=+1.7 kJ / mol.

[0062] The ID of the binuclear intermediate is Heck-Pd2-007. The output XYZ file starts with 45 to indicate 45 atoms, the second line indicates the ID, and subsequent lines record the atom symbols and coordinates, such as Pd10.240.380.20, Pd20.500.380.20, etc. The structural properties recorded in the file are: bonding mode Pd2-(dppf)-C(aryl)-Br, where two Pd atoms are connected by a dppf ligand, simultaneously binding an aryl C and Br, with a coordination number of 3-4 (the coordination numbers of the two Pd atoms are 3 and 4 respectively), and the ligand type is bidentate. The baseline energy for energy calculation is E. 总' The energy (E) for the binuclear catalyst Pd2-dppf c The sum of the energies of 4-bromostyrene (E2) and ethyl acrylate (E3) differs from the baseline energy of the mononuclear system; its relative energy is the energy of the binuclear intermediate (E...). 中' Subtract E 总' The result is +0.9 kJ / mol. This is because the Pd-Pd metal bond in the binuclear system has a stabilizing effect, the intermediate energy increases less, and the bond change cost is (2×22.3+31.2)-(38.0+23.0+14.9)=75.8-75.9=+0.9 kJ / mol.

[0063] For mononuclear systems, an intermediate with a relative energy of +28 kJ / mol and a bond change cost of +29.1 kJ / mol was eliminated based on a threshold of 20 kJ / mol. For binuclear systems, due to the synergistic effect between metals, the threshold was set at 18 kJ / mol, and intermediates with a relative energy of +17 kJ / mol and a bond change cost of +16.8 kJ / mol were retained. In the structure screening, three intermediates with a Pd coordination number of 3 and two intermediates with a coordination number of 5 were eliminated due to structural abnormalities. The distance between the C atom of the benzene ring of the PPh3 ligand and the O atom of the ethyl acrylate ester group was calculated. The distance of six intermediates was 2.4 Å, which is less than the C and O van der Waals radii and 3.2 Å, respectively. These intermediates were deemed to have excessive steric hindrance and were therefore eliminated. No unreasonable bonding structures were found, and 56 intermediates were initially retained. Five repetitive intermediates were then eliminated using molecular fingerprint similarity (threshold 95%), leaving 51 core intermediates, including 28 mononuclear intermediates, 15 binuclear intermediates, and 8 ligand-bridged intermediates. The structures were then optimized using molecular mechanics, such as rotating the dppf ligand biphenyl structure to make it form a 30° angle with the Pd-Pd bond and correcting the Pd-C bond angle to 109.5° to ensure a stable and reasonable structure.

[0064] Step S3. Based on quantum chemical calculation and machine learning models, extract the local metal features of the intermediate and the global features of the substrate, quantize the feature data to the 0-1 range to obtain a standardized vector, remove abnormal samples, and obtain the feature dataset in which the samples and intermediates have a one-to-one mapping relationship.

[0065] To transform the structure and properties of intermediates into computer-processable quantitative features, a feature selection algorithm (such as variance threshold calculation) is first initiated. This algorithm combines quantum chemical calculations (to ensure the chemical significance and accuracy of the features) with machine learning models (to ensure the computability and practicality of the features) to extract two types of features from the intermediates. These two types of features correspond to the active core of the catalyst and the overall characteristics of the substrate, respectively.

[0066] The first category of features consists of local metallic characteristics of transition metal complexes. These features are directly related to catalyst activity and reaction mechanisms and are obtained through PBE functional calculations using density functional theory (DFT). A total of 18 key indicators are identified, including the electron cloud distribution of the metal core, the electronegativity contribution of coordinating atoms, the bond order and bond energy of the metal-ligand bond, the number and distribution of empty orbitals of the metal atom, the bond angle of the metal-ligand bond, and the charge distribution of the coordinating atoms. The electron cloud distribution of the metal core affects the substrate bonding ability; the electronegativity contribution of coordinating atoms determines the direction of electron transfer and the polarity of the metal-ligand bond; the bond order and bond energy of the metal-ligand bond reflect bond stability and directly affect the difficulty of reaction initiation; the empty orbitals of the metal atom determine the coordination number and substrate spatial fit; the metal-ligand bond angle affects ligand orientation and substrate binding convenience; and the charge distribution of coordinating atoms reflects electron-donating ability and regulates the catalytic activity of the metal center. These local metallic characteristics represent the state of the catalyst active center and can be used to determine the feasibility and rate of the reaction.

[0067] The second category of features is the global characteristics of organic substrates. These features reflect the overall structure and reaction potential of the substrate and are extracted using molecular topological analysis tools. A total of 24 indicators are identified, including the topological index of the molecule, the type and number of functional groups, charge distribution uniformity, molecular dipole moment, bond length and bond angle distribution, molar mass and volume, and the length of the conjugated system. The topological index characterizes the complexity of the molecular structure and affects its compatibility with catalysts; the type and number of functional groups determine the reaction type and product diversity, representing the core active sites; charge distribution uniformity locates the reactive regions, and inhomogeneities are likely to become the breakthrough points for the reaction; the molecular dipole moment affects solubility and dispersibility in polar solvents; bond length and bond angle distribution reflect the geometric stability of the molecule, providing a reference for bond breaking and recombination; molar mass and volume affect the flexibility of molecular motion and the probability of collision with the catalyst; and the length of the conjugated system enhances electron delocalization, improving reactivity and selectivity.

[0068] After feature extraction, quantization and cleaning are required to ensure the data format is uniform and complete. The quantization process maps all feature values ​​to the 0-1 range, transforming them into standardized numerical vectors. Specifically, the value range of each feature is first determined, and then the standardized vector is calculated using the linear mapping formula: Standardized value = (Actual value - Minimum value) / (Maximum value - Minimum value). This standardization process eliminates the influence of differences in the magnitude of different features, ensuring that all parameters can be compared and calculated on the same dimension.

[0069] The number of missing features for each intermediate is counted. If the proportion of missing features to the total number of features exceeds 5% (i.e., ≥3 missing features), it is considered an outlier and removed. Missing features may stem from structural resolution failures (such as the inability to accurately calculate the electronic state of certain complex ligands) or data transmission errors. Such samples will interfere with subsequent classification results and must be removed.

[0070] After feature extraction, standard quantization, and removal of missing / outliers, a feature dataset is obtained. Each sample contains 42 feature indicators (18 local metal features + 24 global substrate features), and each is associated with an intermediate ID. The feature dataset retains the essential chemical characteristics of the intermediates (such as metal activity and substrate functional groups) and has the standardized format required by machine learning algorithms. Subsequent classification algorithms can analyze these features to identify the structural commonalities of intermediates (such as mononuclear vs. binuclear, monodentate vs. bidentate ligands), and DFT functional selection can determine the structural complexity of the intermediate based on the features (e.g., intermediates with long conjugated system lengths require more precise functionals).

[0071] In this invention, for 51 core intermediates of the Pd-catalyzed Heck reaction, dual-dimensional features of the transition metal complex and the organic substrate were extracted.

[0072] Regarding the local features of the metal, calculations were performed using the PBE functional of the DFT: The electron cloud density of the d orbitals of the Pd atom in the mononuclear Pd intermediate is 0.82, corresponding to a normalized value of (0.82-0.5) / (1.0-0.5) = 0.32 / 0.5 = 0.64; the electronegativity contribution of the P atom is 0.35, corresponding to a normalized value of 0.35, ranging from 0 to 1; the Pd-C bond order is 1.7, corresponding to a normalized value of (1.7-1.0) / (2.0-1.0) = 0.7; the Pd-C bond energy is 35.7 kJ / mol. The corresponding normalized value is approximately 0.52; the number of empty orbitals in a Pd atom is 4, and the corresponding normalized value is (4-2) / (6-2)=2 / 4=0.5; the common number of empty orbitals in Pd is 2-6; the Pd-P bond angle is 109.5°, and the corresponding normalized value is (109.5-90) / (120-90)=19.5 / 30=0.65; the charge of a P atom is -0.3, and the corresponding normalized value is (-0.3+1.0) / (0+1.0)=0.7.

[0073] The electron cloud overlap of the Pd-Pd bond in the binuclear Pd intermediate is 0.68, corresponding to a normalized value of 0.68; the electronegativity contribution of the P atom in the dppf ligand is 0.33, corresponding to a normalized value of 0.33; the Pd-P bond energy is 22.3 kJ / mol, corresponding to a normalized value of (22.3-20) / (50-20)=2.3 / 30≈0.08; the d orbital occupancy of the Pd atom is 9.3, corresponding to a normalized value of (9.3-8.0) / (10.0-8.0)=1.3 / 2=0.65. All local features are quantized into 0-1 interval vectors.

[0074] Regarding global characteristics, extracted using molecular topology analysis tools, the Wiener topological index of 4-bromostyrene is 1.85, corresponding to a standardized value of (1.85-1.0) / (3.0-1.0)=0.85 / 2=0.425; the number of functional groups is 2, corresponding to a standardized value of (2-1) / (5-1)=1 / 4=0.25, assuming a maximum number of functional groups of 5; the charge distribution uniformity is 0.78, corresponding to a standardized value of 0.78, ranging from 0 to 1, with higher values ​​indicating more uniform distribution; the dipole moment is 1.2D, corresponding to a standardized value of 0.12; the C-Br bond length is 1.91Å, corresponding to a standardized value of (1.91-1.8) / (2.0-1.8)=0.11 / 0.2=0.55; and the conjugated system length is 3, with a standardized value of (3-1) / (5-1)=2 / 4=0.5). The Balaban topological index of ethyl acrylate is 1.62, corresponding to a normalized value of (1.62-1.0) / (3.0-1.0)=0.62 / 2=0.31; the number of functional groups is 2, corresponding to a normalized value of 0.25; the charge distribution uniformity is 0.83, corresponding to a normalized value of 0.83; the dipole moment is 1.5D, corresponding to a normalized value of 0.15; and the C=C bond length is 1.34Å, corresponding to a normalized value of (1.34-1.3) / (1.4-1.3)=0.04 / 0.1=0.4). Functional group identification was performed using the SMILES deep learning model, accurately identifying bromoaryl, vinyl, ester, and carbon-carbon double bonds without any false positives.

[0075] The feature loss of 51 intermediates was detected. In this embodiment, one ligand-bridged intermediate was identified as an abnormal sample and removed because its complex structure, containing mixed dentate ligands and an irregular conjugated system, made it impossible to calculate the electronegativity contribution of coordinating atoms, molecular topological index, and charge distribution uniformity. The missing feature rate was 3 / 42 ≈ 7.1% > 5%. The remaining 50 intermediates had ≤ 2 missing features and were retained. For the missing "molecular volume" feature, the average molecular volume of similar intermediates (e.g., the average volume of mononuclear intermediates is 125 ų) was used to fill the gap. The final result is a high-dimensional feature dataset containing 50 samples and 42 feature indicators. The feature vector of each sample is associated with an intermediate ID (such as Heck-Pd-018, Heck-Pd2-007). For example, the feature vector of Heck-Pd-018 is [0.64, 0.35, 0.7, 0.52, 0.5, 0.65, 0.7, ..., 0.425, 0.25, 0.78, 0.12, 0.55, 0.5] (a total of 42 0-1 values), which are used for subsequent intermediate classification.

[0076] Step S4. Call the classification prediction model based on the random forest algorithm, divide the intermediates into 9 basic categories according to the relevant information of the number of transition metal nuclei and the tooth size of ligands, evaluate the error and efficiency of several DFT functionals in the energy calculation of various intermediates, select the optimal functional, calculate and label the weight coefficient of each functional according to the reciprocal of the error, and construct a density functional approximate consensus set.

[0077] First, the intermediates are classified. Then, a pre-trained classification prediction model is called to classify transition metal complexes based on the number of metal nuclei and the tooth size of ligands.

[0078] The classification prediction model is built based on the random forest algorithm. Random forest is an ensemble learning model composed of multiple decision trees. It can handle high-dimensional data and has strong resistance to overfitting. Its training data comes from a large amount of intermediate feature data of catalytic systems. These data include intermediate features with different metal nucleus numbers and ligand tooth sizes, as well as corresponding calculation errors.

[0079] Metal nuclei are classified into three categories based on the actual number of metal atoms they contain: mononuclear, binuclear, and multinuclear. Mononuclear complexes contain one metal atom, such as Pd(0)-PPh3. These complexes have relatively simple structures and are computationally easier. Binuclear complexes contain two metal atoms, such as Pd2-dppf. These complexes have metal-metal bonds, resulting in more complex electronic interactions. Multinuclear complexes contain three or more metal atoms, such as Pd3 clusters. These complexes have extremely complex structures and electronic effects, requiring higher precision calculations.

[0080] Ligand dentation is classified into three categories based on the number of atoms in the ligand molecule that directly bind to a metal atom: monodentate, bidentate, and polydentate. Monodentate refers to a ligand molecule with only one atom directly bound to a metal atom, such as PPh3 ligand which binds to Pd through one P atom. These ligands have low steric hindrance and simple bonding modes. Bidentate refers to a ligand molecule with two atoms directly bound to a metal atom, such as dppf ligand which binds to Pd through two P atoms, forming a chelate ring. This structure has high stability, but ring strain must be considered in calculations. Polydentate refers to a ligand molecule with three or more atoms directly bound to a metal atom, such as EDTA ligand. These ligands have tighter binding to the metal and a more uniform electron distribution, but their structure is more complex.

[0081] The two dimensions combine to form nine basic categories: monodentate-monodentate, monodentate-bidate, monodentate-multidentate, binadentate-monodentate, binadentate-bidate, binadentate-multidentate, multinucleate-monodentate, multinucleate-bidate, and multinucleate-multidentate. For intermediates with special structures, such as complexes containing mixed-dentate ligands, like a Pd atom simultaneously binding one monodentate PPh3 ligand and one bidentate dppf ligand, the classification prediction model automatically identifies its core features. If the bidentate ligand is dominant, it is classified into a subcategory of the corresponding basic category, such as the binadentate-bidate-mixed-ligand subcategory. This ensures that the classification covers general cases without overlooking special structures, avoiding subsequent functional selection bias due to general classification.

[0082] After classification, DFT functional screening is performed to select several optimal functionals for each class of intermediates, forming a density functional approximate consensus set.

[0083] DFT functionals are mathematical expressions used in density functional theory to approximate electron exchange correlation energies. Different functionals have vastly different applicable scenarios, and their error magnitude, adaptability, and computational efficiency vary in different system types. Therefore, it is necessary to evaluate the performance of functionals separately for each type of intermediate to ensure that the most suitable computational method is selected.

[0084] The DFT functional performance evaluation module was activated, integrating computational data from 10 commonly used DFT functionals: B3LYP, B3LYP-D3, PBE, PBE0, M06-2X, TPSS, ωB97X-D, CAM-B3LYP, B97D3, and M11. Among these, some functionals are dispersion-corrected and suitable for weakly interacting systems; some are hybrid functionals with higher accuracy than fundamental functionals; some are elementary functionals, suitable for complex systems or balancing accuracy and efficiency; and some are long-range corrected functionals suitable for highly conjugate systems.

[0085] The evaluation criteria referenced results specifically designed for high-precision quantum chemical calculations (these methods are extremely computationally expensive and are used as standard databases, typically for model validation and calibration). Two metrics were employed in the evaluation analysis. The first metric was energy calculation error. This involved calculating the error of each functional's energy calculation for a given type of intermediate. The calculation method was to subtract the high-precision quantum chemical calculation value from the functional's calculated value. A smaller error indicated more accurate energy calculations by the functional, better reflecting the thermodynamic stability of the intermediate. Specifically, ten DFT functionals to be evaluated were applied to selected typical intermediate samples, and the energy value for each intermediate was calculated individually. The energy calculation results for each sample using each functional were recorded. For each sample, for each functional, the single-term error for each intermediate sample is calculated using the formula: single-term error = calculated value of functional - high-precision reference value. For the same type of intermediate, the average (or median) absolute value of the single-term error of all samples for each functional is calculated as the average energy error of the functional for this type of intermediate, thus eliminating the influence of random errors in individual samples. The standard deviation of the error of each functional is calculated. The smaller the standard deviation, the more stable the functional is in the same type of intermediate, avoiding the risk of subsequent calculations due to large error fluctuations.

[0086] The second metric is computational efficiency. This involves calculating the computation time for each functional on the corresponding intermediate. This calculation must be performed under identical hardware configurations. Shorter computation times indicate higher efficiency and reduce overall task time. Specifically, during the computation of each intermediate sample for each functional, the complete time from computation start to output result is recorded synchronously. Non-computational time such as data read / write operations is excluded; only the time consumed by the core computation process is recorded to ensure the accuracy of the time data. The average computation time for all samples of the same intermediate type for each functional is taken as the average computation time for that functional on that type of intermediate. Using the shortest time among all functionals as a benchmark, the average time for each functional is converted into a normalized efficiency value, i.e., shortest time / average time. The closer the efficiency value is to 1, the higher the computational efficiency.

[0087] The evaluation must balance error and efficiency. For a certain type of intermediate, all functionals are sorted by average energy error from smallest to largest, and assigned scores decreasing by 1. The smaller the error, the higher the ranking and the more points it receives. Similarly, the average time of each functional for this type of intermediate is sorted from shortest to largest (or normalized efficiency value from largest to smallest), and assigned scores decreasing by 1. The smaller the error, the higher the ranking and the more points it receives. In this way, a score can be obtained for each functional for this type of intermediate, and the functional with the largest score is selected as the optimal functional.

[0088] If a functional has extremely small error but takes too long to process, and another functional has slightly larger error but takes very little time to process, then the choice should be based on the importance of the intermediate. Core intermediates should prioritize functionals with smaller errors, while secondary intermediates can be selected with slightly larger errors while still considering efficiency.

[0089] Based on the evaluation results, 3-5 optimal functionals are selected for each type of intermediate to form a density functional approximate consensus set. The density functional approximate consensus set includes the selected optimal DFT functionals and their corresponding weight coefficients. The consensus set for intermediates with relatively simple structures focuses on functionals with small errors and high efficiency, while the consensus set for intermediates with complex structures or special bonding effects focuses on high-precision functionals that fit their structural characteristics. Simultaneously, weight coefficients are assigned to each functional in the consensus set. The weights are calculated based on the reciprocal of the functional's error. For a given type of intermediate, the formula for calculating the weight coefficient of each functional is: Weight coefficient of a functional = (1 / error of that functional) / (1 / sum of errors of all functionals). The smaller the error, the larger 1 / error, the higher the weight, and the greater the proportion of that functional's result in subsequent calculations, effectively reducing the limitations of a single functional. Weight calculation requires first obtaining the error value of each functional, calculating the reciprocal of each functional's error, and then using the sum of the reciprocals of the errors of all functionals as the denominator and the reciprocal of the error of a single functional as the numerator to obtain the weight coefficient of that functional. In subsequent transition state energy calculations, the calculation results are weighted and fused according to the weight coefficients of each functional to improve the accuracy of the calculation.

[0090] In this invention, for 50 intermediates of the Pd-catalyzed Heck reaction, the classification prediction model is used to classify them according to the number of metal nuclei and the dentency of the ligands: 28 intermediates containing 1 Pd atom + monodentate PPh3 ligand are classified into the mononuclear-monodentate basic category; 15 intermediates containing 2 Pd atoms + bidentate dppf ligand are classified into the binuclear-bidentate basic category; and 7 intermediates containing 2 Pd atoms + bidentate dppf + monodentate auxiliary ligand are classified into the binuclear-bidentate-mixed ligand category due to the dominant structure of the bidentate ligand.

[0091] In the DFT functional evaluation phase, for single-core-single-tooth intermediates, under a hardware configuration of 16 CPU cores and 1 NVIDIA A100 GPU, the average energy error (based on CCSD(T), i.e., coupled cluster theory, the result is used as a reference) and average computation time of 10 functionals are calculated, as shown below.

[0092] B3LYP-D3: The average energy error is 3.2 kJ / mol, and the average calculation time is 1.5 hours; M06-2X: The average energy error is 2.8 kJ / mol, and the average calculation time is 2.2 hours; PBE: The average energy error is 4.1 kJ / mol, and the average computation time is 0.8 hours; PBE0: The average energy error is 3.8 kJ / mol, and the average calculation time is 1.2 hours; TPSS: Average energy error is 4.5 kJ / mol, and average computation time is 1.0 hour; The other five functionals (such as ωB97X-D) had average energy errors exceeding 5 kJ / mol or average calculation times exceeding 3 hours, and were therefore excluded.

[0093] Considering both average energy error and efficiency, three functionals—B3LYP-D3, M06-2X, and PBE—were selected to form the consensus set. Weighting coefficients were calculated based on the reciprocal of the average energy error. The weighting coefficient of B3LYP-D3 is (1 / 3.2) / (1 / 3.2+1 / 2.8+1 / 4.1)≈0.3125 / (0.3125+0.3571+0.2439)≈0.35; The weighting coefficient for M06-2X is approximately 0.3571 / 0.9135, which is approximately 0.40. The weighting coefficient of PBE is approximately 0.2439 / 0.9135, which is approximately 0.25.

[0094] For binuclear-bidorine intermediates, the evaluation found that B3LYP-D3 (average energy error 3.5 kJ / mol, average computation time 1.8 hours), PBE0 (average energy error 3.8 kJ / mol, average computation time 1.5 hours), and TPSS (average energy error 4.0 kJ / mol, average computation time 1.3 hours) had the best computational fit for metal-metal bonds and formed a density functional approximation consensus set with weighting coefficients of 0.36, 0.33, and 0.31, respectively.

[0095] The binuclear-bident-mixed ligand classification includes mixed ligands, and additionally incorporates the M06-2X functional (with an average energy error of 3.1 kJ / mol and an average computation time of 2.5 hours), forming a density functional approximation consensus set for the four functionals. The weighting coefficients for B3LYP-D3 are 0.32, M06-2X is 0.35, PBE0 is 0.23, and TPSS is 0.10, ensuring the accuracy of calculations for intermediates with special structures.

[0096] Step S5. Based on the DFT approximation error value and the number of intermediate atoms, intermediates are divided into three complexity scenarios: low, medium, and high. Differentiated transition state generation methods are adopted so that intermediates in the three different complexity scenarios obtain transition states using different methods. In the low complexity scenario, the initial transition state is generated by interpolation, the structure is optimized by calling the functional with the highest weight, and the imaginary frequency characteristics are verified by frequency analysis to determine the qualified transition state. In the medium complexity scenario, candidate structures are generated by molecular dynamics simulation, and high-quality structures are obtained through energy calculation and structure screening and iterative optimization. Finally, the low-energy structure that meets the imaginary frequency condition is selected as the optimal transition state. In the high complexity scenario, functional cross-validation is added on the basis of the medium complexity method. Abnormal structures are eliminated by multifunctional parallel computation and variance analysis, and high-reliability transition states are selected after iteration.

[0097] Specifically, the target intermediate energy obtained through the high-precision quantum chemical calculation method CCSD(T) is used as the reference energy value. This method, because it can maximally reproduce the complex interactions between electrons, considers the calculation result as the true energy reference, and its own error is negligible (0). It is the core reference for measuring the calculation deviation of all DFT functionals. The energy calculation value of the functional for the intermediate is extracted from the density functional approximation consensus set in step S4. The density functional approximation consensus set includes 3-5 optimal DFT functionals selected from all types of intermediates, as well as the energy calculation results of these 3-5 optimal DFT functionals for each type of intermediate. The energy calculation result of each functional for the intermediate minus the reference energy value of that intermediate equals the one-way error of the functional. The sign of the one-way error only represents the direction of the functional's calculated value relative to the true value; the smaller the absolute value of the one-way error, the higher the calculation accuracy of the functional. Based on this, the weight coefficients corresponding to various functionals for this type of intermediate, recorded in the density functional approximation consensus set, are extracted. The individual error of each functional is multiplied by its corresponding weight coefficient, and then all products are summed to obtain the final total error (i.e., the DFT approximation error value). In actual calculations, the absolute value or algebraic value of the individual error can be used for weighting according to subsequent needs. The core logic of both methods is the same: to make the contribution of high-precision functionals to the total error greater through weight allocation, so that the total error is closer to the true deviation, thereby obtaining the DFT approximation error value.

[0098] Simultaneously, using the number of atoms as a direct quantitative indicator of structural complexity, and referencing the atomic number distribution characteristics of mononuclear, binuclear, and ligand-bridged intermediates, the atomic numbers of various intermediates are divided into three categories with different levels of complexity. In this invention, Complexity scenarios are generated based on the DFT approximation error value and the number of intermediate atoms to comprehensively assess the difficulty of transition state generation. The DFT approximation error value is obtained from the classification prediction model and is expressed in kJ / mol; the number of intermediate atoms is statistically analyzed through molecular structure analysis to reflect the spatial complexity of the system. Scenario generation includes low-complexity, medium-complexity, and high-complexity scenarios. Low-complexity scenarios are defined as DFT approximation errors ≤ 5 kJ / mol and intermediate atom count ≤ 50; medium-complexity scenarios are defined as DFT approximation errors 5-15 kJ / mol or intermediate atom count 50-100; and high-complexity scenarios are defined as DFT approximation errors > 15 kJ / mol or intermediate atom count > 100.

[0099] The computational methods for low-complexity scenarios have small deviations and simple structures, with a single transition state configuration, making them suitable for rapid generation. Medium-complexity scenarios either have large computational deviations or complex structures, and the transition states may have multiple configurations, requiring additional sampling exploration. High-complexity scenarios have extremely complex structures or large computational deviations, with multiple transition state configurations and large differences in stability, requiring rigorous verification.

[0100] After scene segmentation, differentiated transition state generation methods are adopted for different scenes, as detailed below. Transition state generation methods for low-complexity scenarios include: The system receives the 3D structural data of the reactants and products corresponding to the reaction, automatically identifies the atoms of the same elements participating in the reaction in both structures, sets a uniform weight ratio for the reactant and product structures from 0 to 1, calculates the average coordinates of each pair of corresponding atoms based on this ratio, and generates a series of intermediate structures. From this series of intermediate structures, structures with energies between those of the reactants and products are selected as the initial transition state. This ensures that the atomic arrangement of the initial transition state closely resembles the transition state of the actual reaction, reducing the workload of subsequent optimization.

[0101] After receiving the initial transition state structural data, the system calls the highest-weighted functional from the preset density functional approximation consensus set, loading the appropriate basis set, energy, and force convergence threshold, among other basic parameters. It calculates the total energy and forces on each atom of the initial transition state structure. It then determines whether the atomic forces meet the convergence threshold criteria: the absolute value of the force on critical atoms is ≤ the set force threshold, the absolute value of the force on non-critical atoms is ≤ 1.5 times the set force threshold, and the root mean square (RMS) of the overall structural forces is ≤ the set RMS threshold. If the convergence threshold is met, adjustment stops. If not, the coordinates of the corresponding atoms in the initial transition state structure are fine-tuned according to the force direction of each atom. The total energy and forces on each atom of the structure after coordinate adjustment are recalculated, and the atomic force judgment and coordinate adjustment steps are repeated. This continues until the structural energy reaches its peak in the reaction coordinate direction and its minimum in other spatial coordinate directions, and the bond lengths and bond angles of all atoms in the structure fall within a chemically reasonable range. The chemically reasonable range of bond lengths and bond angles is determined by the following criteria: bond lengths match the average value ± 3 standard deviations of the corresponding bonding type in authoritative structural databases (Cambridge Structure Database (CSD) and Protein Structure Database (PDB) in this invention). For bonds to be broken / formed in transition states, the range can be relaxed to the average value ± 5 standard deviations. At the same time, all bond lengths must meet the absolute threshold of the corresponding bond type. Bond angles match the characteristic range of the hybridization type of the central atom. For ring structures, the bond angle can be relaxed to the baseline of the corresponding hybridization type ± 10°.

[0102] After structural optimization, frequency analysis is performed on the optimized transition state structure. The vibrational frequency data of the optimized transition state structure are output, with negative frequencies defined as imaginary frequencies, and the corresponding atomic vibrational states are presented. The transition state structure is judged to be acceptable based on the following criteria: the existence of only one imaginary frequency, and the direction of the atomic vibration corresponding to this imaginary frequency aligns with the stretching direction of the bond to be broken and the contraction direction of the bond to be formed in the reaction. If these conditions are met, the structure is deemed acceptable; otherwise, it is readjusted. If the structure is deemed unacceptable, the interpolation ratio of reactants and products is readjusted to generate a new initial transition state structure.

[0103] For the newly generated initial transition state structure, repeat the steps of structure optimization, frequency analysis and qualification judgment until a qualified transition state is obtained.

[0104] Transition state generation methods for medium-complexity scenarios include: The system receives basic structural data of the intermediate and configures simulation parameters such as simulation temperature and time step for the corresponding molecular dynamics simulation process. Based on the molecular mechanical force field, it calculates the interactions between atoms in the intermediate's basic structure. It simulates the random trajectories of each atom at a set temperature and saves a structural snapshot at regular simulation intervals.

[0105] The core purpose of simulation parameter configuration is to ensure that the simulation can generate diverse and effective candidate structures by setting parameters reasonably. The parameter configuration standard is that the simulation temperature can make the atoms vibrate moderately without destroying the basic molecular framework, and the time step matches the atomic vibration period to ensure simulation accuracy.

[0106] From the saved structural snapshots, 20 snapshots with significant structural differences were selected as candidate transition state structures. These 20 candidate transition state structures were then imported in batches, and all functionals in the density functional approximation consensus set were called to calculate the energy of each candidate transition state structure. Based on the weight coefficients of each functional, the multiple energy calculations for the same candidate transition state structure were weighted and fused to obtain the comprehensive energy of the structure. The 20 candidate transition state structures were then sorted in ascending order of comprehensive energy.

[0107] The criterion for judging obvious structural differences is that the maximum value of the coordinate deviation of corresponding atoms in any two snapshots is greater than or equal to the set coordinate deviation threshold.

[0108] Examine the ligand binding and interatomic spacing of each candidate transition state structure to determine whether the structure is reasonable. Structures with abnormal coordination numbers or unreasonable interatomic spacing are eliminated, and the top 5 high-quality candidate structures with the highest overall energy are retained.

[0109] The criterion for judging whether a structure is reasonable is that the coordination number of the central metal atom is within the characteristic coordination number range of this type of intermediate, and the distance between any two non-bonding atoms is greater than or equal to the sum of their van der Waals radii. If this condition is met, the structure is retained; otherwise, it is discarded.

[0110] These five high-quality candidate structures are used as the initial structures for a new round of molecular simulations. The process of generating snapshots through simulation, calculating the overall energy, and selecting high-quality structures is repeated iteratively. The iteration stops when the set number of iterations is reached, or when the overall energy difference between the high-quality structures retained in two consecutive rounds is less than or equal to a set energy difference threshold. Otherwise, the iteration continues.

[0111] Among the final selected high-quality candidate structures, some or all may satisfy the condition of having only one imaginary frequency corresponding to the response coordinate. Frequency analysis is performed on each of these high-quality candidate structures to screen out all structures that satisfy the condition of having only one imaginary frequency corresponding to the response coordinate. The structures that meet the imaginary frequency condition are then sorted from low to high based on their overall energy. The structures with the lowest overall energy are then output as the optimal transition states.

[0112] Transition state generation methods for highly complex scenarios include: The transition state generation method for high-complexity scenarios adds a functional cross-validation mechanism to the transition state generation method for medium-complexity scenarios.

[0113] After iterating through the same method of generating transition states in medium-complexity scenarios—simulating snapshot generation, calculating comprehensive energy, and screening high-quality structures—high-quality candidate structures are obtained through each round of screening.

[0114] Invoke all functionals in the density functional approximation consensus set, allocate independent computing resources to each functional, and perform energy calculations on each high-quality candidate structure simultaneously to ensure that the calculations of multiple functionals are completed synchronously in order to control the total time consumption.

[0115] After the calculation is completed, energy data of the same high-quality candidate structure under different functionals are extracted. Variance calculation is then performed on the extracted energy data. First, the average value of the energy data set is calculated, then the sum of squares of the differences between each energy value and the average value is calculated, and finally, the variance is obtained by dividing by the number of functionals. The stability and consistency of the energy data are determined by a variance result ≤ a set variance threshold. If this threshold is met, the structure is considered reasonable; otherwise, the structure may be abnormal. Structures with variance exceeding the threshold are removed, and only structures with acceptable variance and optimal overall energy are retained as the initial structure for the next iteration.

[0116] The process of initial structure simulation, screening for high-quality structures, and functional cross-validation is repeated until the set number of iterations is reached. After the set number of iterations is reached, frequency analysis is performed on the remaining candidate structures to screen out structures that meet the following criteria: a single imaginary frequency and corresponding reaction coordinates (imagine frequency is the atomic vibration mode corresponding to the imaginary frequency, where the atom to be broken stretches along the bond axis and the atom to be formed contracts along the bond axis, and the consistency between the vibration direction vector and the preset reaction coordinate vector is greater than or equal to the set consistency threshold; structures that meet these criteria are retained, and those that do not are discarded). Several highly reliable transition states that have undergone multifunctional verification are output.

[0117] In any scenario, frequency analysis is an essential verification step. The number and direction of imaginary frequencies directly determine the rationality of the transition state. If two or more imaginary frequencies appear, it indicates that the structure has additional unstable sites (such as excessive ligand distortion), requiring re-optimization. If the direction of the imaginary frequencies does not correspond to the reaction coordinates (e.g., the target is C-Br bond breaking, but the imaginary frequencies correspond to ligand rotation), it indicates that a transition state of other reactions (such as ligand rearrangement) has been generated, and the sampling direction needs to be adjusted. Only transition states that pass the frequency analysis can proceed to the subsequent reaction pathway construction stage.

[0118] Step S6. First, generate a complete path based on the structural matching between the intermediate and the transition state in step S5. Use the atomic coordinate RMSD value as the matching criterion to form a continuous sequence from reactants to products and record energy change data. Then, compare the original paths pairwise. First, filter and eliminate redundant paths according to the energy barrier difference threshold of key steps. Then, retain different types of paths by judging the mechanism difference through the bond change mode sequence to obtain the core path. Finally, call the experimental database to match the data at a set frequency, classify the deviation level between the path and the experimental data, and carry out differentiated iterative optimization to form a closed loop to continuously improve the accuracy of the path and provide a reliable foundation for the subsequent construction of the catalytic reaction network.

[0119] First, a complete reaction pathway is generated, based on the structural matching of the intermediate and the transition state.

[0120] The reaction pathway is a continuous process from reactants to products, with each step involving a transition from one structure to a transition state and then to the next structure. The product of the preceding structure must be the same substance as the reactant of the following structure; that is, the structures must be completely matched, with identical atomic types, numbers, and bonding relationships.

[0121] The process iterates through all intermediates and transition states obtained in step S5, using an atomic coordinate RMSD value ≤ 0.05 nm as the criterion for structural matching. First, the precursor and product structures corresponding to each transition state are extracted, and the atomic coordinate RMSD values ​​between the product of the preceding structure and the reactant of the following structure are calculated. If the RMSD value meets the criterion, the two structures are considered structurally matched, and these two adjacent steps are joined together. This process is repeated for all structures, and successfully matched structures are sequentially linked in the order of preceding structure, transition state, and subsequent structure. The links are gradually connected starting from the reactant and extending to the product, forming a complete reaction pathway sequence (fragments that cannot be linked with other structures are directly discarded to ensure pathway continuity). Energy change data for each step of each reaction pathway are recorded, including the relative energy of the intermediate and the energy barrier of the transition state. The relative energy reflects the stability of the intermediate, and the energy barrier is the transition state energy minus the energy of the preceding intermediate; the higher the energy barrier, the more difficult the reaction to occur. This creates a complete pathway information in two dimensions: structure and energy, avoiding pathway breaks or structural inconsistencies.

[0122] Traverse all generated original paths and compare each path pairwise. By comparing the difference in reaction energy barriers between the two reaction paths, if the difference in reaction energy barriers between the two reaction paths is less than a set threshold, retain the one with the lower relative energy. Then compare the bond changes between the two paths. If there are differences in the bond breaking type, bond formation type, or bond position in any step, retain both.

[0123] In this invention, the threshold for the difference in reaction energy barriers is set at 3 kJ / mol, and reaction paths with a difference in reaction energy barriers ≥ 3 kJ / mol are retained. When comparing two paths, the critical steps of each path are first identified. The critical steps are the rate-determining steps with the highest energy barriers in the reaction (the slowest step in a chemical reaction that determines the overall reaction rate) or core transformation steps. Then, the energy barrier difference between the critical steps of the two paths is calculated. If the difference is < 3 kJ / mol, and the energy barrier values ​​and reaction mechanisms of all other steps in the two paths are the same, it indicates that the reaction difficulty and trend of the two paths are basically the same, and they are redundant paths. Only the one with lower energy and more likely to occur is retained. If the difference is ≥ 3 kJ / mol, it indicates that the reaction difficulty of the two paths is significantly different, and both should be retained to reflect the reaction probability under different conditions.

[0124] Then, the bond change information of all steps in each path is extracted to clarify the specific type and location of the bonds that need to be broken and the bonds formed in each step. Then, all the bond change information of each path is integrated in the order of steps to form a bond change pattern sequence specific to that path. Finally, the bond change pattern sequences of two paths are compared. If the bond breaking type, bond formation type or bond position of any step in the sequence is different, it is determined that the reaction mechanism of the two paths is different. Paths with such different mechanisms should be retained to cover different reaction path types such as main reaction path and side reaction path.

[0125] By analyzing the differences in reaction energy barriers and bond changes along reaction pathways, 20-30 core pathways are selected from dozens of original pathways, avoiding redundant calculations while ensuring no key mechanisms are overlooked.

[0126] After screening the reaction pathways, the reaction pathways are calibrated.

[0127] The calibration process consists of three steps.

[0128] The first step involves accessing experimental databases (such as the Reaxys database, NIST Chemistry WebBook, etc.). These databases contain tens of thousands of experimental measurement data for catalytic reactions, covering intermediate structures, transition state characteristic peaks, reaction barriers, etc. The corresponding experimental data are then quickly matched based on the reaction type and catalyst type.

[0129] The second step involves classifying the degree of deviation by calculating the discrepancy between the generated path and the experimental data, categorizing them according to two criteria. Small deviations are defined as an average deviation of ≤0.1 nm for structural atomic coordinates and ≤5 kJ / mol for energy barrier deviation, indicating a high degree of agreement between the calculated results and the experiment, and a reliable path. Medium deviations are defined as an average deviation of 0.1 to 0.2 nm for structural atomic coordinates or a deviation of 5 to 10 kJ / mol for energy barrier deviation, indicating some deviation and requiring parameter adjustment and optimization. Large deviations are defined as an average deviation of >0.2 nm for structural atomic coordinates or a deviation of >10 kJ / mol for energy barrier deviation, indicating a significant difference between the calculated results and the experiment, potentially indicating a problem with the model or method, and requiring in-depth optimization.

[0130] The third step involves iterative optimization based on deviation levels, employing differentiated treatment for different deviations to form a closed-loop iteration that continuously improves accuracy. When the deviation is small, the path is directly retained, while the cost function weights are updated, increasing the feature weights corresponding to this type of path. Subsequent generation of similar paths prioritizes referencing these features, optimizing the path generation logic. When the deviation is moderate, the process backtracks to the transition state generation stage to analyze the source of the deviation. If it's a structural deviation (i.e., the average deviation of the structural atomic coordinates exceeds the standard), the number of transition state sampling rounds is increased to explore more configurations. If it's an energy deviation (i.e., the energy barrier exceeds the standard), the functional in the density functional approximation consensus set is replaced to improve energy calculation accuracy, and the transition state is regenerated and the path reconstructed. When the deviation is large, the classification prediction model is retrained in reverse, using the intermediate features corresponding to the path and experimental data as new training samples to update the model parameters. Then, the intermediate classification, density functional approximation consensus set construction, and subsequent processes are re-executed to correct the deviation at its root.

[0131] Reaction path calibration is a crucial step in improving accuracy. By comparing with experimental data, calculation biases are identified and corrected. The calibration frequency is set to trigger once every 20 core paths generated. Calibration is performed after generating 20 paths, and then again every 20 newly added paths, ensuring timely correction of biases. This closed-loop iterative optimization continuously reduces calculation bias, making the generated paths increasingly closer to experimental realities, providing a reliable foundation for the final construction of the catalytic reaction network.

[0132] Step S7. Iteratively converge the reaction paths that have undergone multiple calibration processes in Step S6. Integrate the reaction paths with smaller deviations in multiple consecutive path calibrations that meet the rated threshold into a catalytic reaction network and label them. Calculate the accuracy of the DFT approximate dynamic selection of the catalytic reaction network and the median of the overall reaction barrier calculation error. If both the accuracy of the DFT approximate dynamic selection of the catalytic reaction network and the median of the overall reaction barrier calculation error are satisfied, output the catalytic reaction network and send a notification report. If neither of these thresholds is satisfied, return to the path calibration in Step S6 for optimization and adjustment.

[0133] First, perform iterative convergence judgment. In three consecutive path calibrations, if the proportion of paths with smaller deviations in a single calibration is greater than or equal to 80%, it indicates that the calculation result is reliable. If this proportion requirement is met for three consecutive calibrations, it indicates that the entire calculation system has reached a stable and reliable state, with no random deviations affecting it, and the iterative convergence has been achieved, and the final network integration can be initiated. If the requirement is not met, path calibration and iterative optimization need to be performed again until the target is met for three consecutive calibrations, in order to avoid the overall network deviation caused by local optimization.

[0134] The system integrates and optimizes the core reaction pathway, intermediate structural data, transition state data, and DFT calculation parameters (intermediate structural data includes atomic coordinates of the geometric structure, electron cloud distribution of the electronic structure, and relative energy of the energy data. Transition state data includes structural parameters, imaginary frequency values, and energy barrier data. DFT calculation parameters include functional type, weighting coefficients, and solvent dielectric constant). Intermediates are used as nodes, each labeled with its intermediate ID and relative energy, with colors used to distinguish energy levels. Reaction steps are used as edges, each labeled with its reaction type, energy barrier value, and bond change mode. Reactants are set as starting nodes, and products as ending nodes, forming a complete network link from reactants to products. This visually presents the reaction mechanism and energy change trends, constructing a catalytic reaction network.

[0135] After the network link is formed, the performance of the network link is evaluated to determine the accuracy of the DFT approximate dynamic selection and the median of the overall reaction barrier calculation error.

[0136] The accuracy of the DFT approximate dynamic selection is an indicator of how well the reaction pathways calculated by the dynamically selected DFT functional match the experimental data. First, the total number of core pathways to be evaluated is determined. Then, the number of pathways that meet the matching criteria of structural RMSD ≤ 0.1 nm and energy barrier deviation ≤ 5 kJ / mol is selected. The result of dividing the number of matching pathways by the total number of core pathways and multiplying by 100% is the accuracy of the DFT approximate dynamic selection. In this invention, an accuracy of ≥ 90% is required.

[0137] The median of the overall reaction energy barrier calculation error is obtained by statistically analyzing the energy barrier errors of all transition states, i.e., subtracting the experimental values ​​from the calculated functional values ​​and taking the median. Choosing the median avoids the influence of extreme values ​​and is more objective than the average, thus reflecting the overall accuracy of the energy calculation. In this invention, the median is required to be less than or equal to 5 kilojoules per mole. If the median is greater than 5 kilojoules per mole, the functional weights of the density functional approximation consensus set need to be adjusted or the functional needs to be replaced, and the transition state energy needs to be recalculated.

[0138] Both the accuracy of the DFT approximate dynamic selection and the median error of the overall reaction barrier calculation must meet the required performance indicators before proceeding to the next step. If either indicator fails to meet the requirements, the process will backtrack to the path calibration in step S6 for optimization and adjustment until both indicators meet the standards. This ensures that the final constructed catalytic reaction network has reliable accuracy and practicality.

[0139] The network link that simultaneously meets both the accuracy of the DFT approximate dynamic selection and the median error of the overall reaction barrier calculation is the catalytic reaction network diagram output by this invention. The output format of the catalytic reaction network diagram includes a visualization file, a data table, and a calculation report. The output catalytic reaction network diagram can be associated with all output files via task ID and sent to the user-specified storage location. Simultaneously, a task completion notification is generated, which can be sent via email, SMS, or system message to ensure timely receipt and use of the results. The notification content includes the task ID, total execution time, core performance indicators, file storage path, and contact person.

[0140] The visualization files support 3D structure display and path navigation. The 3D display allows rotation and zoom to view the molecular configuration of intermediates and transition states. The path navigation dynamically demonstrates the structural changes from reactants to products according to the reaction sequence. The format supports pdb or xyz formats that can be opened by commonly used molecular visualization software.

[0141] The data tables present all core parameters in Excel or CSV format, including intermediate ID, atom coordinates, relative energy, transition state ID, imaginary frequency, energy barrier, DFT functional type, and weights. The tables are sorted by path, intermediate, and transition state associations, facilitating subsequent data analysis or importing into other computational software.

[0142] The calculation report, in PDF format, details each step of the calculation process, including initial parameter settings, key thresholds for intermediate enumeration and selection, selection criteria for DFT functionals, bias analysis of path calibration, and performance evaluation results. The report includes data sources and calculation methods, facilitating users' ability to trace calculation details and reproduce results, and also complies with the standardized requirements for publishing scientific research findings.

[0143] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for constructing a transition state and catalytic network, characterized in that, Includes the following steps: Step S1. Obtain the basic data of the target catalytic system uploaded by the user, load the built-in database of elementary reaction rules containing multiple bond breaking and formation modes, parse the data and extract microstructure information, preprocess the data and unify the data format, establish a dedicated task ID and parameter index table, generate logs that record initial parameters, start time and system resource allocation and back them up to local and cloud simultaneously. Step S2. Use graph theory algorithm to traverse the bond interaction combinations between catalyst and substrate, combine elementary reaction rules to screen breakable bond pairs to enumerate reaction intermediates, calculate the relative energy and bond change cost of intermediates, and screen and remove duplicates of enumerated intermediates based on preset energy threshold, metal stable coordination number range, atomic spatial steric hindrance conditions and electron cloud overlap principle to obtain a set of screened intermediates. Step S3. Based on quantum chemical calculation and machine learning models, extract the local metal features of the intermediate and the global features of the substrate, quantize the feature data to the 0-1 range to obtain a standardized vector, remove abnormal samples, and obtain the feature dataset in which the samples and intermediates have a one-to-one mapping relationship. Step S4. Call the classification prediction model based on the random forest algorithm, divide the intermediates into 9 basic categories according to the relevant information of the number of transition metal nuclei and the tooth size of ligands, evaluate the error and efficiency of several DFT functionals in the energy calculation of various intermediates, select the optimal functional, calculate and label the weight coefficient of each functional according to the reciprocal of the error, and construct a density functional approximate consensus set. Step S5. Based on the DFT approximation error value and the number of intermediate atoms, intermediates are divided into three complexity scenarios: low, medium, and high. Differentiated transition state generation methods are adopted so that intermediates in the three different complexity scenarios obtain transition states using different methods. For the low complexity scenario, the initial transition state is generated by interpolation, the structure is optimized by calling the functional with the highest weight, and the imaginary frequency characteristics are verified by frequency analysis to determine the qualified transition state. For the medium complexity scenario, candidate structures are generated by molecular dynamics simulation, and high-quality structures are obtained through energy calculation and structure screening and iterative optimization. Finally, the low-energy structure that meets the imaginary frequency condition is selected as the optimal transition state. For the high complexity scenario, functional cross-validation is added on the basis of the medium complexity method. Abnormal structures are eliminated through multifunctional parallel computation and variance analysis, and the transition state is selected after iteration. Step S6. First, generate a complete path based on the structural matching between the intermediate and the transition state in step S5. Use the atomic coordinate RMSD value as the matching criterion to form a continuous sequence from reactants to products and record energy change data. Then, compare the original paths pairwise. First, filter and eliminate redundant paths according to the energy barrier difference threshold of key steps. Then, retain different types of paths by judging the mechanism difference through the bond change mode sequence to obtain the core path. Finally, call the experimental database to match the data at a set frequency, classify the deviation level between the path and the experimental data, and take differentiated iterative optimization to form a closed loop. Step S7. Iteratively converge the reaction paths that have undergone multiple calibration processes in Step S6. Integrate the reaction paths with smaller deviations in multiple consecutive path calibrations that meet the rated threshold into a catalytic reaction network and label them. Calculate the accuracy of the DFT approximate dynamic selection of the catalytic reaction network and the median of the overall reaction barrier calculation error. If both the accuracy of the DFT approximate dynamic selection of the catalytic reaction network and the median of the overall reaction barrier calculation error are satisfied, output the catalytic reaction network and send a notification report. If neither of these thresholds is satisfied, return to the path calibration in Step S6 for optimization and adjustment.

2. The method for constructing a transition state and catalytic network according to claim 1, characterized in that, In step S1, the core data received includes a CIF format crystal structure file uploaded by the user that records the space group and three-dimensional coordinates of the catalyst crystal, an XYZ format file that records the molecular structure, and a substrate SMILES format string that expresses the molecular topology. Material property data must include at least the substrate purity parameter to identify and remove impurities whose content is below a set threshold. The reaction condition parameters cover the temperature range, pressure threshold and solvent type, and are transformed into calculation boundary parameters through density functional theory. Data preprocessing includes using a built-in parser to extract atomic coordinates, atomic connections, and electronic states from the structure file, correcting minor deviations in atomic coordinates caused by measurement errors, and converting all data into a unified format that the system can recognize. Establishing the associated index involves building a parameter index table to associate rule IDs, intermediate pre-numbers, and computing node information, and generating logs that record initial parameters, task start times, and system resource allocations, which are then synchronously backed up to local and cloud storage.

3. The method for constructing a transition state and catalytic network according to claim 1, characterized in that, In step S2, when constructing the molecular graph model, atoms are labeled with atom-related information as nodes with attributes, and chemical bonds are labeled with chemical bond-related information as edges with parameters. When screening for breakable bond pairs, first screen the bonds that are directly connected between the metal-ligand of the catalyst and the functional group of the substrate from the molecular diagram, and then pair them according to the combination forms in the rules, retaining only the combinations whose bond energies are within the allowable range of the rules. The intermediate's unique ID format is reaction type-metal type-serial number, with the serial number increasing in the order of generation. The intermediate output is an XYZ format molecular structure file. The first line of the file indicates the total number of atoms, the second line indicates the unique ID, and each subsequent line records the atom symbol and the X, Y, and Z coordinates in sequence. The energy threshold is set according to the system type, and the energy cost of bond change is calculated by subtracting the total bond energy of newly formed chemical bonds from the total bond energy of breaking all chemical bonds.

4. The method for constructing a transition state and catalytic network according to claim 1, characterized in that, In step S2, when checking the stable coordination number of metals, each metal is judged according to its inherent stable coordination number range. If it exceeds the range, it is judged as a structural anomaly and is removed. The shortest distance between atoms is calculated as the distance between the ligand atom and the substrate atom. If this distance is less than the sum of the van der Waals radii of the two atoms, the steric hindrance is deemed too great and the atom is rejected. The bond formation rationality check is based on the principle of electron cloud overlap, and bond formation structures that violate this principle are directly eliminated; Molecular fingerprint similarity calculation is used for deduplication. When the molecular fingerprint similarity between two intermediates is greater than or equal to a set threshold, the individual with lower energy is retained and duplicate individuals are removed.

5. The method for constructing a transition state and catalytic network according to claim 1, characterized in that, In step S3, the local characteristics of the metal are obtained through density functional theory calculations, including at least the electron cloud distribution of the metal nucleus, the electronegativity contribution of the coordinating atoms, the bond order and bond energy of the metal-ligand bond, the number and distribution of empty orbitals of the metal atoms, the bond angle of the metal-ligand bond, and the charge distribution of the coordinating atoms. The global characteristics of the substrate are extracted using molecular topology analysis tools and include at least several indicators such as the topological index of the molecule, the type and number of functional groups, the uniformity of charge distribution, the molecular dipole moment, the bond length and bond angle distribution, the molar mass and volume, and the length of the conjugated system. If the proportion of missing features in a certain intermediate exceeds a preset threshold, the intermediate is identified as an abnormal sample and removed. For the few missing features in the retained samples, the average feature value of the same type of intermediate is used to fill them in, and finally a high-dimensional feature dataset is formed that is associated with the intermediate ID one by one.

6. The method for constructing a transition state and catalytic network according to claim 1, characterized in that, In step S4, the classification prediction model is constructed based on the random forest algorithm, and its training data comes from a large amount of intermediate feature data of the catalytic system; Metal core counts are classified into three categories: single-core, dual-core, and multi-core. Ligand dentation is classified into three types: monodentate, dipaldentate, and multidentate. Nine basic categories are formed by combining two dimensions: the number of metal nuclei and the dentation of ligands. For intermediates containing mixed dentation ligands, they are classified into subcategories of the corresponding basic categories based on their dominant structural features. For each type of intermediate, we evaluated a variety of functionals, including B3LYP, B3LYP-D3, PBE, PBE0, M06-2X, TPSS, B97X-D, CAM-B3LYP, B97D3, and M11. The evaluation metrics included the average error of energy calculations with reference to high-precision quantum chemical calculation methods and the average computation time under the same hardware configuration. The optimal functional is selected by comprehensively considering error and efficiency, ranking and scoring the functionals, and then choosing the one with the highest score. Constructing a density functional approximate consensus set involves selecting 3-5 optimal functionals for each type of intermediate and assigning a weight coefficient to each functional based on the inverse of its error. The calculation formula is: weight coefficient of a functional = (1 / error of the functional) / (1 / sum of errors of all selected functionals).

7. The method for constructing a transition state and catalytic network according to claim 1, characterized in that, In step S5, the criteria for scene division are as follows: The low-complexity scenario is one where the DFT approximation error is ≤5kJ / mol and the number of intermediate atoms is ≤50; Medium-complexity scenarios are characterized by DFT approximation errors between 5 and 15 kJ / mol or intermediate atom counts between 50 and 100. High-complexity scenarios are those where the DFT approximation error is >15 kJ / mol or the number of intermediate atoms is >100; For low-complexity scenarios, an initial transition state is generated by interpolating the coordinates of reactant and product structures. The density functional approximates the functional with the highest weight in the consensus set for structural optimization. Frequency analysis is used to verify whether there is a unique imaginary frequency with the same direction as the reaction coordinates to determine a qualified transition state. For medium-complexity scenarios, candidate structures are generated through molecular dynamics simulations, their combined energies are calculated, and high-quality structures are selected for iterative optimization. Finally, low-energy structures that meet the imaginary frequency condition are selected as the optimal transition state through frequency analysis. For high-complexity scenarios, a functional cross-validation mechanism is added to the medium-complexity method. That is, the energy of candidate structures is calculated in parallel using all functionals in the consensus set, and abnormal structures with unstable energy data are eliminated through variance analysis. After iteration and frequency analysis, a highly reliable transition state is output.

8. The method for constructing a transition state and catalytic network according to claim 1, characterized in that, In step S6, a certain set value of the root mean square deviation of atomic coordinates is used as a judgment condition to splice the product of the previous structure and the reactant of the next structure that match the structure of the intermediate and the transition state, thereby forming a continuous sequence from reactant to product, and recording the relative energy and energy barrier data of each step. The selection of core paths includes: First, comparing the original paths pairwise. If the difference in energy barriers of the key steps of the two paths is less than a set threshold and the energy barriers of other steps are the same as the mechanism, they are considered redundant paths, and only the one with the lower energy is retained. Second, comparing the bond change pattern sequence of the retained paths. If the bond breaking type, bond formation type or bond position of any step is different, they are determined to be paths with different mechanisms and are all retained. In path calibration, the differences between the generated path and the experimental data in terms of the average deviation of the structural atomic coordinates and the energy barrier deviation are calculated, and the deviations are divided into three levels: small, medium and large. Differentiated iterative optimization strategies are adopted for different deviation levels, such as updating feature weights, increasing transition state sampling, replacing DFT functionals or retraining classification prediction models. The calibration frequency is set to be triggered after a certain number of core paths are generated.

9. The method for constructing a transition state and catalytic network according to claim 8, characterized in that, The differentiated iterative optimization strategy is as follows: For paths with smaller deviations, update the cost function weights to prioritize their corresponding features; For paths with moderate deviation, if the deviation is structural, increase the number of transition state sampling rounds; if the deviation is energy, replace the functional in the density functional approximation consensus set. For paths with large deviations, the classification prediction model is retrained, and the relevant data is used as new samples to update the model parameters.

10. The method for constructing a transition state and catalytic network according to claim 1, characterized in that, In step S7, the criterion for determining iterative convergence is: in three consecutive path calibrations, the proportion of paths with smaller deviations is greater than or equal to 80%. Using intermediates as nodes and reaction steps as edges, nodes are labeled with ID and relative energy and color-coded to distinguish energy levels, while edges are labeled with reaction type, energy barrier value and bond change mode, forming a complete network link from reactants to products. The performance evaluation includes at least the accuracy of the DFT approximate dynamic selection and the median of the overall reaction barrier calculation error. If either indicator fails to meet the standard, the process is backtracked to the path calibration step in step S6 for optimization.