Molecular dynamics simulation atom bonding state recognition method and system based on conditional random field

Abstract

The invention provides a molecular dynamics simulation atomic bonding state identification method and system based on a conditional random field, and the method comprises the steps: obtaining trajectory data generated through molecular dynamics simulation, and the trajectory data comprises atomic coordinates, atomic types and charge information which change with time; preprocessing the obtained trajectory data, screening a target atom set, and identifying potential bonding atom pairs frame by frame; based on the bonded atom pairs, extracting physical characteristics of the bonded atom pairs in each frame, including an inter-atomic Euclidean distance and a charge product; mapping the physical characteristics into a discrete observation state sequence; constructing a linear chain conditional random field model, and initializing model parameters to define conditional probability distribution between an observation state sequence and a hidden real keying state sequence; and inputting the observation state sequence into a linear chain conditional random field model, solving a hidden state sequence enabling the global probability to be maximum under conditional probability distribution by utilizing a Viterbi decoder, and obtaining an optimal keying state sequence of the atom pair changing along with time.

Description

Technical Field

[0001] This invention relates to the fields of computational chemistry, materials genome engineering, and cheminformatics, and more particularly to a method and system for identifying atomic bonding states in molecular dynamics simulations based on conditional random fields. Background Technology

[0002] Molecular dynamics simulations, as a crucial tool for studying the behavior of matter at the atomic scale, trace the spatiotemporal evolution trajectory of each atom in a system by solving Newton's equations of motion. This simulation method can accurately capture changes in the potential energy of interatomic interactions, providing crucial data support for understanding microscopic mechanisms such as material phase transitions and chemical reaction kinetics. With the development of high-performance computing technology, modern molecular dynamics simulations can now handle complex systems containing millions of atoms and achieve long-term simulations on the nanosecond to microsecond scale. The resulting trajectory data contains complete temporal information on multidimensional physical quantities such as atomic coordinates, velocities, and energies. Among these, accurate identification of interatomic bonding states is a fundamental step in analyzing chemical reaction networks, calculating thermodynamic parameters, and predicting material properties.

[0003] Current technologies primarily employ geometric criteria based on fixed distance thresholds for bonding analysis. While computationally simple, this method suffers from significant drawbacks in practical applications. First, the static cutoff radius cannot accommodate distance fluctuations caused by atomic thermal motion. Especially under high temperature and pressure conditions, minute oscillations in interatomic spacing near the critical value can be misinterpreted as frequent bonding-breaking events, generating numerous spurious chemical bond change signals. Second, existing methods utilize a discretized frame analysis mode for trajectory data processing, completely severing the temporal continuity of the chemical bond evolution process, making it difficult to accurately identify key intermediate states such as metastable transition states. Furthermore, relying solely on geometric distance parameters ignores the influence of quantum chemical effects such as charge distribution and electron cloud overlap on bonding states, resulting in low accuracy in identifying special interactions such as hydrogen bonds and coordination bonds. As the scale and duration of simulation systems continue to increase, traditional serial processing algorithms also face bottlenecks at the engineering implementation level, including low computational efficiency and excessive memory consumption. Summary of the Invention

[0004] To address the technical problems of existing technologies, such as significant thermal fluctuation noise interference, lack of time continuity constraints, and low efficiency in processing massive amounts of data, this invention provides a method and system for identifying atomic bonding states based on molecular dynamics simulation using conditional random fields (CRFs). This invention transforms the bonding discrimination problem into a time series labeling problem, utilizing the global sequence modeling capability of CRFs and combining the dual physical characteristics of distance and charge to achieve accurate and robust identification of dynamic bonding processes. Furthermore, it overcomes the bottleneck of large-scale data processing by employing shared memory parallel technology.

[0005] The technical means employed in this invention are as follows:

[0006] A method for identifying atomic bonding states in molecular dynamics simulations based on conditional random fields, comprising:

[0007] S1. Obtain trajectory data generated by molecular dynamics simulation, including time-varying atomic coordinates, atom type, and charge information;

[0008] S2. The acquired trajectory data is preprocessed to filter the target atom set, and potential bonding atom pairs are identified frame by frame based on the spatial neighborhood search algorithm.

[0009] S3. Based on the identified bonding atom pairs, extract their physical features in each frame, including the Euclidean distance between atoms and the charge product;

[0010] S4. Map the extracted physical features into a discrete sequence of observation states;

[0011] S5. Construct a Linear-Chain CRF model and initialize the model parameters to define the conditional probability distribution between the observed state sequence and the hidden true bonding state sequence.

[0012] S6. Input the observed state sequence into the constructed linear chain conditional random field model, and use the Viterbi decoder to solve for the hidden state sequence that maximizes the global probability under the conditional probability distribution, thereby obtaining the optimal bonding state sequence of the atom pair as time changes.

[0013] Further, step S1 includes:

[0014] S11. Read the trajectory file generated by the molecular dynamics simulation to obtain the original frame sequence;

[0015] S12. Extract the unique identifier, atom type, three-dimensional spatial coordinates and charge information of all atoms frame by frame in chronological order from the original frame sequence to generate the atomic structure and charge information corresponding to each frame.

[0016] S13. Perform periodic boundary condition correction on the three-dimensional spatial coordinates to ensure that all atomic coordinates are within the same simulation space.

[0017] S14. Generate structured trajectory data by arranging the corrected coordinates, atom types, and charge information in chronological order.

[0018] Further, step S2 includes:

[0019] S21. Introduce an atom type mapping mechanism. Based on the preset chemical meaning rules, map the atom types in the trajectory data to target atom tags with reactive properties, eliminate solvent molecules or inert gas atoms that do not participate in the reaction, and screen out the key target atoms that participate in the reaction.

