A molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification

Abstract

This invention belongs to the field of molecular dynamics simulation post-processing technology and discloses a molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification. Addressing the problems of high memory load in long trajectory processing, low automation of mechanism classification, and disconnected verification in existing technologies, this invention first reads trajectory data to construct a dynamic molecular connectivity graph, performing temporal smoothing and substrate intelligent masking filtering on bonding relationships; decoupling reaction events through atomic mapping bipartite graphs and eliminating non-physically correlated fluxes; automatically identifying species adsorption states using a surface chemistry sensing algorithm; then automatically clustering reaction mechanisms using neural networks and geometric descriptors; finally, selecting representative conformations and generating smooth transition state paths using an image-dependent potential algorithm and outputting the structure file. This invention achieves dual denoising from signal to event, significantly improving the analysis efficiency and automation of trajectory analysis for complex reaction systems.

Description

Technical Field

[0001] This invention relates to the field of molecular dynamics simulation post-processing technology, specifically to a molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification. Background Technology

[0002] Molecular dynamics (MD) simulations of atomic trajectories are crucial tools for studying chemical reaction mechanisms at the microscopic scale. With the widespread application of reaction force fields (such as ReaxFF) and machine learning potential functions, the scale of simulation systems has expanded to the tens of thousands or even hundreds of thousands of atoms, and simulation durations have extended from picoseconds to nanoseconds. Faced with such massive and high-dimensional atomic trajectory data, how to efficiently, accurately, and with chemical semantics extract reaction mechanisms has become a bottleneck restricting computational chemistry data mining.

[0003] In existing technologies, reaction identification based on MD simulation trajectories mainly relies on statistical methods based on geometric criteria: by setting fixed bond lengths or bond order cutoff values, the presence of a chemical reaction is determined by monitoring changes in the atomic adjacency matrix over time. Some advanced methods (such as ReacNetGenerator) introduce Hidden Markov Models (HMMs) to temporally smooth atomic connection states, eliminating instantaneous bond length oscillations caused by thermal vibrations; simultaneously, graph traversal algorithms are used to identify molecular species and generate reaction networks. However, the aforementioned existing technologies still have the following significant drawbacks in practical applications, and these drawbacks have not yet been systematically resolved: 1. Memory redundancy and IPC communication bottlenecks in massively parallel computing.

[0004] Existing parallel analysis solutions based on Python objects (such as ReacNetGenerator) suffer from extremely high CPU load when processing long trajectories with hundreds of thousands of frames due to frequent object serialization and deserialization required for inter-process communication (IPC). Furthermore, repeated copying of the entire trajectory object by child processes often leads to memory explosions. Current technologies lack a low-level architecture that can achieve high-performance multi-process scheduling while ensuring memory safety through zero-copy sharing (child processes directly read the trajectory image from physical memory without generating additional copies) and binary stream compression (child processes compress the calculated molecular topology data into a binary stream before sending it back to the main process, greatly reducing the bandwidth load of IPC).

[0005] 2. Misjudgment of transient reversible events and data redundancy.

[0006] While existing methods can filter out random fluctuations in bond lengths near the cutoff value using Hidden Markov Models (HMMs), they cannot handle "dynamic recrossing" phenomena—where molecules rapidly cross energy barriers to form products, only to quickly revert to reactants due to metastable states. From a signal processing perspective, this represents two valid state transitions; however, from a reaction network construction perspective, such high-frequency, non-unidirectional reversible cycles introduce a large number of redundant events, significantly increasing the complexity of the reaction network and the analytical burden. Current techniques lack automated cleaning mechanisms for this type of "event-level" logic.

[0007] 3. Substrate interference and adsorption state identification are difficult in heterogeneous catalytic systems.

[0008] In solid-liquid or solid-gas heterogeneous catalysis simulations (such as Fischer-Tropsch synthesis on metal surfaces), the changes in interatomic geometric distances during the adsorption and desorption of reactant molecules on the catalyst surface are remarkably similar to the formation of chemical bonds. Existing methods (including ReacNetGenerator) fail to distinguish between the "environmental substrate" and "active molecules," often misinterpreting physical adsorption as a chemical reaction. Furthermore, due to the lack of dynamic perception of the adsorption state of species, they cannot automatically distinguish between classic heterogeneous catalysis mechanisms such as LH (where all participating species are in the adsorbed state) or ER (where gas-phase species react directly with adsorbed species), resulting in analytical results lacking deep chemical mechanism semantics.

[0009] 4. Low accuracy in identifying concurrent reactions.

[0010] Existing techniques typically aggregate events based on time frames, assuming that all topological changes occurring within the same time step belong to the same chemical reaction. In large-scale systems, multiple spatially distant and logically unrelated independent reactions (e.g., a decomposition reaction on the left side of the system and a polymerization reaction on the right side) often occur within the same frame. Existing processing methods incorrectly merge these unrelated local changes into a single "global event" containing dozens of molecules, leading to distorted reaction equations and an inability to reconstruct the true elementary reaction steps.

[0011] 5. The level of automation in classifying reaction mechanisms is low.

[0012] Existing analytical tools primarily output statistical data on atomic connections (such as the number of C-H bond breaks), lacking the ability to identify the chemical properties of reactions. Even after generating reaction networks, researchers still need to manually examine numerous structural snapshots to confirm reaction types (such as substitution, elimination, and addition), and cannot automatically cluster events with the same mechanism. This severely limits the efficiency of high-throughput mechanism screening.

[0013] 6. The subsequent quantum chemical verification is difficult to connect.

[0014] MD simulation trajectories are often used to discover potential reaction pathways, but confirming their validity often requires combining them with density functional theory (DFT) to search for transition state energy barriers (such as the Nudged Elastic Band (NEB) method). Currently, the process of extracting initial and final state structures from MD simulation trajectories, performing one-to-one atom matching, and interpolating to generate NEB input files still mainly relies on manual operation or scattered scripts, lacking integrated and automated means, which limits the high-throughput flow from classical mechanics simulation to quantum mechanics verification.