[0020] S22. Load the coordinate information of the target atom set into the shared memory area to construct a globally accessible zero-copy data view;

[0021] S23. Using the KD-Tree space partitioning algorithm or the Cell-List grid hash indexing algorithm, the maximum physical truncation distance is used as the search radius for the target atoms in each frame. Atom pairs with a distance less than the truncation radius within each frame are quickly retrieved to form candidate atom pairs.

[0022] S24. Based on the potential bonding chemistry rules, the candidate atom pairs are filtered by type combination, and the atom pairs that conform to the potential bonding rules are retained to obtain the potential bonding atom pairs corresponding to each frame.

[0023] Further, step S3 includes:

[0024] S31. Read the coordinates and charge information of the potential bonding atom pairs frame by frame;

[0025] S32. Define the distance threshold range, calculate the Euclidean distance between atoms, and classify the distance between atoms according to the statistical law of bond length distribution of different bond types. The categories include close distance, medium distance, long distance and truncation distance.

[0026] S33. Define the electrostatic interaction range, calculate the charge product of atomic pairs to characterize the Coulomb interaction, and classify the charge product of atomic pairs according to the sign and magnitude of the product. The categories include strong attraction, weak attraction, neutral, weak repulsion and strong repulsion.

[0027] Further, step S4 includes:

[0028] S41. Assign corresponding hierarchical scores to the categories of interatomic distance. S dist The categories of the charge product of atomic pairs are assigned corresponding graded scores. S charge ;

[0029] S42, Set Distance Weights w d and charge weight w c ;

[0030] S43, Based on graded scores S dis Tiered scoringS charge Distance weight w d and charge weight w c Design a linear scoring function as follows:

[0031] ;

[0032] The linear scoring function maps continuous physical quantities to discrete observation labels, thereby generating an observation state sequence corresponding to each frame, including bonding tendency observations, bond breaking tendency observations, or fuzzy observations.

[0033] Further, step S5 includes:

[0034] S51. Define the hidden state space: Set the hidden state to contain only two discrete states: bonded and unbonded.

[0035] S52. Construct node feature functions: For each frame t, establish node feature functions between the observed state and the hidden state to describe the relationship between the observed state and the hidden state, that is, the support of the current observation for the current state.

[0036] S53. Constructing edge feature functions: For adjacent frames, establish edge feature functions for hidden state transitions to describe the relationship between adjacent hidden state transitions, i.e., the tendency of state transitions.

[0037] S54. Initialize model parameters: Based on prior knowledge of physicochemical processes, initialize the weights of the node feature functions and edge feature functions to determine the conditional probability distribution between the observed state sequence and the hidden bonding state sequence; where:

[0038] For node features, strong bonding observation signals are assigned positive weights pointing to the bonding state;

[0039] For edge features, the state is given a high positive weight, while the state is given a low or negative weight when it jumps frequently, in order to smooth out thermal fluctuation noise.

[0040] Further, step S6 includes:

[0041] S61. Input the discrete observation state sequence into the Viterbi decoder;

[0042] S62. For the first frame, calculate the initial scores for the two hidden states of unbonded and bonded states respectively. The scores are determined by the node feature function and its weights.

[0043] S63. For each subsequent frame, combine the cumulative score of the previous frame, the observation value of the current frame, the node feature function weight, and the edge feature function weight, calculate and retain the optimal state of the previous frame that maximizes the cumulative score and its corresponding index of the optimal state of the previous frame.

[0044] S64. After reaching the last frame, select the hidden state with the highest cumulative score as the endpoint, and backtrack back to the first frame along the index of the best state of the previous frame to form a complete hidden state sequence.

[0045] S65. Use the hidden state sequence obtained by backtracking as the bonding determination result of the atom pair in each frame to generate the optimal bonding state sequence that changes over time.

[0046] Furthermore, after step S6, the following steps are also included:

[0047] S7. Perform sparse matrix storage and statistical analysis on the optimal bonding state sequence.

[0048] Further, step S7 includes:

[0049] S71. Store the frame-by-frame binary bonding state matrix in a compressed sparse row format, recording only the non-zero elements in the bonding state and their row and column indices.

[0050] S72. Traverse the sparse matrix frame by frame, calculate the ratio of the number of bonds of each type of chemical bond in each frame to the total number of candidate pairs, and obtain the evolution curve of the bonding rate over time.

[0051] S73. Scan the continuous bonding segments in the sparse matrix, count the number of frames for each segment, and generate a bond lifetime distribution histogram.

[0052] S74. Accumulate the number of bonding events according to the type of atom pair and draw a bonding frequency heatmap to identify highly active bond types;

[0053] S75. Save the bonding rate curve, bond lifetime distribution histogram, and bonding frequency heatmap to a non-volatile storage medium to complete the statistical analysis.

[0054] This invention also provides a molecular dynamics simulation atomic bonding state identification system based on the above-mentioned molecular dynamics simulation atomic bonding state identification method, comprising:

[0055] The trajectory acquisition module is used to acquire trajectory data generated by molecular dynamics simulations, including time-varying atomic coordinates, atom types, and charge information;

[0056] The preprocessing module is used to preprocess the acquired trajectory data, filter the target atom set, and identify potential bonding atom pairs frame by frame based on the spatial neighborhood search algorithm.

[0057] The feature extraction module is used to extract the physical features of the identified bonding atom pairs in each frame, including the Euclidean distance between atoms and the charge product.

[0058] A mapping module is used to map the extracted physical features into a discrete sequence of observed states;