[0015] In summary, existing technologies suffer from several significant drawbacks when processing MD simulation trajectories of complex reaction systems (especially heterogeneous catalysis), including limited architectural performance, high noise levels, susceptibility to misinterpretation of concurrent reactions, lack of mechanistic semantics, and disconnected verification processes. Therefore, there is an urgent need for an analytical method that enables high-performance multi-process analysis, dual purification of signals and events, intelligent decoupling of concurrent reactions, automatic identification of reaction mechanisms, and seamless integration with subsequent quantum chemical calculations. Summary of the Invention

[0016] In view of this, the present invention provides a molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification, in order to solve the prominent problems of existing technologies in processing MD simulation trajectories of complex reaction systems, such as low processing performance, large noise interference, easy misjudgment of concurrent reactions, lack of mechanistic semantics, and disconnection of verification.

[0017] In a first aspect, the present invention provides a molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification, the method comprising: S1: Based on a zero-copy memory sharing and binary stream compression architecture, read molecular dynamics trajectory data from all time frames, extract the three-dimensional coordinates of all atoms from each time frame, calculate the interatomic distance, and determine the initial bonding relationship according to the preset bond length cutoff value; use a hidden Markov model for temporal smoothing, and introduce a substrate smart mask to filter the environmental background connection, generating a dynamic molecular connection graph. S2: Identify stable plateau periods in all time frames, compare the dynamic molecular connectivity graphs of adjacent time frames within the stable plateau period, and identify reactants and new products after chemical reactions; construct an atomic mapping bipartite graph and calculate connected components to isolate independent reaction events, and based on the element conservation filtering algorithm of chemical formula analysis, remove non-physical correlation fluxes between reactants and products that have no common elements to obtain purified chemical reaction events. S3: Arrange the chemical reaction events in chronological order to obtain an event list. Traverse the event list and, by identifying identity reactions, instantaneous reversible noise, and atomic index overlap, remove redundant chemical reaction events with stability below a preset stability threshold to obtain the final list of chemical reaction events. S4: The list of chemical reaction events is converted into a text sequence, semantic feature vectors are extracted using a neural network, and a reaction semantic map is constructed using a dimensionality reduction algorithm; at the same time, local structural features of the reaction center atoms in each event are extracted using geometric descriptors to construct a structural evolution landscape map; in the low-dimensional landscape map space, events with the same reaction mechanism are automatically clustered. S5: For each type of chemical reaction event with the same reaction mechanism, select a representative conformation, remove the influence of periodic boundary conditions, and generate an adaptive simulation box with a vacuum layer; use a nonlinear interpolation algorithm to generate a nonlinear interpolation path between the initial and final states of the species participating in the chemical reaction event, and output it as a structure file.

[0018] In one optional implementation, the timing smoothing process in step S1 includes: A hidden Markov model is constructed. For atom pairs that meet the preset bond length truncation condition within at least one time frame, the connection state sequence of the atom pairs in different time frames is taken as the observation sequence. The optimal hidden state sequence is calculated using the Viterbi algorithm. When the atom pair is in a non-bonding state, if an isolated bonding signal with a continuous frame number lower than a preset threshold appears in the observation sequence, it is identified as physical vibration noise and removed by the Viterbi optimal path.

[0019] In one optional implementation, the substrate smart mask filtering of the background connection in step S1 includes: During the topology analysis initialization phase, a one-dimensional Boolean mask vector is calculated and generated based on the preset catalyst element type or atom index number. During the atom adjacency search process in each frame, bitwise operation filtering is performed using the one-dimensional Boolean mask vector to automatically skip bonding determination requests involving the interior of the substrate atoms and their interfaces.

[0020] In one optional implementation, the element conservation filtering algorithm in step S2 includes: Analyze the Hill representation chemical formulas of reactant and product entities and extract the element sets; determine if the intersection of the reactant and product sets is empty; if the intersection is empty, forcibly disconnect the associated edges corresponding to the intersection of the reactant and product sets in the atom mapping bipartite graph.

[0021] In an optional implementation, step S3 further includes: Using a preset adsorption cutoff radius scaling factor, the nearest mirror image distance between reactant or product molecules and substrate atoms is calculated in each time frame. If the nearest mirror image distance is less than the adsorption cutoff radius scaling factor, the corresponding species is marked as an adsorbed species; otherwise, it is marked as a gaseous species. By comparing the changes in species states before and after the chemical reaction, the species are automatically identified and classified.

[0022] In one optional implementation: In step S5, the selected representative conformations are preprocessed, the preprocessing including: Cross-boundary stitching: The inverse unit cell matrix transformation method is used to detect and repair broken molecules that cross periodic boundaries, and representative conformations are reduced to geometrically continuous molecular entities. Adaptive cell construction: Calculate the minimum bounding box of the reaction center atomic clusters of representative conformations, and construct independent adaptive orthogonal simulation cells by extending a vacuum layer of at least 10 Å in each of the three-dimensional directions based on the minimum bounding box. Centroid Alignment: Calculate the geometric centroids of the initial and final states of the species, and then use a translation transformation to make the geometric centroids of the initial and final states of the species coincide with the geometric center of the adaptive orthogonal simulation cell.

[0023] In an optional implementation, the text sequence extraction in step S4 includes concatenating the SMILES strings of reactants and new products into text entries; the neural network is a Transformer model based on a bidirectional self-attention mechanism, which extracts reaction semantic feature vectors through self-supervised learning; the geometric descriptor is a SOAP descriptor, used to capture the local three-dimensional configuration evolution within a 5 angstrom radius of the chemical reaction center; and the UMAP algorithm is used to map the high-dimensional semantic feature vectors to a low-dimensional landscape map space.

[0024] In an optional implementation, the output of the structure file in step S5 includes: The output contains a multi-frame trajectory file in XYZ format, which includes key time frames before, during, and after the chemical reaction. Each frame contains the original frame number and reaction ID information from the original molecular dynamics trajectory. The output contains a continuous sequence of configurations generated by image-dependent potential interpolation, which is used by quantum chemistry software to search for initial path guesses of elementary reaction transition states.

[0025] Compared with the prior art, the present invention has the following beneficial effects: This invention achieves dual denoising at both the signal and event levels through Hidden Markov Models (HMMs) and iterative time window algorithms. Combined with physical logic filtering based on element conservation, it effectively addresses core pain points in existing technologies such as thermal vibration noise interference, kinetic redundancy, and non-physical spurious correlations. Simultaneously, this invention utilizes a high-performance zero-copy architecture and binary compressed streaming technology to significantly reduce memory overhead and communication latency during ultra-long trajectory analysis. Through substrate intelligent masking and surface sensing mechanisms, it achieves automatic identification and spatial decoupling of multiphase catalytic interfaces (LH / ER). Furthermore, it uses atomic mapping bipartite graphs to decompose concurrent spatial reactions at the same time point into independent elementary reactions. In addition, this invention combines neural network semantic and geometric feature-based dual-perspective clustering to achieve automatic classification of reaction types. Through periodic boundary condition repair, centroid alignment, and automated IDPP interpolation processes, it establishes a seamless connection from MD kinetic trajectories to quantum chemical verification, significantly improving the processing efficiency, accuracy, and full automation level of complex reaction system analysis. Attached Figure Description

[0026] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0027] Figure 1 This is a flowchart illustrating a molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification according to an embodiment of the present invention. Figure 2 These are evolutionary landscape maps and automatic event clustering distribution maps according to embodiments of the present invention; Figure 3 This is a clustering effect diagram of a reaction semantic map based on a neural network according to an embodiment of the present invention; Figure 4 These are population evolution curves of chemical bonds and chemical species over time, according to embodiments of the present invention. Figure 5 This is a statistical distribution diagram of reaction types after noise reduction according to an embodiment of the present invention; Figure 6 This is a Sankey diagram of the net flux of the reaction pathway according to an embodiment of the present invention; Figure 7 This is a complex reaction network diagram according to an embodiment of the present invention; Figure 8 This is an interactive network graph triggered based on a specific node according to an embodiment of the present invention. Detailed Implementation

[0028] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0029] This invention provides a systematic solution to the problems of noise redundancy, high difficulty in multiphase interface identification, and cumbersome mechanism classification in large-scale MD trajectory analysis, aiming to improve the accuracy of extracting chemical reaction mechanisms from trajectories and the efficiency of data processing.

[0030] To illustrate the technical solution of the present invention in more detail, the following embodiments use the molecular dynamics simulation trajectory of the Fischer-Tropsch synthesis reaction on the surface of an iron (Fe) catalyst as an example. The experimental system runs on Python 3.8 or later, and the hardware environment is a computing workstation equipped with a multi-core CPU.

[0031] Figure 1 This is a flowchart of a molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification according to an embodiment of the present invention. The flowchart includes the following steps: S1. Constructing a dynamic molecular connectivity graph: Based on a zero-copy memory sharing and binary stream compression architecture, molecular dynamics trajectory data at all time frames are read, the three-dimensional coordinates of all atoms are extracted from each time frame, the distance between atoms is calculated, and the initial bonding relationship is determined according to the preset bond length truncation value; a hidden Markov model is used for temporal smoothing, and a substrate intelligent mask is introduced to filter the environmental background connectivity to generate a dynamic molecular connectivity graph.

[0032] Optionally, the temporal smoothing process includes: constructing a hidden Markov model; for atom pairs that satisfy a preset bond length truncation condition within at least one time frame, taking the connection state sequence of the atom pairs in different time frames as the observation sequence; using the Viterbi algorithm to calculate the optimal hidden state sequence; when the atom pair is in a non-bonding state, if an isolated bonding signal with a continuous frame number lower than a preset threshold appears in the observation sequence, it is determined to be physical vibration noise and removed through the Viterbi optimal path.

[0033] Optionally, the substrate intelligent mask filtering environment background connection includes: in the topology analysis initialization stage, calculating and generating a one-dimensional Boolean mask vector according to the preset catalyst element type or atom index number; in the atom adjacency search process of each frame, performing bit operation filtering using the one-dimensional Boolean mask vector to automatically skip bonding determination requests involving the interior of substrate atoms and their interfaces.

[0034] In this embodiment, the molecular dynamics simulation trajectory of the Fischer-Tropsch synthesis reaction on the surface of an iron (Fe) catalyst is used as an example. The simulation system includes a plate model composed of thousands of iron atoms and hundreds of carbon monoxide (CO) and hydrogen (H2) gas molecules.

[0035] 1. Read molecular dynamics trajectory data in standard formats (such as XYZ or LAMMPS Dump) across all time frames. For each frame, the experimental system constructs an all-atom connectivity matrix using the natural truncation radius algorithm. For any pair of atoms in the all-atom connectivity matrix, compare their interatomic distance with a preset bond length truncation value. If the interatomic distance is ≤ the bond length truncation value, mark it as 1; otherwise, mark it as 0. Thus, each frame yields an initial atomic adjacency matrix representing the bonding relationships between atoms.

[0036] The steps for constructing an all-atomic connectivity matrix include: Cutoff radius pre-calculation: To avoid redundant calculations, the experimental system uses the natural cutoff radius algorithm to calculate the covalent radius of each element in the system during the initialization phase. By introducing a global bond length scaling factor, a cutoff value matrix based on element-pair definitions is generated.