[0059] The model building and parameter initialization module is used to build a linear-chain conditional random field (CRF) model and initialize the model parameters to define the conditional probability distribution between the observed state sequence and the hidden true bonding state sequence.

[0060] The decoding module is used to input the observed state sequence into the constructed linear chain conditional random field model, and use the Viterbi decoder to solve for the hidden state sequence that maximizes the global probability under the conditional probability distribution, thereby obtaining the optimal bonding state sequence of the atom pair as time changes.

[0061] The results analysis module is used to perform sparse matrix storage and statistical analysis on the optimal bonding state sequence.

[0062] Compared with the prior art, the present invention has the following advantages:

[0063] 1. By introducing the state transition characteristics of CRF, this invention can effectively filter out spurious bonding / bonding signals caused by thermal fluctuations under high-temperature conditions. Experiments show that in 3000K high-temperature simulation data, the bond lifetime distribution identified by the method of this invention is more in line with Arrhenius's law, eliminating a large number of non-physical glitches with extremely short lifetimes (<100fs) found in traditional methods.

[0064] 2. This invention no longer relies solely on distance, but incorporates charge information. This is particularly important for ionic or polar covalent bond systems, as it avoids misjudging spatially close but charge-repelling atomic pairs as bonds.

[0065] 3. This invention employs shared memory technology to avoid data serialization and communication overhead between multiple processes, thereby improving data reading and preprocessing speed by 3-5 times. Combined with sparse matrix storage, it makes it possible to process tens of millions of frames of fully atomic trajectories on ordinary workstations.

[0066] In summary, the framework of this invention is not limited to specific atom types. By adjusting the threshold parameters in the feature extraction module, it can be widely applied to various MD simulation scenarios such as metal oxidation, organic synthesis, and biomacromolecule folding. Attached Figure Description

[0067] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the 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 based on these drawings without creative effort.

[0068] Figure 1 This is a schematic diagram of the overall process and code architecture of the method of the present invention.

[0069] Figure 2 This is a timing flowchart of the parallel computing and CRF analysis of the present invention.

[0070] Figure 3 This is a schematic diagram of the CRF model structure and Viterbi decoding principle in this invention.

[0071] Figure 4 This is a class diagram of the core functional modules of the system of the present invention. Detailed Implementation

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

[0073] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0074] like Figure 1 As shown, this invention provides a method for identifying atomic bonding states based on molecular dynamics simulation using conditional random fields, comprising:

[0075] S1. Obtain trajectory data generated by molecular dynamics simulation, including time-varying atomic coordinates, atom type, and charge information;

[0076] S2. The acquired trajectory data is preprocessed to filter the target atom set, and potential bonding atom pairs are identified frame by frame based on the spatial neighborhood search algorithm.

[0077] S3. Based on the identified bonding atom pairs, extract their physical features in each frame, including the Euclidean distance between atoms and the charge product;

[0078] S4. Map the extracted physical features into a discrete sequence of observation states;

[0079] S5. Construct a linear chain conditional random field model and initialize the model parameters to define the conditional probability distribution between the observed state sequence and the hidden true bonding state sequence.

[0080] S6. Input the observed state sequence into the constructed linear chain conditional random field model, and use the Viterbi decoder to solve for the hidden state sequence that maximizes the global probability under the conditional probability distribution, thereby obtaining the optimal bonding state sequence of the atom pair as time changes.

[0081] In a specific implementation, as a preferred embodiment of the present invention, step S1 includes:

[0082] S11. Read the trajectory file generated by the molecular dynamics simulation (such as the LAMMPS dump format or the GROMACS xtc format) to obtain the original frame sequence;

[0083] S12. Extract the unique identifier, atom type, three-dimensional spatial coordinates and charge information of all atoms frame by frame in chronological order from the original frame sequence to generate the atomic structure and charge information corresponding to each frame.

[0084] S13. Perform periodic boundary condition correction on the three-dimensional spatial coordinates to ensure that all atomic coordinates are within the same simulation space.

[0085] S14. Generate structured trajectory data by arranging the corrected coordinates, atom types, and charge information in chronological order.

[0086] In a specific implementation, as a preferred embodiment of the present invention, step S2 includes:

[0087] S21. To reduce memory consumption, an atom type mapping mechanism is introduced. Based on the preset chemical meaning rules, the atom types in the trajectory data are mapped to target atom tags with reactive properties. Solvent molecules or inert gas atoms that do not participate in the reaction are eliminated, and key target atoms that participate in the reaction (such as core metal atoms, surface active site atoms, environmental oxidant atoms, etc.) are screened out.

[0088] S22. Load the coordinate information of the target atom set into the shared memory area to construct a globally accessible zero-copy data view; In this embodiment, considering the large size of the MD trajectory data file, the present invention designs a parallel I / O architecture based on Python multiprocessing and shared memory, such as... Figure 2 As shown, multi-process parallel processing and shared memory technology are adopted, namely: the main process loads trajectory data into the system's shared memory area; multiple worker processes directly access trajectory data slices through shared memory and execute atomic compression, filtering, and atomic pair discovery tasks in parallel; each worker process summarizes the processing results to the main process to build a global atomic pair index.

[0089] S23. Using the KD-Tree spatial partitioning algorithm or the Cell-List grid hash indexing algorithm, the maximum physical cutoff distance is used as the search radius for target atoms in each frame. Atom pairs with a distance less than the cutoff radius within each frame are quickly retrieved to form candidate atom pairs. In this embodiment, to avoid performing fully connected computations on all N(N-1) / 2 atom pairs, a coarse-screening strategy based on the maximum physical cutoff distance is implemented. In each frame, the spatial neighborhood is quickly searched using the KD-Tree or Cell-List algorithm, retaining only atom pairs with a distance less than a broad threshold (e.g., 3.5 Å or 4.0 Å) and atom type combinations that conform to potential bonding rules (e.g., metal-oxygen, metal-fluorine, etc.). These selected atom pairs are defined as candidate atom pairs, and their frame index sequence is recorded to form a sparse adjacency list structure, thereby reducing the subsequent computational complexity from O(N) to O(N) / 2. 2 The value is reduced to the O(N) level.