[0037] Efficient Neighbor List Retrieval: To avoid the traditional O(N) error... 2 To address the performance bottleneck caused by brute-force calculation of inter-atomic distances, the experimental system employs a neighbor list algorithm based on a cell list. This algorithm only performs distance verification on spatially adjacent atomic pairs.

[0038] Bonding determination and masking: For each retrieved atom pair (u,v), calculate its Euclidean distance d. u,v If d u,v If the bond length truncation value of the corresponding element is less than or equal to the initial connection, then an initial connection is determined to exist. The key optimization point is: during the determination process, a bitwise operation with O(1) time complexity is performed synchronously for filtering. If the reactant u or product v is marked as a base atom, then the connection request is directly executed.

[0039] Sparse topological representation: The adjacency relationships generated in each frame are stored as a sparse edge list, recording only the atomic index pairs that form bonds. This sparse representation significantly reduces memory usage compared to a complete adjacency matrix.

[0040] 2. To eliminate transient bonding oscillations caused by thermal vibrations, a temporal smoothing process is performed on the temporal connection states of each pair of atoms based on the Hidden Markov Model (HMM). This algorithm forces state transitions to have a certain degree of persistence by setting extremely high state persistence probabilities (e.g., 0.99) and noise tolerance (e.g., 0.1).

[0041] Specifically, the time series smoothing process based on the Hidden Markov Model is as follows: Observation sequence extraction: For specific atom pairs in the trajectory, extract their bonding states within T frames (0 for disconnection, 1 for connection), and construct the observation sequence O = {o1,o2,…,o...} T}

[0042] Logarithmic probability matrix definition: To prevent numerical underflow at the computer's underlying level, the experimental system defines a state transition probability matrix (describing the persistence of bonding states) and an emission probability matrix (describing the probability of the observed signal deviating from the true physical state) in logarithmic space.

[0043] Dynamic programming initialization: At the first frame (t=0), the path probability table dp is initialized based on the initial distribution and the observation value o0 of the first frame.

[0044] Recursive optimization (recursive process): Utilizing Numba hardware-accelerated recursive logic, for each subsequent frame t, calculate the maximum probability of transitioning from each state in the previous frame to the current possible state (bonding or breaking), and simultaneously record the preceding state path that generates this maximum probability. The formula for calculating the maximum probability is expressed as: In the formula, Let represent the maximum probability when the experimental system is in state s at frame t; This indicates the possible state of the experimental system at frame t-1; Indicates from state The probability of transitioning to state s; This indicates the observed value in state s. The probability of; This represents the observation value in frame t, used to assess the probability of observing this value under different conditions by combining the emission probability.

[0045] Backtracking to find the optimal path: After processing the last frame, start from the termination state with the highest probability, backtrack according to the recorded path information, and obtain the globally optimal hidden state sequence (i.e. the real bonding / breaking sequence).

[0046] Topology Reconstruction Update: The connection relationship of the atom pair in each frame is redefined based on the hidden sequence generated by the Viterbi algorithm, thereby completely filtering out spurious bonding signals caused by atomic thermal motion crossing the truncation threshold.

[0047] 3. Spatial-level decoupling: Basis filtering and surface sensing based on O(1) Boolean mask.

[0048] To avoid misinterpreting physical adsorption as a chemical reaction, the experimental system needs to efficiently distinguish between catalyst substrate atoms and reactant molecules. In this embodiment, for a slab model containing thousands of Fe atoms, a Boolean masking mechanism with O(1) time complexity is introduced to optimize computational performance in large-scale systems.

[0049] Basis attribute pre-labeling and O(1) mask construction: During the initialization phase, the experimental system pre-constructs a one-dimensional Boolean mask vector in memory, equal in length to the total number of atoms in the system, based on the catalyst element type (Fe) or atom index number. The one-dimensional Boolean mask vector maps the basis attribute of each atom to a Boolean value (True is marked as a basis and False as a reactant molecule during the experiment). In each subsequent frame of topology reconstruction loop, the algorithm does not need to traverse the basis list of thousands of records, but instead performs subscript addressing by the unique index number of the atom, thereby instantly determining the atom attribute in constant time. This fast determination mechanism at the bit operation level enables the experimental system to automatically skip all invalid connection determinations involving the interior of the basis in large-scale high-frequency loops.

[0050] Connected component identification and species stripping: Perform a breadth-first search on the initial atomic adjacency matrix after filtering the substrate to find all molecular entities.

[0051] Gas-phase species identification: Connected components that do not contain any base atoms are identified as gas-phase reactants (such as free CO, H2 and intermediates such as CH3 and OH).

[0052] Adsorbed species identification: Non-substrate atom connected components that have geometric connections with at least one substrate atom (based on a preset adsorption_scale scaling factor) are defined as adsorbed molecules.

[0053] Surface chemical mechanism prediction: The experimental system monitors the evolution of the adsorption state of species before and after the reaction in real time. LH (Langmuir-Hinshelwood) mechanism: If all initial species involved in the reaction are in an adsorbed state, they are automatically marked as LH mechanism pre-options; ER (Eley-Rideal) mechanism: If it involves gaseous species directly attacking surface adsorbates, it is automatically classified as an ER mechanism pre-option.

[0054] After the above decoupling process, the resulting dynamic molecular connectivity map achieves precise focusing on catalytic reaction events, completely eliminating noise inside the substrate and physical adsorption interference, providing highly pure input data for subsequent "event purification".

[0055] S2. Decoupling of reaction events based on atomic mapping bipartite graph: Identify stable plateau periods in all time frames, compare the dynamic molecular connectivity graphs of adjacent time frames within the stable plateau periods, and simultaneously identify reactants and new products after chemical reactions; construct atomic mapping bipartite graphs and calculate connected components to isolate independent reaction events, and based on the element conservation filtering algorithm of chemical formula analysis, eliminate non-physical correlation fluxes between reactants and products that have no common elements to obtain purified chemical reaction events.

[0056] Optionally, the element conservation filtering algorithm includes: parsing the Hill expression chemical formulas of reactant and product entities and extracting the element set; determining whether the intersection of the reactant set and the product set is empty; if the intersection is empty, forcibly disconnecting the associated edges corresponding to the intersection of the reactant set and the product set in the atom mapping bipartite graph.

[0057] In this embodiment, the processing logic is as follows: Preliminary separation of independent reaction events (atomic mapping bipartite graph decoupling): By comparing the dynamic molecular connectivity graphs of adjacent stable plateau phases, the set of reactant entities U that disappeared in the previous state and the set of new product entities V that were generated in the subsequent state are identified. An atomic mapping bipartite graph is constructed with the set of reactant entities U as the left node and the set of new product entities V as the right node. If reactant u i ∈U and product v j If there is atomic index sharing between ∈V, then an associated edge e is established. ij .

[0058] Connectivity component decoupling: Calculate the connected components of the atomic mapping bipartite graph. Each independent connected component represents a spatially and logically isolated chemical reaction event. This step successfully isolates multiple concurrent, unrelated elementary reactions within the system.

[0059] Deep purification of atomic flow paths (element conservation filtering): For each chemical reaction event identified above, an element conservation filtering algorithm based on chemical formula analysis is further introduced to refine the determination of atomic flow logic within the event: regular expressions are used to parse the Hill expression chemical formulas of all species within the event (such as C). x H y O z Extract the element set contained in each species (e.g., the element set corresponding to species H is {H}, and the element set corresponding to species CO2 is {C,O}). For each atom mapping association edge e within a chemical reaction event...ij (from reactant u) i Pointing to product v j ), calculate the intersection of the two sets of species elements I=Set(u i )∩Set(v j If the intersection I is empty, the mapping path is determined to be invalid at the physical level (i.e., the atomic components in the reactants can never be converted into the product), and the system will forcibly remove the false associated traffic.

[0060] In Fischer-Tropsch synthesis simulations, this dual purification logic ensures the rigor of the analytical results. Example scenario: In a reaction involving H + HOCO -> H2 + CO2, although H and CO2 are within the same reactive component, since the set of H elements {H} and the set of CO2 elements {C,O} have no intersection, this invention automatically eliminates interfering connections from H to CO2, retaining only the paths H to H2 and HOCO -> CO2 that conform to the atomic conservation logic.

[0061] Spatial decoupling example: In large-scale Fischer-Tropsch synthesis simulation systems, multiple reactions may occur concurrently at different active sites on the catalyst surface within the same time window. For example, within a certain time window, H2 dissociation and adsorption occur on the left side of the system, while CO bond breaking and activation simultaneously occur on the far right side. Traditional methods incorrectly merge these spatially isolated atomic changes into a chaotic reaction containing both H2 and CO. This invention, by calculating the connected components of the atom mapping bipartite graph, can accurately decouple the H2→2H reaction on the left side. ∗ With CO→C on the right ∗ +O ∗ It can be broken down into two parallel elementary reaction events.

[0062] S3. Event-level logic purification based on iterative time window: The chemical reaction events are arranged in chronological order to obtain an event list. The event list is traversed, and redundant chemical reaction events with stability lower than a preset stability threshold are eliminated by identifying identity reactions, instantaneous reversible noise, and atomic index overlap, so as to obtain the final list of chemical reaction events.

[0063] Optionally, it also includes: using a preset adsorption cutoff radius scaling factor to calculate the nearest mirror image distance between reactant or product molecules and substrate atoms in each time frame; if the nearest mirror image distance is less than the adsorption cutoff radius scaling factor, the corresponding species is marked as an adsorbed species, otherwise it is marked as a gaseous species; by comparing the changes in species states before and after the chemical reaction, it automatically identifies and classifies them.

[0064] In this embodiment, the processing logic is as follows: First, an adsorption cutoff radius scaling factor is preset based on the catalyst surface characteristics. For each molecular species in each time frame, the minimum mirror image distance between all its atoms and the nearest substrate atom is calculated. If this distance is less than the adsorption cutoff radius scaling factor multiplied by the van der Waals radius of the major atoms in the molecule, the molecule is marked as an adsorbed species; otherwise, it is marked as a gaseous species.

[0065] Subsequently, by comparing the state changes of species before and after the chemical reaction, the system automatically classifies the reaction events: the transformation of gaseous molecules into adsorbed molecules is an adsorption reaction; the transformation of adsorbed molecules into gaseous molecules is a desorption reaction; and the interconversion between adsorbed molecules is a surface reaction. Based on this, the system outputs a list of reaction events with state labels, providing richer semantic information for catalytic reaction network analysis.

[0066] The experimental system does not record every minute geometric fluctuation. Instead, it uses a sliding window defined by the preset `min_stability` parameter to search for stable plateaus where species composition remains constant across all time frames. A transition event is recorded only when the duration of a change in species topology exceeds a preset threshold (5 frames in this example, corresponding to 50 femtoseconds). The initially identified events are arranged into an event list in chronological order, and a recursive cleanup is performed using the preset `ReactionCleaner` module: identifying and deleting spurious events where the reactant and product SMILES strings and chemical formulas are completely identical (such as conformational torsions that do not change the topological connections within the molecule); for two adjacent events (Event 1: A→B; Event 2: B→A), the system performs an atom index consistency check. If the atom numbers involved in the two events completely overlap and the time interval between their occurrences is less than the preset stability threshold, they are determined to be non-reactive kinetic re-entry noise.

[0067] Multi-level iterative convergence algorithm: After deleting a pair of inverse events, events that were originally not adjacent (such as events i−1 and i+2) will become new adjacent due to "collapse," which may generate secondary inverse pairs. The experimental system adopts a multi-round iterative mechanism, repeatedly executing the detection and deletion logic until the event linked list reaches a fully converged state at the topological level.

[0068] Data Analysis in this Example: In the Fischer-Tropsch synthesis trajectory analysis (reaction temperature 500-800 K) of this example, a total of 680 initial state transition records were identified. After recursive processing by the preset ReactionCleaner module, 342 records (50.3%) were identified as kinetic re-crossing noise and automatically removed, ultimately retaining 338 physically meaningful and valid reaction events. Threshold Scientificity: The preset stability threshold of 50 femtoseconds is determined based on the thermal fluctuation scale of the Fischer-Tropsch system. On the Fe catalyst surface, the effective diffusion and bonding of adsorbed species are typically on the picosecond scale, while the thermal fluctuation rebound (re-crossing) based on the energy barrier top is typically within 100 femtoseconds. Using 50 femtoseconds as the cutoff standard can accurately remove high-frequency physical noise while ensuring that no real concerted reaction steps are missed. Figure 4 As shown, the purified species evolution curve exhibits a clear kinetic trend: the number of reactant H2 molecules (blue curve) decreases smoothly, while the number of surface-active H atoms (red curve) increases accordingly, eliminating the violent sawtooth fluctuations originally caused by the rapid reversible dissociation-recombination equilibrium on the surface. This indicates that event-level logic purification successfully filters out non-reactive noise caused by kinetic re-transitions, making the reaction network more realistic and concise.

[0069] S4. Dual-view mechanism identification and automatic clustering: The list of chemical reaction events is converted into a text sequence, semantic feature vectors are extracted using a neural network, and a reaction semantic map is constructed using a dimensionality reduction algorithm; at the same time, local structural features of the reaction center atoms in each event are extracted using geometric descriptors to construct a structural evolution landscape map; in the low-dimensional landscape map space, events with the same reaction mechanism are automatically clustered.

[0070] Optionally, the text sequence extraction includes concatenating the SMILES strings of reactants and new products into text entries; the neural network is a Transformer model based on a bidirectional self-attention mechanism, which extracts reaction semantic feature vectors through self-supervised learning; the geometric descriptor is a SOAP descriptor, used to capture the local three-dimensional configuration evolution within a 5 angstrom radius of the chemical reaction center; and the UMAP algorithm is used to map the high-dimensional semantic feature vectors to a low-dimensional landscape graph space.

[0071] In this embodiment, the 338 purified valid response events are converted into SMILES text sequences, input into a pre-trained BERT architecture to extract semantic features, and then visualized using the UMAP algorithm for dimensionality reduction. For example... Figure 2 The diagram shows the evolution of chemical reaction events in structural space and the automatic clustering distribution diagram, as shown below. Figure 3 As shown, discrete reaction events are clustered into several clear mechanistic clusters, combined with... Figure 5The denoised reaction type statistical distribution diagram shows that this embodiment successfully resolved the distribution of key elementary steps in the Fischer-Tropsch synthesis reaction network: Adsorption and activation process: The experimental system identified high-frequency dissociative adsorption (89 times) and molecular adsorption (50 times), mainly corresponding to the surface dissociation of H2 and the initial adsorption of CO, which provided active intermediates for subsequent reactions.

[0072] Heterogeneous catalysis mechanism determination: Based on a surface sensing algorithm, the experimental system accurately distinguished and identified Surface L-H Hydrogenation (LH hydrogenation mechanism, 29 times) and Surface LH Transfer (LH transfer mechanism, 14 times). This signifies the co-evolution among surface-adsorbed species (such as CH4). x The hydrogenation process of H+ is the dominant mechanism in this system.

[0073] Surface dissociation and coupling: Surface dissociation (23 times) and surface L-H coupling (7 times) were captured, which intuitively reflect the activation and breaking of C−O bonds and the initial carbon chain growth behavior.

[0074] Micro-competitive paths: The system can even capture extremely low-frequency Surface ER Transfer (ER transfer mechanism, 2nd time) and Gas-phase Homolysis (gas-phase homogeneous splitting, 2nd time), demonstrating its ability to deeply explore secondary paths in complex systems.

[0075] As shown in Figure 6, the system generates a net flux diagram of the Fischer-Tropsch synthesis reaction pathway based on the element conservation filtering logic. Figure 6 It intuitively reproduces the actual atomic flow logic within the system: clearly demonstrating the flow from the initial reactant (CO) ∗ H2 ∗ ) to key intermediate (C ∗ H ∗ CH ∗ ), and eventually diverted to hydrocarbon products (CH4) ∗ (etc.) and oxygen-containing byproducts (H2O) ∗ The complete chain of the reaction is shown. The line width represents the net flux after deducting the reverse reaction, visually revealing the main pathway of the reaction. High-precision identification of isomers: The figure shows the identification via CHO. -1 with CHO -2 nodes distinguish isomers with the same chemical formula but different topological structures, demonstrating the system's ability to perform refined analysis of complex species. Physical authenticity verification: All flow paths strictly adhere to the law of conservation of atoms, ensuring the physical rigor of the extracted reaction mechanisms in complex multi-component systems.