[0090] S24. Based on the potential bonding chemistry rules, the candidate atom pairs are filtered by type combination, and the atom pairs that conform to the potential bonding rules are retained to obtain the potential bonding atom pairs corresponding to each frame.

[0091] In a specific implementation, as a preferred embodiment of the present invention, step S3 includes:

[0092] S31. Read the coordinates and charge information of the potential bonding atom pairs frame by frame;

[0093] S32. Define the distance threshold range, calculate the Euclidean distance between atoms, and classify the distance between atoms according to the statistical law of bond length distribution of different bond types (such as Al-O bond and Al-F bond). The categories include close distance, medium distance, long distance and truncation distance.

[0094] S33. Define the electrostatic interaction region and calculate the charge product of atomic pairs.q i q j The product of charges is used to characterize Coulomb interactions. Based on the sign and magnitude of the product, the product of charges of atomic pairs is classified into categories including strong attraction, weak attraction, neutral, weak repulsion, and strong repulsion.

[0095] In a specific implementation, as a preferred embodiment of the present invention, step S4 includes:

[0096] S41. Assign corresponding hierarchical scores to the categories of interatomic distance. S dist The categories of the charge product of atomic pairs are assigned corresponding graded scores. S charge ;

[0097] S42, Set Distance Weights w d and charge weight w c ;

[0098] S43, Based on graded scores S dis Tiered scoring S charge Distance weight w d and charge weight w c Design a linear scoring function as follows:

[0099] ;

[0100] The linear scoring function maps continuous physical quantities to discrete observation labels (e.g., label "1" represents a strong bonding signal, label "2" represents a fuzzy / transitional signal, and label "0" represents a bond breaking signal), thereby generating an observation state sequence corresponding to each frame, including bonding tendency observations, bond breaking tendency observations, or fuzzy observations.

[0101] In this embodiment, a comprehensive score is calculated by combining distance classification and charge classification, and the current frame is marked as a bonding tendency observation, a bond breaking tendency observation, or a fuzzy observation based on the comprehensive score.

[0102] In a specific implementation, as a preferred embodiment of the present invention, step S5 includes:

[0103] S51. Define the hidden state space: Set the hidden state to contain only two discrete states: bonded and unbonded.

[0104] S52. Construct node feature functions: such as Figure 3As shown, for each frame t, a node feature function f_k(t, y_t, x_t) is established between the observed state x_t and the hidden state y_t to describe the relationship between the observed state and the hidden state, that is, the support of the current observation for the current state.

[0105] S53. Constructing edge feature functions: such as Figure 3 As shown, for adjacent frames (t−1, t), a side feature function g_m(y_t−1, y_t) for hidden state transition is established to describe the relationship between adjacent hidden states, i.e., the tendency of state transition;

[0106] S54. Initialize model parameters: Based on prior knowledge of physicochemical processes, initialize the weights of the node feature functions and edge feature functions to determine the conditional probability distribution between the observed state sequence and the hidden bonding state sequence; where:

[0107] For node features, strong bonding observation signals are assigned positive weights pointing to bonding states; for edge features, state persistence (bonding-bonding, bond breaking-bonding) is assigned higher positive weights, while frequent state switching (bonding-bonding) is assigned lower or negative weights to smooth out thermal fluctuation noise. In this embodiment, to suppress thermal fluctuation noise, the invention specifically assigns higher positive weights to state persistence (i.e., 1→1 or 0→0), while assigning lower or even negative weights to state switching (i.e., 1→0 or 0→1). This design introduces a hysteresis loop effect similar to that in physics into the probabilistic model, ensuring that bonding states only flip when the observed signal undergoes a sustained and significant change.

[0108] In a specific implementation, as a preferred embodiment of the present invention, step S6 includes:

[0109] S61. Input the discrete observation state sequence into the Viterbi decoder;

[0110] S62. For the first frame, calculate the initial scores for the two hidden states of unbonded and bonded states respectively. The scores are determined by the node feature function and its weights.

[0111] S63. For each subsequent frame, combine the cumulative score of the previous frame, the observation value of the current frame, the node feature function weight, and the edge feature function weight, calculate and retain the optimal state of the previous frame that maximizes the cumulative score and its corresponding index of the optimal state of the previous frame.

[0112] S64. After reaching the last frame, select the hidden state with the highest cumulative score as the endpoint, and backtrack back to the first frame along the index of the best state of the previous frame to form a complete hidden state sequence.

[0113] S65. Use the hidden state sequence obtained by backtracking as the bonding determination result of the atom pair in each frame to generate the optimal bonding state sequence that changes over time.

[0114] In this embodiment, as Figure 3 As shown, the Viterbi algorithm from the dynamic programming domain is used to find the optimal path among all possible hidden state paths that maximizes the weighted sum of the global feature functions (i.e., maximizes the global probability). Unlike traditional greedy algorithms, this step does not make a judgment based solely on the information of the current frame, but comprehensively considers all observation information from the past, present, and future. Even if, at a certain instant (within a few frames), the distance between atoms temporarily increases due to thermal motion, as long as the bonding signal in the preceding and following context remains strong, the transition feature constraints of the CRF model will automatically correct this abnormal state, smoothing it into a bonding state, thereby effectively removing flicker noise.

[0115] In a specific implementation, as a preferred embodiment of the present invention, after step S6, the following is also included:

[0116] S7. Perform sparse matrix storage and statistical analysis on the optimal bonding state sequence.

[0117] In a specific implementation, as a preferred embodiment of the present invention, step S7 includes:

[0118] S71. Store the frame-by-frame binary bonding state matrix in a compressed sparse row format, recording only the non-zero elements in the bonding state and their row and column indices.

[0119] S72. Traverse the sparse matrix frame by frame, calculate the ratio of the number of bonds of each type of chemical bond in each frame to the total number of candidate pairs, and obtain the evolution curve of the bonding rate over time.

[0120] S73. Scan the continuous bonding segments in the sparse matrix, count the number of frames for each segment, and generate a bond lifetime distribution histogram.

[0121] S74. Accumulate the number of bonding events according to the type of atom pair and draw a bonding frequency heatmap to identify highly active bond types;

[0122] S75. Save the bonding rate curve, bond lifetime distribution histogram, and bonding frequency heatmap to a non-volatile storage medium to complete the statistical analysis.

[0123] In this embodiment, the decoding output of the bonding state sequence is a huge binary matrix. Considering the sparsity of chemical reactions (most atom pairs are not connected most of the time), this invention uses a compressed sparse row (CSR) format to store the results, compressing the storage space by several orders of magnitude. Based on this matrix, the system further calculates the bond lifetime distribution, bonding rate evolution curve, and reaction network topology.

[0124] Corresponding to the molecular dynamics simulation atomic bonding state identification method in this application, this application also provides a molecular dynamics simulation atomic bonding state identification system, including:

[0125] The trajectory acquisition module is used to acquire trajectory data generated by molecular dynamics simulations, including time-varying atomic coordinates, atom types, and charge information;

[0126] The preprocessing module is used to preprocess the acquired trajectory data, filter the target atom set, and identify potential bonding atom pairs frame by frame based on the spatial neighborhood search algorithm.

[0127] The feature extraction module is used to extract the physical features of the identified bonding atom pairs in each frame, including the Euclidean distance between atoms and the charge product.

[0128] A mapping module is used to map the extracted physical features into a discrete sequence of observed states;

[0129] The model building and parameter initialization module is used to build a linear chain conditional random field model and initialize the model parameters to define the conditional probability distribution between the observed state sequence and the hidden true bonding state sequence.

[0130] The decoding module is used to input the observed state sequence into the constructed linear chain conditional random field model, and use the Viterbi decoder to solve for the hidden state sequence that maximizes the global probability under the conditional probability distribution, thereby obtaining the optimal bonding state sequence of the atom pair as time changes.

[0131] The results analysis module is used to perform sparse matrix storage and statistical analysis on the optimal bonding state sequence.

[0132] The embodiments of the present invention are described simply because they correspond to those in the embodiments above. For any similarities, please refer to the descriptions in the embodiments above, which will not be elaborated here.

[0133] Example

[0134] This embodiment uses the reaction molecular dynamics simulation data of an energetic composite material of aluminum nanoparticles (ANP) and polyvinylidene fluoride (PVDF) as an example to demonstrate the complete workflow of the system.

[0135] This system is developed based on a Python environment. Its core dependencies include NumPy for numerical computation, Pandas for data structure operations, SciPy.sparse for sparse matrix storage, and Multiprocessing for parallel computing. The recommended hardware configuration is a multi-core CPU workstation (32 cores or more) with at least 128GB of memory. Due to the use of shared memory technology, memory size is a key bottleneck for processing large-scale trajectories.

[0136] System software architecture such as Figure 4 As shown, it consists of the following core class modules:

[0137] StepwiseAnalyzer: As the main control module, it is responsible for scheduling the entire pipeline and managing the dependencies between each step.

[0138] DataCompressor: Responsible for cleaning, type mapping, and compression of trajectory data.

[0139] FeatureExtractor: Responsible for calculating the geometric and charge characteristics of atomic pairs.

[0140] BondingStateManager: Responsible for managing atomic pair indexes, state matrix storage, and CSR transformations.

[0141] PureCRFProcessor (CRF Processor): Built-in AtomPairCRF class, implementing feature function definition and Viterbi decoding algorithm.

[0142] OutputManager: Responsible for logging, intermediate result caching (HDF5 / Pickle), and final chart output.

[0143] Based on the architecture described above, the complete workflow of the system is as follows:

[0144] Step 1: Parallel data preprocessing based on shared memory:

[0145] When processing large-scale MD trajectories (such as dump.lammpstrj), file I / O and memory usage are the main bottlenecks. This embodiment employs the following optimization strategies:

[0146] Atom Type Mapping: Fine-grained atom classification is defined through the AtomTypeMapper class. For example, in an aluminum / fluorine / oxygen system with core-shell structured nano-aluminum particles, the original atom type IDs (1, 2, 3...) are mapped to chemically meaningful labels, such as: CA1 (Core Aluminum) corresponds to type ID [1, 4], representing aluminum atoms located inside the particles; sAl (Surface Aluminum) corresponds to type ID [2, 5], representing active aluminum atoms located on the surface; sO (Surface Oxygen) corresponds to type ID [3, 6], representing oxygen atoms in the surface oxide layer; eO (Environmental Oxygen) corresponds to type ID [7], representing free oxygen in the environment. This mapping not only facilitates subsequent analysis but also allows background gases or irrelevant atoms to be directly filtered out during the preprocessing stage, significantly reducing the amount of data.