[0076] In addition, this experimental system also includes a visualization module that integrates a node-triggered interactive full-path tracing strategy. The HTML-format reaction map generated by the experimental system incorporates a bidirectional topological traversal algorithm. Users can click on any species node in the map with the mouse, and the experimental system will capture the interactive command in real time and execute cascading highlighting.

[0077] Specifically, when a specific intermediate (such as CH*) is selected, the experimental system will automatically trace back all its upstream source paths and simultaneously track all its downstream destinations and evolutions. At the same time, it will perform grayscale masking on irrelevant interference paths in the graph, such as... Figure 7 , 8 As shown, this interactive traceability mechanism allows researchers to precisely locate and analyze the complete life cycle of specific chemical components in a "what you see is what you get" manner within large-scale and complex reaction networks, significantly improving the intuitiveness and efficiency of mechanistic analysis.

[0078] In this embodiment, the experimental system utilizes a pre-trained BERT deep learning architecture to achieve automatic understanding and feature mapping of chemical reaction logic. The specific processing steps are as follows: Textual Encoding of Chemical Reactions: First, the experimental system transforms each purified chemical reaction event into a structured text entry. The SMILES strings of reactants and products are used as basic vocabulary and concatenated into a reaction sequence (e.g., [CLS]H2.CO [SEP] COH2[SEP]) using special connectors. This mapping method transforms the abstract molecular topological evolution into a semantic sequence recognizable by Natural Language Processing (NLP).