[0147] Shared memory pool construction: After the main process reads the trajectory file, it does not directly pass the huge list of frames to the child process. Instead, it creates a multiprocessing.shared_memory.SharedMemory block. The parsed trajectory data (containing id, type, x, y, z, q fields) is written to this memory block as a NumPy structured array.

[0148] Parallel compression workflow (Worker Function): The system starts multiple worker processes via ProcessPoolExecutor. Each worker receives a task packet containing a frame index range (e.g., frames 0-100) and the name of a shared memory block. The worker function _compress_frames_worker_mp performs the following operations:

[0149] Operation 1: Connect to an existing shared memory block by name.

[0150] Operation 2: Read a data slice of a specified frame directly from memory (Zero-copy read).

[0151] Operation 3: Use Boolean indexing to filter out the atomic rows belonging to cAl, sAl, sO, eO.

[0152] Operation 4: Output a lightweight DataFrame containing atomic categories and specific types.

[0153] This design avoids the huge overhead of Python inter-process communication (Pickle serialization), allowing preprocessing speed to increase linearly with the number of CPU cores.

[0154] Step 2: Candidate Pair Discovery Based on Distance Truncation

[0155] To construct the input for the CRF model, it is necessary to first determine which atom pairs are likely to form bonds. With hundreds of millions of atom pairs in the entire system, direct calculation is not feasible. This embodiment employs a frame-by-frame neighborhood search strategy, as follows:

[0156] Set a relatively wide physical cutoff radius (Max Cutoff), for example, 3.5 Å. This value is much larger than the conventional chemical bond length (about 1.8-2.0 Å) in order to ensure that all possible bonding events and the "connected strands" state after bond breakage are covered, and to avoid omissions.

[0157] In each frame, the KD-Tree algorithm is used to search for atom pairs with a distance of less than 3.5 Å.

[0158] Type filtering: Only combinations that conform to chemical logic are retained. For example, the system configuration defines a list of valid bond types: [('sO1-cAl1'), ('sO2-sAl2'), ('F-Al')...]. If a pair of atoms is found to be sAl and sAl (metallic bond, of interest in non-aluminum / fluorine / oxygen systems) or eO and eO, it is discarded directly in the current step.

[0159] Global index construction: BondingStateManager maintains a global dictionary `atom_pairs`, with keys being tuples of type `atom1_id`, `atom2_id`, and values ​​being unique integer indices. Each atomic pair has a unique ID regardless of how many frames it appears in. Simultaneously, a list of frames `frames_list` is recorded for each atomic pair, preparing for subsequent sparse computation.

[0160] Step 3: Multidimensional physical feature extraction and discretization

[0161] This step serves as a bridge connecting the physical world and the probabilistic model. The FeatureExtractor class calculates features for each candidate atom pair in each frame's state, specifically including:

[0162] Distance Categorization: Different chemical bonds have different equilibrium bond lengths, therefore the feature thresholds are dynamically configured. Taking the F-Al bond as an example, the system sets the threshold ranges as follows: For the very_close category, the range is 0 ≤ d < 1.8 Å; for the close category, the range is 1.8 ≤ d < 2.2 Å; for the medium category, the range is 2.2 ≤ d < 2.7 Å; for the far category, the range is 2.7 ≤ d < 3.5 Å; and for the beyond_cutoff category, the range is d ≥ 3.5 Å.

[0163] Charge categorization: Calculating the charge product Q=q 1 q 2 Based on the sign and magnitude of the product, the charge products of atomic pairs are categorized as follows: Q < -0.5 (a typical characteristic of ionic bonds) is classified as strong_attractive; -0.5 ≤ Q < -0.1 is classified as attractive; -0.1 ≤ Q < 0.1 is classified as neutral; 0.1 ≤ Q < 0.5 is classified as repulsive; and Q ≥ 0.5 is classified as strong_repulsive.

[0164] Integrated Scoring and Observation Mapping: To transform continuous features into discrete observations required by the CRF, the system defines the scoring logic:

[0165] ;

[0166] For example: very_close gets +5 points, far gets -1 point; strong_attractive gets +5 points, repulsive gets -1 point. Based on the total score, observation labels are generated. x t If the score is ≥ 5 (e.g., very close distance and strong charge attraction), mark it as 1 (strong bonding tendency). If the score is ≤ -2 (e.g., far distance or strong repulsion), mark it as '0' (strong bond breaking tendency). Other intermediate scores are marked as 2 (fuzzy / transitional state). For frames where candidate pairs do not appear (i.e., distance > 3.5 Å), directly fill the observation value with 0.

[0167] Step 4: Linear Chain CRF Model Construction and Viterbi Decoding (The core algorithm of this invention is implemented by the PureCRFProcessor and AtomPairCRF classes)

[0168] Linear chain CRF model construction: The model contains two hidden states, State=1 (bonded) and State=0 (not bonded). The weights of the potential functions are defined as follows (based on the get_feature_weight method):

[0169] Define node feature weights:

[0170] The weight of obs_bond_strong (observation is 1, state is 1) is defined as +5.0, which means that when the physical features strongly suggest bonding, the model will determine that bonding is highly probable.

[0171] The weight of obs_unbond_strong (observation is 0, state is 0) is defined as +5.0.

[0172] The weight of obs_bond_medium (observation 2, state 1) is defined as +2.0, indicating that the fuzzy signal has some support for bonding, but it is not as significant as the strong signal.

[0173] Define the weight of obs_bond_weak (observation is 0, state is 1) as -1.0. This means that the observation shows a broken bond, but the state is determined to be a bond. This is a penalty, but it is allowed under strong transition constraints (i.e., noise-resistant).

[0174] Define edge / transition weights:

[0175] The weight of transition_bond_bond(1→1) is defined as +2.5, which greatly rewards the continuity of the state.

[0176] The weight of transition_unbond_unbond(0→0) is defined as +2.5.

[0177] The weight of transition_bond_unbond(1↔0) is defined as -1.0, which penalizes the change of state.

[0178] This weighting configuration is mathematically equivalent to introducing an energy barrier. The Viterbi path will only switch states when the energy accumulated by a new observation signal (such as a series of break signals) is sufficient to offset the penalty for state transition (-1.0) and exceed the reward for state preservation (+2.5). This is the fundamental reason why this invention can eliminate flicker noise.

[0179] Viterbi Decoding: Observation Sequence for Each Atom Pair The decoding process is as follows:

[0180] Initialization (t=0): Calculate the scores for initial states of 0 and 1.

[0181] Recursion: For t=1 to T-1, calculate the maximum score path to the current state j.

[0182] ;

[0183] Simultaneously record the path pointer .

[0184] Backtracking: Starting from the state with the highest score at time T-1, trace the optimal path backward. The final sequence This is the optimal bonding state sequence for the atom pair throughout the entire simulation process.

[0185] Step 5: Sparse Storage and Result Analysis

[0186] After decoding, the system faces storage challenges. Assuming 10,000 candidate atom pairs and 300,000 frames, directly storing an int8 array would require approximately 3GB of memory. For larger systems, memory would quickly run out. This embodiment utilizes the BondingStateManager class to convert the results into a CSR (Compressed Sparse Row) matrix format. Because chemical bonds are sparse in both space and time (any atom typically bonds with only a few atoms, and the bonding time is finite), sparse matrices can reduce storage requirements to 1%-5% of the original. The final result is saved as an HDF5 file (.h5), containing:

[0187] data, indices, indptr: the three essential elements of a CSR matrix;

[0188] pair_info: A mapping table between atomic pair IDs and types.

[0189] Based on the decoded state matrix, the system automatically performs the following statistics:

[0190] Bonding Ratio: Counts the number of each type of bond (such as sAl-sO) in each frame and generates a curve that evolves over time to analyze the reaction rate.

[0191] Bond Lifetime Distribution: This involves scanning consecutive "1" segments in the state sequence. For example, the sequence 001111100 corresponds to a bond with a lifetime of 5 frames. The system statistically analyzes the frequency distribution of all bond lifetimes to distinguish between transient collisions (extremely short lifetimes) and steady-state chemical bonds (long lifetimes).

[0192] Repairing specific atom pairs: The code used by the method also includes debugging and repair logic for specific objects of interest, such as (2692, 6316), allowing users to manually intervene or delve into the reaction behavior of specific sites.

[0193] In summary, this invention, by constructing a complete data stream—from parallel reading, multi-feature fusion, CRF probabilistic inference to sparse analysis—achieves high-precision and high-efficiency identification of complex bonding processes in molecular dynamics simulations, providing a powerful analytical tool for the field of computational chemistry.

[0194] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for identifying atomic bonding states based on molecular dynamics simulation using conditional random fields, characterized in that, include: S1. Obtain trajectory data generated by molecular dynamics simulation, including time-varying atomic coordinates, atom type, and charge information; S2. The acquired trajectory data is preprocessed to filter the target atom set, and potential bonding atom pairs are identified frame by frame based on the spatial neighborhood search algorithm. S3. Based on the identified bonding atom pairs, extract their physical features in each frame, including the Euclidean distance between atoms and the charge product; S4. Map the extracted physical features into a discrete sequence of observation states; S5. Construct a linear chain conditional random field model and initialize the model parameters to define the conditional probability distribution between the observed state sequence and the hidden true bonding state sequence. S6. Input the observed state sequence into the constructed linear chain conditional random field model, and use the Viterbi decoder to solve for the hidden state sequence that maximizes the global probability under the conditional probability distribution, thereby obtaining the optimal bonding state sequence of the atom pair as time changes.

2. The method for identifying atomic bonding states based on conditional random fields in molecular dynamics simulations according to claim 1, characterized in that, Step S1 includes: S11. Read the trajectory file generated by the molecular dynamics simulation to obtain the original frame sequence; S12. Extract the unique identifier, atom type, three-dimensional spatial coordinates and charge information of all atoms frame by frame in chronological order from the original frame sequence to generate the atomic structure and charge information corresponding to each frame. S13. Perform periodic boundary condition correction on the three-dimensional spatial coordinates to ensure that all atomic coordinates are within the same simulation space. S14. Generate structured trajectory data by arranging the corrected coordinates, atom types, and charge information in chronological order.

3. The method for identifying atomic bonding states based on conditional random fields in molecular dynamics simulation according to claim 1, characterized in that, Step S2 includes: S21. Introduce an atom type mapping mechanism. Based on the preset chemical meaning rules, map the atom types in the trajectory data to target atom tags with reactive properties, eliminate solvent molecules or inert gas atoms that do not participate in the reaction, and screen out the key target atoms that participate in the reaction. S22. Load the coordinate information of the target atom set into the shared memory area to construct a globally accessible zero-copy data view; S23. Using the KD-Tree space partitioning algorithm or the Cell-List grid hash indexing algorithm, the maximum physical truncation distance is used as the search radius for the target atoms in each frame. Atom pairs with a distance less than the truncation radius within each frame are quickly retrieved to form candidate atom pairs. S24. Based on the potential bonding chemistry rules, the candidate atom pairs are filtered by type combination, and the atom pairs that conform to the potential bonding rules are retained to obtain the potential bonding atom pairs corresponding to each frame.