[0079] Bidirectional attention mechanism: The core of the BERT deep learning architecture lies in its multi-layer Transformer encoder module. Unlike traditional unidirectional models, this architecture uses a bidirectional self-attention mechanism to simultaneously capture the contextual relationships between chemical entities in the reaction sequence.

[0080] In the Fischer-Tropsch synthesis scenario: the Transformer model can automatically identify CO molecules in the reactants and surface carbides C in the products. ∗ There is a very strong semantic association between them (i.e., a "precursor-product" relationship), and this path is given a higher weight.

[0081] Hidden Layer Semantic Vector Extraction: After the text sequence is input into the Transformer model, it undergoes multiple nonlinear transformations. Finally, the [CLS] bit vector representing global features in the penultimate layer of the Transformer model is extracted as the semantic fingerprint of the reaction event. This vector is a high-dimensional dense vector (e.g., 768-dimensional), which not only contains species composition information but also implicitly encodes deep semantic features such as bond breaking type and chemical environment changes.

[0082] Self-supervised learning and pre-training: The Transformer model underwent large-scale self-supervised pre-training on a library of tens of millions of known chemical reaction equations (such as the USPTO dataset), mastering the basic chemical transformation logic. Therefore, when faced with the evolution of Fischer-Tropsch synthesis intermediates that were not seen in this embodiment, the Transformer model can still accurately identify potential hydrogenation, coupling, and other mechanism patterns thanks to its powerful transfer learning capabilities.

[0083] Through the semantic mapping of the BERT deep learning architecture described above, this invention transforms the originally discrete and fragmented atomic motion trajectories into a computationally computable and clusterable reaction semantic map within a low-dimensional manifold space.

[0084] In summary, the above statistical results (covering the main reaction pathways) and Figure 3 The reaction semantic clusters in the data exhibit a high degree of consistency, verifying the invention's ability to automatically analyze and classify different chemical mechanisms in complex catalytic networks.

[0085] S5. Transition state path generation: For each type of chemical reaction event with the same reaction mechanism, a representative conformation is selected, the influence of periodic boundary conditions is removed, and an adaptive simulation box with a vacuum layer is generated according to the actual size of the molecular cluster. An interpolation algorithm is used to generate an intermediate configuration sequence between the initial state of the reactant entity and the final state of the new product entity, and the output is a structure file.

[0086] Optionally, the selected representative conformation is preprocessed, the preprocessing including: cross-boundary stitching: using the inverse unit cell matrix transformation method to detect and repair broken molecules crossing periodic boundaries, restoring the representative conformation to a geometrically continuous molecular entity; adaptive unit cell construction: calculating the minimum bounding box of the reaction center atomic clusters of the representative conformation, and based on the minimum bounding box, extending the vacuum layer by at least 10 Å in each of the three-dimensional directions to construct an independent adaptive orthogonal simulated unit cell; centroid alignment: calculating the geometric centroids of the initial and final states of the species, and aligning the geometric centroids of the initial and final states of the species with the geometric center of the adaptive orthogonal simulated unit cell through translation transformation.

[0087] Optionally, the output of the structure file in step S5 includes: outputting a multi-frame trajectory file in XYZ format containing key time frames before, during, and after the chemical reaction, and each frame contains the original frame number and reaction ID information in the original molecular dynamics trajectory; outputting a continuous configuration sequence generated by image-dependent potential interpolation algorithm, wherein the continuous configuration sequence is used by quantum chemistry software to search for the initial path guess of the transition state of the elementary reaction.

[0088] In this embodiment, representativeness screening is performed as follows: For each type of reaction event with the same reactants and products (same SMILES), a median screening strategy is adopted. The system calculates the geometric radius of gyration for each data sample and selects the data sample with the most typical geometric conformation and the least affected by boundary conditions as the representative conformation. For example, in the generated reaction_events_clean.csv data, the calculated radius of gyration of event RXN-0007 (Frame 169) is at the median, so it is marked as the representative conformation, and a valid path file is generated under the directory 02_Structures / auto_neb_inputs / RXN-0007.... For subsequent repeated reactions of the same type (such as RXN-0011, RXN-0019, etc.), their status is marked as "Skipped (Representative Exists)," automatically skipping redundant calculations, thereby saving storage space and subsequent DFT calculation resources.

[0089] Adaptive simulation box optimization: For the selected representative reactions, the influence of periodic boundary conditions is automatically removed, and an adaptive simulation box with a vacuum layer is automatically generated based on the actual size of the molecular clusters, placing the clusters in the center. In this embodiment, an adaptive simulation box containing a 10 angstrom vacuum layer is established.

[0090] Nonlinear interpolation: The Image Dependent Pair Potential (IDPP) algorithm is used to generate a series of smoothly transitioning intermediate structure images between the initial state of the reactant entity and the final state of the new product entity, avoiding nuclear overlap. In this embodiment, the IDPP algorithm generates a transitional initial state prediction path file containing 8 interpolated images, which can be directly used in VASP software for NEB calculation.

[0091] Although embodiments of the invention have been described in conjunction with the accompanying drawings, those skilled in the art can make various modifications and variations without departing from the spirit and scope of the invention, and such modifications and variations all fall within the scope defined by the appended claims.

Claims

1. A molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification, characterized in that, The method includes: S1: Based on a zero-copy memory sharing and binary stream compression architecture, read molecular dynamics trajectory data from all time frames, extract the three-dimensional coordinates of all atoms from each time frame, calculate the interatomic distance, and determine the initial bonding relationship according to the preset bond length cutoff value; use a hidden Markov model for temporal smoothing, and introduce a substrate smart mask to filter the environmental background connection, generating a dynamic molecular connection graph. S2: Identify stable plateau periods in all time frames, compare the dynamic molecular connectivity graphs of adjacent time frames within the stable plateau period, and identify reactants and new products after chemical reactions; construct an atomic mapping bipartite graph and calculate connected components to isolate independent reaction events, and based on the element conservation filtering algorithm of chemical formula analysis, remove non-physical correlation fluxes between reactants and products that have no common elements to obtain purified chemical reaction events. S3: Arrange the chemical reaction events in chronological order to obtain an event list. Traverse the event list and, by identifying identity reactions, instantaneous reversible noise, and atomic index overlap, remove redundant chemical reaction events with stability below a preset stability threshold to obtain the final list of chemical reaction events. S4: The list of chemical reaction events is converted into a text sequence, semantic feature vectors are extracted using a neural network, and a reaction semantic map is constructed using a dimensionality reduction algorithm; at the same time, local structural features of the reaction center atoms in each event are extracted using geometric descriptors to construct a structural evolution landscape map; in the low-dimensional landscape map space, events with the same reaction mechanism are automatically clustered. S5: For each type of chemical reaction event with the same reaction mechanism, select a representative conformation, remove the influence of periodic boundary conditions, and generate an adaptive simulation box with a vacuum layer; use a nonlinear interpolation algorithm to generate a nonlinear interpolation path between the initial and final states of the species participating in the chemical reaction event, and output it as a structure file.

2. The molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification according to claim 1, characterized in that, The timing smoothing process described in step S1 includes: A hidden Markov model is constructed. For atom pairs that meet the preset bond length truncation condition within at least one time frame, the connection state sequence of the atom pairs in different time frames is taken as the observation sequence. The optimal hidden state sequence is calculated using the Viterbi algorithm. When the atom pair is in a non-bonding state, if an isolated bonding signal with a continuous frame number lower than a preset threshold appears in the observation sequence, it is identified as physical vibration noise and removed by the Viterbi optimal path.

3. The molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification according to claim 1, characterized in that, The substrate smart mask filtering environment background connection mentioned in step S1 includes: During the topology analysis initialization phase, a one-dimensional Boolean mask vector is calculated and generated based on the preset catalyst element type or atom index number. During the atom adjacency search process in each frame, bitwise operation filtering is performed using the one-dimensional Boolean mask vector to automatically skip bonding determination requests involving the interior of the substrate atoms and their interfaces.

4. The molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification according to claim 1, characterized in that, The element conservation filtering algorithm described in step S2 includes: Analyze the Hill representation chemical formulas of reactant and product entities and extract the element sets; determine if the intersection of the reactant and product sets is empty; if the intersection is empty, forcibly disconnect the associated edges corresponding to the intersection of the reactant and product sets in the atom mapping bipartite graph.

5. The molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification according to claim 1, characterized in that, Step S3 also includes: Using a preset adsorption cutoff radius scaling factor, the nearest mirror image distance between reactant or product molecules and substrate atoms is calculated in each time frame. If the nearest mirror image distance is less than the adsorption cutoff radius scaling factor, the corresponding species is marked as an adsorbed species; otherwise, it is marked as a gaseous species. By comparing the changes in species states before and after the chemical reaction, the species are automatically identified and classified.

6. The molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification according to claim 1, characterized in that, In step S5, the selected representative conformations are preprocessed, and the preprocessing includes: Cross-boundary stitching: The inverse unit cell matrix transformation method is used to detect and repair broken molecules that cross periodic boundaries, and the representative conformation is reduced to a geometrically continuous molecular entity. Adaptive cell construction: Calculate the minimum bounding box of the reaction center atomic clusters of representative conformations, and construct independent adaptive orthogonal simulation cells by extending a vacuum layer of at least 10 Å in each of the three-dimensional directions based on the minimum bounding box. Centroid Alignment: Calculate the geometric centroids of the initial and final states of the species, and then use a translation transformation to make the geometric centroids of the initial and final states of the species coincide with the geometric center of the adaptive orthogonal simulation cell.

7. The molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification according to claim 1, characterized in that, The text sequence extraction in step S4 includes concatenating the SMILES strings of reactants and new products into text entries; the neural network is a Transformer model based on a bidirectional self-attention mechanism, which extracts reaction semantic feature vectors through self-supervised learning; the geometric descriptor is a SOAP descriptor, used to capture the local three-dimensional configuration evolution within a 5 angstrom radius of the chemical reaction center; and the UMAP algorithm is used to map the high-dimensional semantic feature vectors to a low-dimensional landscape map space.

8. The molecular dynamics trajectory reaction identification method based on spatiotemporal decoupling and event purification according to claim 1, characterized in that, The output of the structure file in step S5 includes: The output contains a multi-frame trajectory file in XYZ format, which includes key time frames before, during, and after the chemical reaction. Each frame contains the original frame number and reaction ID information from the original molecular dynamics trajectory. The output contains a continuous sequence of configurations generated by image-dependent potential interpolation, which is used by quantum chemistry software to search for initial path guesses of elementary reaction transition states.

Images