4. The method for identifying atomic bonding states based on conditional random fields in molecular dynamics simulation according to claim 1, characterized in that, Step S3 includes: S31. Read the coordinates and charge information of the potential bonding atom pairs frame by frame; S32. Define the distance threshold range, calculate the Euclidean distance between atoms, and classify the distance between atoms according to the statistical law of bond length distribution of different bond types. The categories include close distance, medium distance, long distance and truncation distance. S33. Define the electrostatic interaction range, calculate the charge product of atomic pairs to characterize the Coulomb interaction, and classify the charge product of atomic pairs according to the sign and magnitude of the product. The categories include strong attraction, weak attraction, neutral, weak repulsion and strong repulsion.

5. The method for identifying atomic bonding states based on conditional random fields in molecular dynamics simulations according to claim 1, characterized in that, Step S4 includes: S41. Assign corresponding hierarchical scores to the categories of interatomic distance. S dist The categories of the charge product of atomic pairs are assigned corresponding graded scores. S charge ; S42, Set Distance Weights w d and charge weight w c ; S43, Based on graded scores S dis Tiered scoring S charge Distance weight w d and charge weight w c Design a linear scoring function as follows: ； The linear scoring function maps continuous physical quantities to discrete observation labels, thereby generating an observation state sequence corresponding to each frame, including bonding tendency observations, bond breaking tendency observations, or fuzzy observations.

6. The method for identifying atomic bonding states based on conditional random fields in molecular dynamics simulations according to claim 1, characterized in that, Step S5 includes: S51. Define the hidden state space: Set the hidden state to contain only two discrete states: bonded and unbonded. S52. Construct node feature functions: For each frame, establish node feature functions between the observed state and the hidden state to describe the relationship between the observed state and the hidden state, that is, the support of the current observation for the current state. S53. Constructing edge feature functions: For adjacent frames, establish edge feature functions for hidden state transitions to describe the relationship between adjacent hidden state transitions, i.e., the tendency of state transitions. S54. Initialize model parameters: Based on prior knowledge of physicochemical processes, initialize the weights of the node feature functions and edge feature functions to determine the conditional probability distribution between the observed state sequence and the hidden bonding state sequence; where: For node features, strong bonding observation signals are assigned positive weights pointing to the bonding state; For edge features, the state is given a high positive weight, while the state is given a low or negative weight when it jumps frequently, in order to smooth out thermal fluctuation noise.

7. The method for identifying atomic bonding states based on conditional random fields in molecular dynamics simulation according to claim 1, characterized in that, Step S6 includes: S61. Input the discrete observation state sequence into the Viterbi decoder; S62. For the first frame, calculate the initial scores for the two hidden states of unbonded and bonded states respectively. The scores are determined by the node feature function and its weights. S63. For each subsequent frame, combine the cumulative score of the previous frame, the observation value of the current frame, the node feature function weight, and the edge feature function weight, calculate and retain the optimal state of the previous frame that maximizes the cumulative score and its corresponding index of the optimal state of the previous frame. S64. After reaching the last frame, select the hidden state with the highest cumulative score as the endpoint, and backtrack back to the first frame along the index of the best state of the previous frame to form a complete hidden state sequence. S65. Use the hidden state sequence obtained by backtracking as the bonding determination result of the atom pair in each frame to generate the optimal bonding state sequence that changes over time.

8. The method for identifying atomic bonding states based on conditional random fields in molecular dynamics simulation according to claim 1, characterized in that, After step S6, the following also includes: S7. Perform sparse matrix storage and statistical analysis on the optimal bonding state sequence.

9. The method for identifying atomic bonding states based on conditional random fields in molecular dynamics simulation according to claim 8, characterized in that, Step S7 includes: S71. Store the frame-by-frame binary bonding state matrix in a compressed sparse row format, recording only the non-zero elements in the bonding state and their row and column indices. S72. Traverse the sparse matrix frame by frame, calculate the ratio of the number of bonds of each type of chemical bond in each frame to the total number of candidate pairs, and obtain the evolution curve of the bonding rate over time. S73. Scan the continuous bonding segments in the sparse matrix, count the number of frames for each segment, and generate a bond lifetime distribution histogram. S74. Accumulate the number of bonding events according to the type of atom pair and draw a bonding frequency heatmap to identify highly active bond types; S75. Save the bonding rate curve, bond lifetime distribution histogram, and bonding frequency heatmap to a non-volatile storage medium to complete the statistical analysis.

10. A molecular dynamics simulation atomic bonding state identification system based on the molecular dynamics simulation atomic bonding state identification method according to any one of claims 1-9, characterized in that, include: The trajectory acquisition module is used to acquire trajectory data generated by molecular dynamics simulations, including time-varying atomic coordinates, atom types, and charge information; The preprocessing module is used to preprocess the acquired trajectory data, filter the target atom set, and identify potential bonding atom pairs frame by frame based on the spatial neighborhood search algorithm. The feature extraction module is used to extract the physical features of the identified bonding atom pairs in each frame, including the Euclidean distance between atoms and the charge product. A mapping module is used to map the extracted physical features into a discrete sequence of observed states; The model building and parameter initialization module is used to build a linear chain conditional random field model and initialize the model parameters to define the conditional probability distribution between the observed state sequence and the hidden true bonding state sequence. The decoding module is used to input the observed state sequence into the constructed linear chain conditional random field model, and use the Viterbi decoder to solve for the hidden state sequence that maximizes the global probability under the conditional probability distribution, thereby obtaining the optimal bonding state sequence of the atom pair as time changes. The results analysis module is used to perform sparse matrix storage and statistical analysis on the optimal bonding state sequence.

Images