Molecular dynamics fast simulation method and system for lithium-silicon alloy

Abstract

This application relates to the field of molecular dynamics simulation technology, and particularly to a method and system for rapid molecular dynamics simulation of lithium-silicon alloys. The method includes the following steps: constructing an initial structural model of the lithium-silicon alloy to be tested based on a preset lithium-silicon ratio and lattice constant; constructing an adaptive cutoff radius lookup table for atomic pairs; performing accelerated reactive molecular dynamics simulation on the initial structural model of the lithium-silicon alloy under preset temperature and simulation conditions using the adaptive cutoff radius lookup table for atomic pair adaptation and evolution, and obtaining the structural data of the lithium-silicon alloy after the simulation; and outputting the atomic topological evolution process, bond order characteristics, or material stability of the lithium-silicon alloy to be tested based on the structural data, thereby eliminating invalid atomic pair calculations in the lithium-silicon alloy simulation and improving the computational efficiency of large-scale reactive molecular dynamics simulations while ensuring the accuracy of bond order calculations and the physical correctness of the reaction process.

Description

Technical Field

[0001] This application relates to the field of molecular dynamics simulation technology, and in particular to a rapid molecular dynamics simulation method and system for lithium-silicon alloys. Background Technology

[0002] Molecular dynamics (MD) simulations are computational methods that study the structural evolution and dynamic behavior of materials at the atomic scale by numerically integrating atomic equations of motion. This method can reveal structural changes and reaction mechanisms at the microscopic level without relying on macroscopic empirical models, and has been widely applied in materials science, chemical reaction kinetics, and biomolecular simulations. In molecular dynamics simulations, the selection of the interatomic interaction potential function directly determines the accuracy and applicability of the simulation. With the improvement of computational power, researchers are constantly increasing their requirements for the scale and timescale of simulation systems.

[0003] Traditional classical force fields (such as Lennard-Jones potential and EAM potential) are usually based on fixed topological structures or simplified forms of many-body interactions, making it difficult to accurately describe reaction processes involving chemical bond breaking and formation. To address this issue, reactive force fields (ReaxFF) have been proposed and are widely used in reactive molecular dynamics simulations in fields such as battery materials, combustion chemistry, and catalytic reactions.

[0004] However, ReaxFF requires calculating the bond order for a large number of atomic pairs in the system at each time step, and further updates bond energies, bond angle energies, and many-body interaction terms based on the bond order. Its overall computational complexity is high, making it one of the main performance bottlenecks in reactive molecular dynamics simulations. This is particularly true in large-scale lithium-silicon alloy systems. Lithium-silicon alloys are key battery materials, primarily used as anodes in lithium-ion batteries and thermal batteries. Due to their high capacity and high-power discharge characteristics, they are widely used in electric vehicles, energy storage systems, and power supplies. Existing techniques often employ molecular dynamics simulations to explore the microscopic properties of lithium-silicon alloys. This allows for the analysis of lithium insertion / extraction phase transitions and lithium diffusion behavior at the atomic level, quantitative analysis of lithium insertion volume expansion and interfacial stress evolution, and clarification of SEI formation and defect regulation mechanisms, providing microscopic theoretical support for the modification and structural optimization of silicon anode materials. In existing molecular dynamics simulations of lithium-silicon alloys, many distant atomic pairs do not exhibit significant interactions, yet they are still indiscriminately included in the bond order calculation process in traditional algorithms, resulting in significant redundant computational overhead.

[0005] Existing acceleration methods typically reduce computation by using a uniform cutoff radius combined with a neighbor list. However, this method applies the same spatial cutoff range to all atom types, failing to adequately consider the differences in interaction strength and bond order decay characteristics between different atom types. Taking lithium-silicon alloy systems as an example, the interaction strength and effective range cutoff range between lithium-lithium (Li-Li), silicon-silicon (Si-Si), and lithium-silicon (Li-Si) systems differ significantly. Existing methods use a one-size-fits-all uniform cutoff, which cannot finely adjust for the different atom pair characteristics in lithium-silicon alloys. Therefore, although this method can guarantee computational accuracy, it still includes a large number of invalid atom pairs in the calculation process, limiting the overall computational efficiency improvement of complex systems such as lithium-silicon alloys.

[0006] Therefore, how to adaptively set the atomic pair cutoff range based on the interaction characteristics between different atomic types, and eliminate invalid atomic pair calculations in the simulation of lithium-silicon alloy systems while ensuring the accuracy of bond-level calculations and the physical correctness of the reaction process, and improve the computational efficiency of large-scale reaction molecular dynamics simulations, is a technical problem that urgently needs to be solved. Summary of the Invention

[0007] To address the shortcomings of existing technologies, this application provides a rapid molecular dynamics simulation method and system for lithium-silicon alloys. The core idea is to replace the traditional unified geometric truncation strategy with a type-dependent truncation mechanism based on bond-level decay characteristics, thereby significantly reducing the computation of invalid atom pairs while ensuring computational accuracy.

[0008] On the one hand, a rapid molecular dynamics simulation method for lithium-silicon alloys is provided, including the following steps: For the lithium-silicon alloy to be tested, an initial structural model of the lithium-silicon alloy is constructed based on the preset lithium-silicon ratio and lattice constant. Construct an adaptive truncation radius lookup table for atomic pairs; Under preset temperature and simulation conditions, the initial structural model of the lithium-silicon alloy is subjected to accelerated reactive molecular dynamics simulation using the atomic pair adaptive cutoff radius lookup table to perform relaxation and evolution of atomic positions and charge states. After the simulation is completed, the structural data of the lithium-silicon alloy are obtained. Based on the lithium-silicon alloy structural data, the atomic topological evolution process, bond order characteristics, or material stability of the lithium-silicon alloy under test are output.

[0009] On the other hand, a rapid molecular dynamics simulation system for lithium-silicon alloys is provided, including: The model building module is configured to: construct an initial structural model of the lithium-silicon alloy under test based on the preset lithium-silicon ratio and lattice constant; The lookup table building module is configured to: build an adaptive truncation radius lookup table for atomic pairs; The molecular dynamics simulation module is configured to: perform accelerated reactive molecular dynamics simulation on the initial structural model of the lithium-silicon alloy using the atomic pair adaptive cutoff radius lookup table under preset temperature and simulation conditions, so as to perform relaxation and evolution of atomic positions and charge states, and obtain lithium-silicon alloy structural data after the simulation is completed. The output module is configured to output the atomic topological evolution process, bond order characteristics, or material stability of the lithium-silicon alloy under test, based on the lithium-silicon alloy structure data.

[0010] The above technical solution has the following advantages or beneficial effects: In the process of rapid molecular dynamics simulation of lithium-silicon alloys, the bond-level threshold screening is moved forward to geometric pre-screening, and an atom-type-dependent adaptive cutoff radius system is constructed. This significantly reduces the proportion of invalid atom pairs entering the bond-level calculation process, and improves the efficiency of neighbor list construction and overall computation without changing the original ReaxFF physical model. Attached Figure Description

[0011] The accompanying drawings, which form part of this application, are used to provide a further understanding of this application. The illustrative embodiments of this application and their descriptions are used to explain this application and do not constitute an undue limitation of this application.

[0012] Figure 1 This is a flowchart of the method in Example 1. Detailed Implementation

[0013] It should be noted that the following detailed descriptions are exemplary and intended to provide further explanation of this application. Unless otherwise specified, all technical and scientific terms used in this invention have the same meaning as commonly understood by one of ordinary skill in the art to which this application pertains.

[0014] Example 1 This embodiment provides a rapid molecular dynamics simulation method for lithium-silicon alloys; like Figure 1 As shown, a fast molecular dynamics simulation method for lithium-silicon alloys includes the following steps: S101: For the lithium-silicon alloy to be tested, an initial structural model of the lithium-silicon alloy is constructed based on the preset lithium-silicon ratio and lattice constant. S102: Construct an adaptive cutoff radius lookup table for atom pairs; S103: Under preset temperature and simulation conditions, the initial structural model of the lithium-silicon alloy is subjected to accelerated reactive molecular dynamics simulation using the atomic pair adaptive cutoff radius lookup table to perform relaxation and evolution of atomic positions and charge states. After the simulation is completed, the structural data of the lithium-silicon alloy are obtained. S104: Based on the lithium-silicon alloy structure data, output the atomic topological evolution process, bond order characteristics, or material stability of the lithium-silicon alloy to be tested.

[0015] Further, in step S101, based on a preset lithium-silicon ratio and lattice constant, an initial structural model of the lithium-silicon alloy is constructed, including: Lithium atoms are considered as Class I atoms, and silicon atoms are considered as Class II atoms; Obtain standard unit cell structure data of lithium-silicon alloy under a preset space group, wherein the preset space group is a pre-defined theoretical space group label used to characterize the combination of three-dimensional spatial symmetry elements of the crystal; the standard unit cell structure data includes: preset space group symmetry, lattice parameters, and fractional coordinates of first-type atoms and second-type atoms; Based on the preset lithium-silicon ratio and lattice constant, the standard unit cell structure data is subjected to a three-dimensional matrix cell expansion operation to obtain a virtual simulated supercell with a target three-dimensional size, and three-dimensional periodic boundary conditions are configured. Based on the atomic occupancy state after cell expansion, a predetermined number of first-type atoms are retained or extracted and deleted within the virtual simulated supercell using a random algorithm based on a random number generator to create atomic vacancies. This ensures that the total number of first-type atoms and the total number of second-type atoms within the cell meet a preset lithium-silicon ratio requirement, resulting in an initial lithium-silicon alloy structural model. The initial lithium-silicon alloy structural model includes: three-dimensional boundary size information of the supercell, topological type labels for each atom, and three-dimensional spatial coordinate data of all retained atoms.

[0016] In this embodiment, in step S101, lithium atoms are considered as first-type atoms and silicon atoms are considered as second-type atoms.

[0017] First, a target lithium-silicon ratio is set (e.g., 3.25 lithium atoms per silicon atom), and the structural data of the standard unit cell (i.e., the smallest repeating unit in the crystal structure) of the lithium-silicon alloy under a preset space group are obtained. This ensures a stable crystalline state. Taking the alloy as an example, its preset space group (i.e., the classification label describing the symmetrical arrangement of atoms inside the crystal) is numbered 220 in the cubic crystal system. Space group. The standard unit cell structure data obtained includes: lattice parameters (i.e., the edge dimensions of the unit cell, lattice constant). (Å); and the initial relative position data of the atoms within the unit cell (containing 60 lithium atoms and 16 silicon atoms).

[0018] Next, the above standard unit cell structure data were analyzed. The three-dimensional matrix cell expansion operation (i.e., like building blocks, replicating the original cell three times in each of the length, width, and height directions to form a large box) is performed. After the cell expansion operation, a virtual simulated supercell with the target three-dimensional dimensions (length, width, and height are all 31.8 angstroms) is obtained, and three-dimensional periodic boundary conditions are configured for it (i.e., assuming that this box extends infinitely in space on all sides). At this time, the initial total number of atoms in this large box is: lithium atoms. One, silicon atom indivual.

[0019] Finally, adjustments are made based on the current number of atoms: to ensure the number of atoms in the box strictly meets the preset ratio requirement (i.e., when there are 432 silicon atoms, only 1404 lithium atoms are needed to achieve the target ratio), the difference in lithium atoms that need to be removed is calculated to be... indivual.

[0020] At this point, a random algorithm is run (i.e., the computer assigns numbers to all lithium atoms like drawing lots, and randomly selects 216 numbers). These 216 selected lithium atoms are then removed from the box, creating 216 atomic vacancies in the previously filled lattice. After the deletion is complete, the final initial structural model of the lithium-silicon alloy is obtained. This model is output as a data file, which directly contains the three-dimensional boundary dimensions of the large box, atom type labels distinguishing lithium and silicon, and the specific three-dimensional spatial coordinates of the remaining 1404 lithium atoms and 432 silicon atoms within the box.

[0021] Further, S102: constructing an adaptive truncation radius lookup table for atom pairs includes: S102-1: Analyze the reactive force field parameter file and extract the two-body parameters between different atom types; S102-2: Utilizing the monotonically decreasing property of the bond order function of the reaction force field, under a preset bond order threshold, the critical distances corresponding to each atomic type pair (including lithium-lithium, lithium-silicon, and silicon-silicon pairs) are pre-calculated using a binary search method, and the critical distances are used as the adaptive cutoff radii for the corresponding atomic type pairs. S102-3: Construct an adaptive truncation radius lookup table for atom pairs based on the adaptive truncation radius.

[0022] Furthermore, S102-1: Analyzing the reaction force field parameter file and extracting the two-body parameters between different atom types, specifically includes: During the initialization phase of the molecular dynamics simulation, the preset ReaxFF reactive force field potential function file is parsed. The file includes the monomeric property parameters of lithium atoms and silicon atoms, as well as the two-body parameters of lithium-lithium, lithium-silicon, and silicon-silicon atom pairs formed by their pairwise combinations. The two-body interaction parameter data block in the file is retrieved based on the preset atom type identifier, and the characteristic bonding distance parameters used to characterize σ bonds, π bonds, and secondary π bonds are read and extracted. , and And the corresponding empirical fitting parameters for bond-level distance decay. to ; The characteristic bonding distance parameter and the empirical fitting parameter are stored as basic physical data in memory variables.

[0023] Furthermore, the two-body parameters refer to the physical parameters used in the ReaxFF reactive force field to describe the two-body interaction and potential energy contribution between any two types of atoms in the system. In the lithium-silicon alloy system, the two-body parameters are the characteristic bonding distance parameters corresponding to lithium-lithium, lithium-silicon, and silicon-silicon atom pairs obtained by analyzing the reactive force field parameter file. , , and empirical fitting parameters for bond-level distance decay to ; When parsing the file, the two-body interaction parameter data block in the file is retrieved according to the preset atom type identifier, and the characteristic bonding distance parameters and empirical fitting parameters corresponding to the lithium-lithium, lithium-silicon, and silicon-silicon atom pairs are read and extracted. The extracted two-body parameters are directly stored in a temporary array or memory structure as basic two-body physical feature inputs. In subsequent steps S102-2, when solving for the adaptive cutoff radius of each atomic type pair through binary search, the parameters are substituted into the corresponding ReaxFF bond-level function for equation construction and calculation.

[0024] Further, S102-2: Utilizing the monotonically decreasing property of the bond order function of the reaction force field, under a preset bond order threshold, the critical distances corresponding to each atom type pair (including lithium-lithium, lithium-silicon, and silicon-silicon pairs) are pre-calculated using a binary search method, and the critical distances are used as the adaptive cutoff radii for the corresponding atom type pairs, specifically including: For any atom type Using ReaxFF key-level functions Based on the monotonically decreasing property, an equivalence relationship is established between the bond level threshold and the geometric distance; ReaxFF bond level function. The specific expression is:

[0025] in, These represent the contributions of the σ bond, π bond, and quadratic π bond correction terms, respectively.

[0026] Each component is typically represented by an exponential decay function related to the interatomic distance, for example:

[0027] in, Represents atoms and The distance between them This refers to the characteristic distance parameter corresponding to the bond type. The first empirical fitting parameters are... The second empirical fitting parameter, The third empirical fitting parameter, The fourth empirical fitting parameter, The fifth empirical fitting parameter, The sixth empirical fitting parameter is used to regulate the decay behavior of the bond order as a function of distance.

[0028] For a given key level threshold Due to the key-level function With geometric distance As it increases, it strictly monotonically decreases, and therefore there must exist a unique critical distance. Satisfy the equation ; The critical distance is calculated by iteratively solving the problem within a set distance search interval using a binary search algorithm. The specific values ​​of are such that the following relationship holds: ; in, This represents a logically necessary and sufficient condition, if and only if the real-time geometric distance between the atomic pairs is... Less than the critical distance of this type of atom pair At that time, its calculated real-time key level Only then can it exceed the set key level threshold. .

[0029] Therefore, the chemical criteria that originally required dynamic calculation of bond order during simulation for screening are equivalently transformed into a geometric criterion that relies solely on the initial geometric distance, and the calculated critical distance is used to... It is directly defined as the adaptive cutoff radius of the atom type pair, realizing the forward shift and simplification of the bond-level screening process.

[0030] For each pair of atom types, a binary search algorithm is used to search within the interval The inner solution satisfies The critical distance is used as the adaptive cutoff radius for the current atom type pair; where, This represents the set lower limit distance for the binary search, used to avoid non-physical minimum distance intervals where atoms are extremely overlapping, and to prevent overflow of the exponential term calculation; This represents the preset ReaxFF reactive force field uniform space geometric cutoff radius, which serves as the upper limit distance for the binary search.

[0031] The implementation process of the binary search method is as follows: traverse all combinations of atom types in the system in the outer loop. And construct uncorrected total bond order expressions for σ bonds, π bonds, and secondary π bonds based on the current two-body parameters of the atom pairs. ; in the initial search interval Then, the iterative loop begins, first calculating the distance to the center point of the current search interval. Subsequently, the key-level function is called to calculate the real-time total key-level value at the distance from the center point. ; Next, based on the monotonically decreasing property of the reaction force field bond order function, the interval reduction determination is performed: if the calculated total bond order value... Greater than the preset key level threshold This indicates that the current distance is still within the bonding range, and the true critical distance is located further to the right. Therefore, the search lower bound is updated to... Conversely, if Less than or equal to the preset key level threshold If the current distance exceeds the effective bonding range, and the true critical distance is located closer to the left, then the search upper limit should be updated. The loop continues to iterate through the above center point division and interval update steps until the difference between the upper and lower limits of the search interval is less than the preset convergence accuracy threshold. At this point, the upper limit of convergence will be... The critical distance for this atomic type pair is determined and assigned to the adaptive cutoff radius variable in the two-body parameter structure of the current atomic type pair.

[0032] The above process is completed during the simulation initialization phase, and the solution results are stored in the two-body parameter table, which can be directly reused in subsequent molecular dynamics time steps without repeated calculations, thus avoiding additional runtime overhead.

[0033] Further, S102-3: Based on the adaptive cutoff radius, constructing an adaptive cutoff radius lookup table for atomic pairs specifically includes: After pre-calculation, in combination of all atom types As an index, the calculated adaptive cutoff radius will be applied to each atom type. Perform structured association mapping to build an adaptive truncation radius lookup table for atomic pairs in memory; The atom pair adaptive truncation radius lookup table includes: one-dimensional and two-dimensional atom type identifier indexes, adaptive truncation radius floating-point numbers, and adaptive truncation radius squared values ​​used to accelerate distance comparisons.

[0034] Taking a system containing both lithium and silicon atoms as an example, the atom pair adaptive cutoff radius lookup table is specifically represented in memory as a symmetric matrix structure with atom type numbers as row and column indices. An example of its specific storage information is as follows: When traversing to the index combination When it is a lithium-lithium (Li-Li) atom pair, this entry stores its proprietary adaptive cutoff radius (e.g., 5.0 angstroms) and square value (25.0 square angstroms). When traversing to the index combination or When the lithium-silicon (Li-Si) atom pair is represented, this entry stores the adaptive cutoff radius (e.g., 4.2 Å) and square value (17.64 square Å) corresponding to the lithium-silicon heterobond. When traversing to the index combination When it is a silicon-silicon (Si-Si) atom pair, this entry stores the adaptive cutoff radius (e.g., 4.8 angstroms) and square value (23.04 square angstroms) of the silicon-silicon atom pair.

[0035] By using this matrix-form adaptive cutoff radius lookup table for atom pairs, subsequent steps can directly read the target cutoff radius based on the current combination of atom types, avoiding repeated bond-level function solving during molecular dynamics simulations.

[0036] Further, in step S103: under preset temperature and simulated conditions: The preset temperature is typically between 300 Kelvin and 800 Kelvin. The simulation conditions include temperature and pressure control using an NVT ensemble (i.e., a constant number of atoms, constant volume, and constant temperature ensemble) or an NPT ensemble (i.e., a constant number of atoms, constant pressure, and constant temperature ensemble), and the time integration step of the kinetic iteration is set to 0.1 femtoseconds to 0.5 femtoseconds.

[0037] Further, S103: using the atom pair adaptive cutoff radius lookup table to perform accelerated reactive molecular dynamics simulations on the initial structural model of the lithium-silicon alloy, specifically refers to: Using the aforementioned adaptive cutoff radius lookup table for atomic pairs, reactive molecular dynamics simulations were performed on the initial structural model of the lithium-silicon alloy in the molecular dynamics software LAMMPS to conduct relaxation and evolution of atomic positions and charge states. After the simulation was completed, the structural data of the lithium-silicon alloy were obtained.

[0038] Further, in S103: the initial structural model of the lithium-silicon alloy is subjected to accelerated reactive molecular dynamics simulation using the atom pair adaptive cutoff radius lookup table, and each time step in the simulation process includes the following sub-steps: S103-1: Based on the combination of atomic types of candidate neighbor atomic pairs, select the corresponding adaptive truncation radius from the atomic pair adaptive truncation radius lookup table, perform initial spatial distance screening on candidate neighbor atomic pairs based on the adaptive truncation radius, retain only atomic pairs that meet the truncation conditions, and store the atomic pairs that meet the truncation conditions in the far neighbor list to obtain the far neighbor list. S103-2: In the chemical bond construction stage, the candidate atom pairs in the distant neighbor list are screened a second time using an adaptive cutoff radius to obtain the effective atom pairs after screening, so as to eliminate the invalid atom pairs that do not have effective bond level contributions. S103-3: Call the bond level calculation function to determine the bond level of the effective atom pairs, and calculate the bond energy, many-body terms and non-bonded interactions based on the determination results. Update the atom states iteratively through time integration until the simulation is completed.

[0039] Further, S103-1: Selecting the corresponding adaptive cutoff radius from the adaptive cutoff radius lookup table of the atomic pairs according to the atomic type combination of the candidate neighbor atomic pairs, specifically includes: At each time step of the simulation, atomic pairs in the LAMMPS native nearest neighbor list are first obtained as candidate neighbor atomic pairs. Subsequently, the type combination of the two atoms in the current candidate neighbor atom pairs is dynamically extracted. In the lithium-silicon alloy system, this specifically includes lithium-lithium atom pairs, silicon-silicon atom pairs, and lithium-silicon atom pairs. Next, using the type number of the two atoms as the second-order row and column index value, the search is performed in the adaptive cutoff radius lookup table of atomic pairs. When the type combination of the current candidate neighbor atomic pairs is a lithium-lithium atomic pair, the corresponding lithium-lithium adaptive cutoff radius in the lookup table is matched and selected. When the type combination of the current candidate neighbor atom pairs is silicon-silicon atom pairs, match and select the corresponding silicon-silicon adaptive cutoff radius in the lookup table; When the type combination of the current candidate neighbor atom pairs is a lithium-silicon atom pair, match and select the corresponding lithium-silicon adaptive cutoff radius in the lookup table; The extracted adaptive cutoff radius is used to perform initial spatial distance screening of the current candidate neighbor atom pairs.

[0040] Further, S103-1: Initially screen candidate neighbor atom pairs by spatial distance using an adaptive cutoff radius, retaining only atom pairs that meet the cutoff condition, and storing these atom pairs in a distant neighbor list to obtain the distant neighbor list, specifically including: During the process of traversing candidate neighbor atom pairs, the current atom is calculated in real time. With atoms The three-dimensional Euclidean distance between them ; Then, the spatial distance The adaptive cutoff radius extracted from the lookup table corresponding to the current combination of atom types. Perform a geometric size comparison, wherein the truncation condition refers to satisfying a mathematical inequality. ; If the distance or squared distance of the current atomic pair satisfies the truncation condition, the current atomic pair is determined to have potential chemical bonding or effective bond level contribution in physics. The current atomic pair is retained as a valid far neighbor, and its global atomic index number and related topological information are written into the far neighbor list data structure used internally by ReaxFF. Otherwise, if the truncation condition is not met, the current atomic pair is determined to be an invalid atomic pair, and it is directly ignored and removed from the list.

[0041] This spatial distance screening process eliminates a large number of redundant atom pairs at long distances before proceeding to complex bond-level energy calculations, thus obtaining a simplified list of distant neighbors.

[0042] Convert the native LAMMPS neighbor list into a distant neighbor list used internally by ReaxFF.

[0043] The LAMMPS native neighbor list refers to a universal nearest neighbor index table built based on a universal nonbonded geometric cutoff radius. It uses the local index of the current local atom as the row axis and stores the global numbers of all candidate neighbor atoms within that spatial range without discrimination.

[0044] The distant neighbor list used internally by ReaxFF is a simplified version of the reaction nearest neighbor list used for subsequent bond-level dynamic determination. In memory, it is represented as a compact array of structures containing the characteristic parameters of the atom pair interactions.

[0045] The LAMMPS native neighbor list is converted into a far neighbor list used internally by ReaxFF. The specific conversion process is as follows: at the beginning of each time step, the memory space of the far neighbor list is first initialized. Then, a double nested loop is used. The outer loop iterates through each local atom in the LAMMPS native neighbor list one by one, and the inner loop reads all the candidate neighbor atoms corresponding to the atom in turn. During the traversal, the spatial distance screening process of step S103-1 is called in real time, and the adaptive cutoff radius corresponding to the current atom type combination is used to determine whether the current atom pair meets the cutoff condition. If the conditions are met, the local number of the neighboring atom, the two-body parameter index, and the real-time distance information are pushed into and stored in the far neighbor list structure, and the far neighbor counter of the atom is incremented by one. Conversely, if the truncation condition is not met, the atom is skipped directly until the traversal is completed and the transformation is finished.

[0046] In traditional implementations, this stage typically uses a uniform cutoff radius to filter neighboring atoms. However, this invention dynamically selects a type-related cutoff radius based on the atom type to pre-screen candidate neighboring atoms, retaining only atom pairs that meet the cutoff criteria in the distant neighbor list, thereby reducing the number of atom pairs that need to be processed in subsequent bond-level calculations.

[0047] Furthermore, for local atoms (atoms located within a region that are actually managed and have their positions and force states updated by the current computing core in the domain decomposition parallel computation) and ghost atoms (boundary atoms managed by spatially adjacent computing regions, but whose position information is copied and communicated to the current computing core's memory for calculating cross-boundary multi-body interactions), different adaptive cutoff radii can be adopted based on their atom type combinations to achieve a more refined neighbor filtering mechanism and reduce the storage and traversal overhead of invalid neighbor atom pairs. For local atoms, a non-bonded cutoff radius (the maximum spatial geometric cutoff distance, uniformly specified by the global control parameters of the reaction force field, including long-range van der Waals and Coulomb interactions) is used; for ghost atoms, a type-related cutoff radius (the adaptive cutoff radius constructed in step S102-3, set separately for lithium-lithium, lithium-silicon, and silicon-silicon atom pairs) is used.

[0048] In the specific implementation of the judgment, the specific cutoff condition formula is as follows: ; in, This represents the real-time three-dimensional geometric distance between the current local atom and the ghost atom (which is specifically a lithium atom or a silicon atom). This indicates the type of atom pair currently interacting. The corresponding adaptive cutoff radius value. Therefore, the neighbor list construction process is based on the unified cutoff radius condition of the traditional ReaxFF force field: ; in, Similarly, it represents the geometric distance between two atoms. This indicates the uniform spatial geometric cutoff radius used by the traditional algorithm to cut all atom pairs in one stroke; Extended to the type-related truncation condition proposed in this invention: in, This represents the real-time geometric distance between the current atom pairs. This indicates that the adaptive cutoff radius is dynamically retrieved from the lookup table based on the current atom combination (such as lithium-lithium, silicon-silicon, or lithium-silicon combination), thereby enabling differentiated spatial filtering of different atom type pairs during the boundary communication and list construction phases of parallel computing.

[0049] A multi-layered screening mechanism is introduced during the neighbor list construction stage, including initial screening based on spatial distance and type-related truncation screening, to retain only atomic pairs that are physically likely to form effective bond-level contributions for subsequent ReaxFF energy calculations.

[0050] Further, in S103-2: during the chemical bond construction stage, an adaptive cutoff radius is used to perform a secondary screening of candidate atom pairs in the distant neighbor list to obtain the screened effective atom pairs, thereby eliminating invalid atom pairs that do not contribute to the effective bond level. Specifically, this includes: In the chemical bond building branch, each candidate atom pair in the distant neighbor list is read sequentially to obtain the real-time spatial distance between the two atoms. The corresponding adaptive cutoff radius is retrieved from the lookup table based on the combination of its atom types (such as lithium-lithium, lithium-silicon, or silicon-silicon). ; Subsequently, the real-time spatial distance is determined. Does it meet the secondary screening criteria? ; If the conditions are met, the candidate atom pair is determined to be a valid atom pair and remains in an active state. The table entry data of the valid atom pair includes the local topological index and real-time distance of the lithium atom or silicon atom. Conversely, if If the distance between the candidate atoms is too great, their electron orbitals will not overlap effectively and their total bond level contribution will be close to 0. Therefore, they are considered invalid atom pairs with no effective bond level contribution.

[0051] In the specific elimination process, the invalid atom pairs are skipped directly during the loop traversal, and they are refused to be substituted into the subsequent complex bond-level determination and exponent term update process, thereby eliminating redundant calculations at the very beginning.

[0052] Furthermore, in S103-2: during the chemical bond construction stage, after performing a secondary screening of candidate atom pairs in the distant neighbor list using an adaptive cutoff radius to obtain the screened effective atom pairs, the process further includes: ReaxFF constructs a list of chemical bonds based on geometric distance and bond order criteria: After eliminating invalid atom pairs using the adaptive cutoff radius, the real-time total bond level is calculated for the remaining valid atom pairs (such as lithium-lithium, lithium-silicon, or silicon-silicon atom pairs that satisfy distance constraints). ;

[0053] Subsequently, the real-time total key level value will be... The value is compared with a preset chemical bond determination threshold, and the result is determined if and only if the real-time total bond order value is [value missing]. If the value is greater than or equal to this threshold, it is determined that the two atoms (such as the lithium atom and the silicon atom) have formed a valid chemical bond physically. For example, the preset chemical bond determination threshold is a pre-set numerical pass line used to filter out invalid chemical bonds, such as a fixed empirical value of 0.1; if and only if the real-time total bond level value When the value is greater than or equal to this threshold, it is determined that the two atoms (such as the lithium atom and the silicon atom) have formed a valid chemical bond physically. For example, if the real-time total bond order value between a specific lithium atom and a silicon atom at a certain instant... The calculated result is 0.45. Since it is greater than or equal to the preset threshold of 0.1, it is determined that the lithium atom and silicon atom have formed a valid chemical bond in the current state, and this is added to the chemical bond list; otherwise, if the real-time total bond level value is lower... If the value drops to 0.08, since it is less than the preset threshold of 0.1, it is determined that no effective chemical bond has been formed between the two atoms or that the original chemical bond has been broken.

[0054] Next, the local numbers of the two bonded atoms, the real-time bond order components, and the topological connections are used as bond records and written sequentially and appended to the chemical bond list in memory, thus completing the dynamic construction of the chemical bond list.

[0055] The purpose and significance of constructing the chemical bond list is to capture the dynamic behavior of bond breaking and bonding in lithium-silicon alloy systems during charging, discharging or alloying processes at the microscopic level, accurately describe the structural evolution and phase transition of materials without pre-defining the topological structure, abstract complex chemical structures into discrete bond connections, and provide core topological basis for subsequent energy and force calculations.

[0056] The constructed list of chemical bonds is directly iterated and called in subsequent calculation steps of the molecular dynamics simulation.

[0057] In traditional implementations, all atomic pairs that enter the distant neighbor list need to participate in the subsequent bond-level judgment process, resulting in a large number of invalid atomic pairs entering the bond-level calculation stage.

[0058] To address this issue, this invention continues to employ a type-related truncation mechanism in the bond-level screening and chemical bond construction stages. This mechanism performs a secondary screening of candidate atom pairs in the distant neighbor list, retaining only atom pairs that may contribute to the effective bond level for subsequent chemical bond determination and search processes. This further reduces the overhead of invalid bond level calculations and traversal. Subsequently, for atom pairs that meet the truncation criteria, the bond-level calculation function is invoked to determine the bond level and construct the corresponding chemical bond list.

[0059] Further, S103-3: The effective atom pairs are subjected to bond-level determination using a bond-level calculation function, and bond energies, many-body terms, and non-bonded interactions are calculated based on the determination results. The atom states are iteratively updated through time integration until the simulation is complete. Specifically, this includes: S103-31: In the energy calculation stage at each time step, extract the effective atom pairs (such as lithium-lithium, lithium-silicon, or silicon-silicon atom pairs that satisfy distance constraints) retained after secondary screening, and call the bond level calculation function to calculate the real-time total bond level value. ; S103-32: Obtain the current real-time total bond order value for each effective atom pair. Then, dynamic physical corrections are made based on the local overcoordination of atoms, and the results are compared with a preset bonding cutoff threshold to determine whether the effective atom pair meets the minimum chemical bond interaction requirements.

[0060] Further, step S103-32 specifically includes the following steps: S103-321: First, based on the bond order of all atoms connected to the target atom in the current simulation step, the real-time total valence bond number of the target atom is accumulated; then, the real-time total valence bond number is subtracted from the inherent ideal valence of the element to which the atom belongs to calculate the overcoordination factor; when the overcoordination factor is greater than zero, the current atom is determined to be in a locally overcoordinating state; finally, using an exponential decay function based on the overcoordination factor, the real-time total bond order value of all effective atom pairs involved by the current atom is dynamically reduced and weakened to obtain a corrected real-time total bond order value, so as to eliminate the virtual high bond order caused by atomic compression and restore the real physical interaction; S103-322: The corrected real-time total bond level value is compared with a preset bond truncation threshold. If the corrected real-time total bond level value is greater than or equal to the bond truncation threshold, it is determined that there is a minimum chemical bond interaction between the two atoms, the current atom pair is retained and output as valid bonding data; if the corrected real-time total bond level value is less than the bond truncation threshold, it is determined that the interaction between the two atoms is extremely weak, the bond level value of the current atom pair is directly cleared to zero and discarded.

[0061] S103-33: Based on the results of the bond level determination, further parallel calculations of bond energy, many-body terms, and non-bonded interactions are performed; S103-331: Bond energy calculation process: Directly using the corrected real-time total bond order value as input, substitute it into the two-body bond energy decay formula that depends on the atom type pair to calculate the total bond energy directly contributed by the bonding behavior. This causes energy to shift and evolve as chemical bonds break or form; Further, S103-331 specifically includes the following steps: S103-3311: Calculation of energy components of multi-physics single bonds: The corrected real-time total bond level value is decomposed into contribution components based on different electron orbital overlap types; The components specifically include: the corrected Key level Corrected Key level And the revised Key level ; By combining the intrinsic bond energy parameters of specific atom type pairs, the energies of different single bond types are calculated using the two-body bond energy decay formula. The specific mathematical expression of the two-body bond energy decay formula is as follows:

[0062] in, This represents the total bond energy directly contributed by the bonding behavior; , , These represent the target atoms within the current simulation step. With atoms Corrected between key, Key and The real-time key-level components of the key, and the sum of the three is equal to the corrected real-time total key-level value; Indicates a specific type of atom in the formation of a standard Intrinsic dissociation energy parameters for single bonds; and Each represents a specific pair of atoms in the formation Key and Intrinsic dissociation energy parameters of the bond; and This represents the empirical energy decay correction index, which depends on a specific pair of atom types. S103-3312: Continuous slip evolution determination: Due to the total bond energy in the two-body bond energy decay formula... It is a continuously differentiable function of each real-time bond level component. The total bond energy is calculated as the two atoms move further apart during thermal motion, causing each real-time bond level component to gradually decay and approach zero. It also smoothly slides and transitions to zero, thereby physically eliminating the energy jump caused by forced truncation at the cutoff radius, and realizing the dynamic and smooth evolution of energy as chemical bonds break or form.

[0063] S103-332: Multibody Calculation Process: Retrieve the chemical bond list, search for all three-atom combinations (such as lithium-silicon-lithium) and consecutive four-atom sequences sharing a central atom in the current initial structural model of the lithium-silicon alloy, and substitute them into the specific three-body valence angle energies of the reaction force field. Equations and torsional energy of a four-body dihedral angle In the equation, the multibody energy terms and the forces on the corresponding atoms are calculated by using the product correction factor related to the bond order. Furthermore, S103-332 specifically includes: for any retrieved atom The triatomic combination consists of atoms To share the central atom, the reaction force field used to substitute it has a specific three-body valence angle energy. The specific mathematical expression of the equation is:

[0064] in, This represents the valence angle energy component generated due to bond angle bending; Indicates the atoms in the current simulation step Real-time physical bond angles between them; Indicates the central atom The ideal equilibrium bond angle is determined by the type of atoms connected to it; and This represents the intrinsic energy constant and angular decay exponent parameters that depend on a specific combination of triatomic types. and These represent product correction factors related to bond order, used to dynamically adjust the contribution weight of bond angle energy based on the current bond saturation. Their calculation formula is as follows: , in, and These are empirical parameters. and These are the corrected real-time total bond order values ​​for the corresponding atomic pairs; For any one of the retrieved atoms The tetrahedral torsional energy of a continuously arranged tetraatomic sequence is... The specific mathematical expression of the equation is:

[0065] in, This represents the dihedral torsional energy component caused by the key shaft torsion. Indicates the number of atoms in the current simulation step. The plane formed by atoms The real-time physical dihedral angular displacement between the planes formed; The intrinsic dihedral torsional barrier parameter corresponding to a specific four-atom sequence; The product correction factor associated with consecutive triple bond orders is calculated as follows:

[0066] in For a specific torsional correction index, , and These are the corrected real-time total bond order values ​​for three adjacent chemical bonds in the four-atom sequence.

[0067] S103-333: Calculation process for non-bonded interactions: For all atomic pairs in the system, their chemical bond properties are stripped, and their Euclidean distances are substituted into the van der Waals potential equation shielded from Coulomb potential and step function constraints. The long-range Coulomb interaction energies of the entire system are then solved. With non-bonded van der Waals interaction energy Specifically: For any two atoms in the entire system and its Euclidean distance Substituting into the shielded Coulomb potential equation, solve for the long-range Coulomb interaction energy of the entire system. The specific mathematical expression of the shielded Coulomb potential equation is as follows:

[0068] in, Represents the long-range Coulomb interaction energy component; This is the constant term obtained by multiplying the dielectric constant by the unit conversion factor; and Each represents an atom and atoms Real-time dynamic charge value at the current simulation step; Indicates dependence on atoms With atoms The shielding parameter for the element type is used to eliminate electrostatic potential divergence and jumps between two atoms at extremely close distances; At the same time, the Euclidean distance Substituting the van der Waals potential equation with step function constraints, solve for the unbonded van der Waals interaction energy of the entire system. The specific mathematical expression of the van der Waals potential equation constrained by the step function is as follows:

[0069] in, This represents the energy component of the nonbonded van der Waals interaction; The van der Waals potential depth parameter representing a specific pair of atom types; This represents the slope constant of the potential energy curve, which depends on a specific pair of atom types. The van der Waals equilibrium distance parameter represents a specific pair of atom types; The step distance function with adaptive shielding correction is introduced, and its specific mathematical expression is as follows: , in, and This represents the preset nonbonding distance correction index and shielding constant.

[0070] After superimposing the above energies and differentiating them to obtain the resultant force on each lithium atom and silicon atom in the lithium-silicon alloy, the atomic state is updated iteratively through time integration, specifically: Within each simulation step, the total bond energy is... Three-body valence angle energy Tetrahedral torsional energy Long-range Coulomb interaction energy With non-bonded van der Waals interaction energy By superimposing potential energy, the total potential energy of the system can be obtained. ; Regarding the total potential energy By taking the negative gradient derivative with respect to the three-dimensional spatial coordinates of each atom, we obtain the resultant force vector of each lithium atom and silicon atom in the system at the current simulation step. : ; Next, Newton's laws of motion are retrieved, and the resultant force vector is... The data is substituted into the time integrator and iteratively solved using the Velocity Verlet time integration algorithm to calculate and update the real-time velocity and three-dimensional spatial coordinate state of each lithium atom and silicon atom in the next simulation step.

[0071] The aforementioned time integration refers to using a numerical integration algorithm to solve Newton's equations of motion within a set discrete time step (such as 0.1 femtoseconds) to convert the net force on the atom at the current moment into acceleration, and then iteratively calculating the new velocity and new three-dimensional spatial coordinate position of each atom at the next moment. At the same time, the algorithm dynamically adjusts the atomic velocity to maintain the set thermodynamic temperature.

[0072] S103-34: Repeatedly execute the discrimination, solution and integration process in each time step of S103-31~S103-33 until the preset total number of simulation time steps is reached, and the simulation is completed and the lithium silicon alloy structure data is output.

[0073] Further, step S104: Based on the lithium-silicon alloy structural data, output the atomic topological evolution process, bond order characteristics, or material stability of the lithium-silicon alloy under test, specifically including: S104-1: Output of the atomic topological evolution process: Extract the text trajectory file output by LAMMPS, which contains the real-time three-dimensional coordinates of lithium atoms and silicon atoms at each time step, import it into three-dimensional visualization software for continuous bonding rendering, and thus output the microscopic topological evolution process of lithium-silicon alloy. S104-2: Output of bond order characteristics: Obtain the bond order text file calculated by ReaxFF in LAMMPS, and read the real-time total bond order of lithium-lithium, lithium-silicon, or silicon-silicon atomic pairs in it to characterize the strength evolution of chemical bonds in lithium-silicon alloys. S104-3: Output of material stability: Extract the total energy, pressure and volume time series data of the system recorded in the LAMMPS thermodynamic output file, evaluate the thermodynamic equilibrium by analyzing energy fluctuations, directly calculate the volume expansion rate of the material using volume changes, and output the lithium ion diffusion coefficient by calculating the root mean square displacement through the extracted coordinates, and finally quantitatively output the material stability.

[0074] The relative energy fluctuation ratio is obtained by calculating the ratio of the time-series root mean square fluctuation of the total energy of the system to the average value of the total energy. When the relative fluctuation ratio is less than or equal to a preset equilibrium threshold, the system is determined to have reached thermodynamic equilibrium. Calculating the volumetric expansion rate of a material using volume change The calculation formula is as follows:

[0075] in, For the initial volume, The average volume after thermodynamic equilibrium; Using the extracted real-time coordinate vectors of lithium atoms Coordinate vector at the initial equilibrium time Calculate the root mean square displacement :

[0076] in, This represents the total number of lithium atoms. Finally, the diffusion coefficient of lithium ions is solved and output using Einstein's relation. : .

[0077] After applying the above optimization scheme, the single-core performance of initial structure models of lithium-silicon alloys of different scales was tested on the experimental platform using LAMMPS software. The experimental data are shown in Tables 1, 2, and 3. The test results show that the reaction force field acceleration method based on adaptive truncation of atom type proposed in this invention can effectively reduce the number of invalid atom pairs calculated in the lithium-silicon alloy system, thereby improving the overall execution efficiency of the ReaxFF force field.

[0078] Specifically, in the distant neighbor list construction and screening stage, approximately 1.08 times performance improvement was achieved due to the screening of redundant neighbor atom pairs; in the distant bond level screening and chemical bond construction stage, approximately 1.2-1.3 times performance improvement was achieved due to the reduction in the computational scale of candidate bond levels. Overall, the results show that this invention can achieve an overall performance improvement of approximately 5%-10% while ensuring simulation accuracy and system physical correctness.

[0079] Furthermore, this invention does not require modification of the original energy expression and physical model of ReaxFF, but only optimizes the neighbor screening and candidate atom pair generation mechanism. Therefore, it has good engineering compatibility and portability. This method can be used as a prerequisite for further optimization and can be directly integrated into the existing ReaxFF molecular dynamics simulation framework.

[0080] Table 1. Performance Comparison Before and After Optimization in the Construction and Screening Stages of the Far Neighbor List in the Lithium-Silicon Alloy Model

[0081] Table 2. Performance Comparison Before and After Optimization in Bond Order Screening and Chemical Bond Building Stages of Lithium-Silicon Alloy Model

[0082] Table 3. Performance Comparison Before and After Using the Optimization Strategy of This Invention in the Lithium-Silicon Alloy Model

[0083] The core feature of the Reactive Force Field (ReaxFF) lies in its use of a continuously changing bond order mechanism to describe the bonding and breaking processes between atoms. Unlike traditional fixed-bond models, ReaxFF dynamically calculates the bond order through interatomic distances and constructs the total energy of the system based on the bond order, including various energy terms such as bond energy, bond angle energy, torsional energy, and non-bonded interactions, thereby achieving a continuous description of the entire chemical reaction process.

[0084] In ReaxFF, the total energy of the system It is usually expressed as a superposition of the following multiple energy contributions: ; in, For bond energy, For lone pair electron energy, and These are overcoordination and undercoordination correction terms, respectively. For valence angle energy, This is the dihedral torsional energy, and and These represent the non-bonded Coulomb interaction energy and the van der Waals interaction energy, respectively. In this energy system, all reactive many-body energy terms depend directly or indirectly on the dynamic bond order between atoms. The most fundamental bond energy term... For example, its magnitude is entirely determined by the instantaneous bond order between atoms, thus ensuring a smooth evolution of energy during bond breaking and formation. Therefore, accurately calculating the dynamic bond order between atoms is a key foundation for solving the energy and forces of the ReaxFF system.

[0085] In ReaxFF, atoms With atoms The overall key level is usually defined as the superposition of multiple contributors:

[0086] in These represent the contributions of the σ bond, π bond, and quadratic π bond correction terms, respectively.

[0087] Each component is typically represented using an exponential decay form related to the interatomic distance, for example:

[0088] in, Represents atoms and The distance between them This refers to the characteristic distance parameter corresponding to the bond type. The first empirical fitting parameters are... The second empirical fitting parameter, The third empirical fitting parameter, The fourth empirical fitting parameter, The fifth empirical fitting parameter, The sixth empirical fitting parameter is used to regulate the decay behavior of the bond order as a function of distance.

[0089] Through the above form, the key level It can vary continuously with the interatomic distance, so that the bonding and breaking processes can be described in the form of a smooth function, avoiding the discontinuity problem caused by topological abrupt changes in traditional force fields.

[0090] Existing ReaxFF reactive force fields typically employ a uniform cutoff radius strategy during neighbor atom screening, meaning that the same spatial cutoff range is used to construct the neighbor list for all atom types. However, significant differences exist in the interaction strength, bond order decay characteristics, and effective bonding distances among different atom types, leading to the following problems with the uniform cutoff strategy: (1) Using the same cutoff range for different atom types lacks specificity.

[0091] Existing implementations typically use a uniform geometric cutoff radius to construct the bond-order neighbor list. Since the effective bonding distances of different atom pairs (such as carbon-carbon, carbon-hydrogen, and hydrogen-hydrogen) vary greatly, the uniform cutoff radius needs to cover the atom pairs with the furthest interaction range. This results in a large number of atom pairs that cannot form effective chemical bonds physically or have extremely low bond orders entering the bond-order calculation process, thus generating a large amount of unnecessary exponential operations and memory access overhead.

[0092] (2) A large number of invalid atom pairs are included in the bond-level calculation process.

[0093] Under the uniform truncation strategy, a large number of atomic pairs that are far apart and do not contribute effectively to bond order are still retained in the neighbor list. These invalid atomic pairs still require complex exponential function calculations and conditional judgments in subsequent bond order calculations, resulting in significant redundant computational overhead. This problem is particularly prominent in large-scale system simulations.

[0094] (3) The uniform cutoff radius limits the improvement of ReaxFF computation efficiency.

[0095] Since ReaxFF requires calculating the bond level and related action terms for each atomic pair based on the neighbor list, the uniform truncation strategy will result in a large neighbor list, increasing the computational load during atomic pair traversal and bond level update, thus limiting further improvement in the overall computational efficiency of ReaxFF.

[0096] Therefore, the technical problem to be solved by this invention is: how to adaptively set the atom pair cutoff range according to the interaction characteristics and bond order decay law between different atom types, so as to reduce the number of invalid atom pairs entering the neighbor list and participating in the bond order calculation while ensuring the accuracy of ReaxFF bond order calculation and the physical correctness of the reaction process, thereby reducing the redundant overhead in the ReaxFF calculation process and improving the computational efficiency of reaction molecular dynamics simulation.

[0097] In the original ReaxFF implementation, a list of distant neighbors is first constructed using a fixed geometric cutoff radius. Then, the bond order is calculated for each atomic pair in the list, and a chemical bond list is constructed by filtering atomic pairs that meet certain criteria using a fixed bond order threshold. However, this process inherently relies on a uniform spatial cutoff and a uniform bond order threshold, failing to consider the differences in bond order decay behavior between different atom types. This results in a large number of non-bonding atomic pairs still entering the bond order calculation stage, generating redundant computational overhead.

[0098] Example 2 This embodiment provides a rapid molecular dynamics simulation system for lithium-silicon alloys, including: The model building module is configured to: construct an initial structural model of the lithium-silicon alloy under test based on the preset lithium-silicon ratio and lattice constant; The lookup table building module is configured to: build an adaptive truncation radius lookup table for atomic pairs; The molecular dynamics simulation module is configured to: perform accelerated reactive molecular dynamics simulation on the initial structural model of the lithium-silicon alloy using the atomic pair adaptive cutoff radius lookup table under preset temperature and simulation conditions, so as to perform relaxation and evolution of atomic positions and charge states, and obtain lithium-silicon alloy structural data after the simulation is completed. The output module is configured to output the atomic topological evolution process, bond order characteristics, or material stability of the lithium-silicon alloy under test, based on the lithium-silicon alloy structure data.

[0099] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A method for fast simulation of molecular dynamics of lithium-silicon alloys, characterized by, Includes the following steps: For the lithium-silicon alloy to be tested, an initial structural model of the lithium-silicon alloy is constructed based on the preset lithium-silicon ratio and lattice constant. Construct an adaptive truncation radius lookup table for atom pairs, including: Analyze the reactive force field parameter file and extract the two-body parameters between different atom types; By utilizing the monotonically decreasing property of the bond order function of the reaction force field, under a preset bond order threshold, the critical distance corresponding to each atomic type pair is pre-calculated using a binary search method, and the critical distance is used as the adaptive cutoff radius of the corresponding atomic type pair. Based on the adaptive cutoff radius, an adaptive cutoff radius lookup table for atomic pairs is constructed; Utilizing the monotonically decreasing property of the bond order function of the reaction force field, under a preset bond order threshold, the critical distance corresponding to each atomic type pair is pre-calculated using a binary search method, and the critical distance is used as the adaptive cutoff radius for the corresponding atomic type pair. Specifically, this includes: For any atomic type pair , the equivalent relationship between the bond order threshold and the geometric distance is established by using the monotonic decreasing property of the ReaxFF bond order function ; the specific expression of the ReaxFF bond order function is: ； wherein respectively represent the contribution of the σ, π and quadratic π bond correction terms; Each component is typically represented using an exponential decay function related to the interatomic distance: ； ； ； wherein represents the distance between atoms , , is a characteristic distance parameter for the corresponding bond type, is a first empirical fitting parameter, is a second empirical fitting parameter, is a third empirical fitting parameter, is a fourth empirical fitting parameter, is a fifth empirical fitting parameter, is a sixth empirical fitting parameter, the empirical fitting parameters being used to tune the decay behavior of the bond order with distance; For a given key level threshold , the key level function is strictly monotonically decreasing with increasing geometric distance , so there is a unique critical distance satisfying the equation ; The critical distance is calculated by iteratively solving the problem within a set distance search interval using a binary search algorithm. The specific values ​​of are such that the following relationship holds: ； in, This represents a logically necessary and sufficient condition, if and only if the real-time geometric distance between the atomic pairs is... Less than the critical distance of this type of atom pair At that time, its calculated real-time key level Only then can it exceed the set key level threshold. ; For each pair of atom types, a binary search algorithm is used to search within the interval The inner solution satisfies The critical distance is used as the adaptive cutoff radius for the current atom type pair; where, This represents the set lower limit distance for the binary search, used to avoid non-physical minimum distance intervals where atoms are extremely overlapping, and to prevent overflow of the exponential term calculation; This represents the predefined geometric cutoff radius of the unified space of the ReaxFF reactive force field, which serves as the upper limit distance for the binary search. Under preset temperature and simulation conditions, the initial structural model of the lithium-silicon alloy is subjected to accelerated reactive molecular dynamics simulation using the atomic pair adaptive cutoff radius lookup table to perform relaxation and evolution of atomic positions and charge states. After the simulation is completed, the structural data of the lithium-silicon alloy are obtained. Based on the lithium-silicon alloy structural data, the atomic topological evolution process, bond order characteristics, or material stability of the lithium-silicon alloy under test are output.

2. The rapid molecular dynamics simulation method for lithium-silicon alloys as described in claim 1, characterized in that, Based on the preset lithium-silicon ratio and lattice constant, an initial structural model of the lithium-silicon alloy is constructed, including: Lithium atoms are considered as Class I atoms, and silicon atoms are considered as Class II atoms; Obtain standard unit cell structure data of lithium-silicon alloy under a preset space group, wherein the preset space group is a pre-defined theoretical space group label used to characterize the combination of three-dimensional spatial symmetry elements of the crystal; the standard unit cell structure data includes: preset space group symmetry, lattice parameters, and fractional coordinates of first-type atoms and second-type atoms; Based on the preset lithium-silicon ratio and lattice constant, the standard unit cell structure data is subjected to a three-dimensional matrix cell expansion operation to obtain a virtual simulated supercell with a target three-dimensional size, and three-dimensional periodic boundary conditions are configured. Based on the atomic occupancy state after cell expansion, a predetermined number of first-type atoms are retained or extracted and deleted within the virtual simulated supercell using a random algorithm based on a random number generator to create atomic vacancies. This ensures that the total number of first-type atoms and the total number of second-type atoms within the cell meet a preset lithium-silicon ratio requirement, resulting in an initial lithium-silicon alloy structural model. The initial lithium-silicon alloy structural model includes: three-dimensional boundary size information of the supercell, topological type labels for each atom, and three-dimensional spatial coordinate data of all retained atoms.

3. The rapid molecular dynamics simulation method for lithium-silicon alloys as described in claim 1, characterized in that, Accelerated reactive molecular dynamics simulations were performed on the initial structural model of the lithium-silicon alloy using the aforementioned atom pair adaptive cutoff radius lookup table. Each time step in the simulation process included the following sub-steps: Based on the combination of atomic types of candidate neighbor atomic pairs, the corresponding adaptive truncation radius is selected from the adaptive truncation radius lookup table of atomic pairs. Based on the adaptive truncation radius, the candidate neighbor atomic pairs are initially screened by spatial distance. Only atomic pairs that meet the truncation conditions are retained. The atomic pairs that meet the truncation conditions are stored in the far neighbor list to obtain the far neighbor list. In the chemical bond construction stage, the candidate atom pairs in the distant neighbor list are screened a second time using an adaptive cutoff radius to obtain the effective atom pairs after screening, so as to eliminate the invalid atom pairs that do not have effective bond level contributions. The effective atom pairs are used to call the bond level calculation function to determine the bond level, and the bond energy, many-body terms and non-bonded interactions are calculated based on the determination results. The atom states are updated iteratively through time integration until the simulation is completed.

4. The rapid molecular dynamics simulation method for lithium-silicon alloys as described in claim 3, characterized in that, Based on the combination of atom types of candidate neighbor atom pairs, the corresponding adaptive cutoff radius is selected from the adaptive cutoff radius lookup table for the atom pairs, specifically including: At each time step of the simulation, atomic pairs in the LAMMPS native nearest neighbor list are first obtained as candidate neighbor atomic pairs. Subsequently, the type combination of the two atoms in the current candidate neighbor atom pairs is dynamically extracted. In the lithium-silicon alloy system, this specifically includes lithium-lithium atom pairs, silicon-silicon atom pairs, and lithium-silicon atom pairs. Next, using the type number of the two atoms as the second-order row and column index value, the search is performed in the adaptive cutoff radius lookup table of atomic pairs. When the type combination of the current candidate neighbor atomic pairs is a lithium-lithium atomic pair, the corresponding lithium-lithium adaptive cutoff radius in the lookup table is matched and selected. When the type combination of the current candidate neighbor atom pairs is silicon-silicon atom pairs, match and select the corresponding silicon-silicon adaptive cutoff radius in the lookup table; When the type combination of the current candidate neighbor atom pairs is a lithium-silicon atom pair, match and select the corresponding lithium-silicon adaptive cutoff radius in the lookup table; The extracted adaptive cutoff radius is used to perform initial spatial distance screening of the current candidate neighbor atom pairs.

5. The rapid molecular dynamics simulation method for lithium-silicon alloys as described in claim 3, characterized in that, Based on an adaptive truncation radius, candidate neighbor atom pairs are initially screened by spatial distance. Only atom pairs that meet the truncation criteria are retained. These atom pairs are then stored in a distant neighbor list, resulting in the distant neighbor list, which specifically includes: During the process of traversing candidate neighbor atom pairs, the current atom is calculated in real time. With atoms The three-dimensional Euclidean distance between them ; Then, the spatial distance The adaptive cutoff radius extracted from the lookup table corresponding to the current combination of atom types. Perform a geometric size comparison, wherein the truncation condition refers to satisfying a mathematical inequality. ; If the distance or squared distance of the current atomic pair satisfies the truncation condition, the current atomic pair is determined to have potential chemical bonding or effective bond level contribution in physics. The current atomic pair is retained as a valid far neighbor, and its global atomic index number and related topological information are written into the far neighbor list data structure used internally by ReaxFF. Otherwise, if the truncation condition is not met, the current atomic pair is determined to be an invalid atomic pair, and it is directly ignored and removed from the list.

6. The rapid molecular dynamics simulation method for lithium-silicon alloys as described in claim 3, characterized in that... In the chemical bond construction stage, an adaptive cutoff radius is used to perform a secondary screening of candidate atom pairs in the distant neighbor list to obtain the effective atom pairs after screening, thereby eliminating invalid atom pairs that do not contribute to the effective bond order. Specifically, this includes: In the chemical bond construction branch, each candidate atom pair in the distant neighbor list is read sequentially to obtain the real-time spatial distance between the two atoms. And retrieve the corresponding adaptive cutoff radius from the lookup table based on its atom type combination. ; Subsequently, the real-time spatial distance is determined. Does it meet the secondary screening criteria? ; If the conditions are met, the candidate atom pair is determined to be a valid atom pair and remains in an active state. The table entry data of the valid atom pair includes the local topological index and real-time distance of the lithium atom or silicon atom. Conversely, if If the distance between the candidate atoms is too great, their electron orbitals will not overlap effectively and their total bond level contribution will be close to 0. Therefore, they are considered invalid atom pairs with no effective bond level contribution.

7. The rapid molecular dynamics simulation method for lithium-silicon alloys as described in claim 3, characterized in that, The effective atom pairs are subjected to bond order determination using a bond order calculation function. Based on the determination results, bond energies, many-body terms, and non-bonded interactions are calculated. The atom states are iteratively updated through time integration until the simulation is complete. Specifically, this includes: In the energy calculation stage at each time step, the effective atom pairs retained after secondary screening are extracted one by one, and the real-time total bond level value is calculated by calling the bond level calculation function. ; After obtaining the current real-time total bond order value for each effective atom pair Then, dynamic physical corrections are made based on the local overcoordination of atoms, and the results are compared with a preset bonding cutoff threshold to determine whether the effective atom pair meets the minimum chemical bond interaction requirements. Based on the results of the bond order determination, the calculations of bond energy, many-body terms, and non-bonded interactions are carried out in parallel. The discrimination, solution and integration process is executed repeatedly in each time step until the preset total number of simulation time steps is reached, and the simulation is completed and the lithium silicon alloy structure data is output.

8. A rapid molecular dynamics simulation system for lithium-silicon alloys, characterized by: include: The model building module is configured to: construct an initial structural model of the lithium-silicon alloy under test based on the preset lithium-silicon ratio and lattice constant; The lookup table building module is configured to: build an adaptive truncation radius lookup table for atomic pairs, including: Analyze the reactive force field parameter file and extract the two-body parameters between different atom types; By utilizing the monotonically decreasing property of the bond order function of the reaction force field, under a preset bond order threshold, the critical distance corresponding to each atomic type pair is pre-calculated using a binary search method, and the critical distance is used as the adaptive cutoff radius of the corresponding atomic type pair. Based on the adaptive cutoff radius, an adaptive cutoff radius lookup table for atomic pairs is constructed; Utilizing the monotonically decreasing property of the bond order function of the reaction force field, under a preset bond order threshold, the critical distance corresponding to each atomic type pair is pre-calculated using a binary search method, and the critical distance is used as the adaptive cutoff radius for the corresponding atomic type pair. Specifically, this includes: For any atom type Using ReaxFF key-level functions Based on the monotonically decreasing property, an equivalence relationship is established between the bond level threshold and the geometric distance; ReaxFF bond level function. The specific expression is: ； in, These represent the contributions of the σ bond, π bond, and quadratic π bond correction terms, respectively. Each component is typically represented using an exponential decay function related to the interatomic distance: ； ； ； in, Represents atoms and The distance between them This refers to the characteristic distance parameter corresponding to the bond type. The first empirically fitted parameters are... The second empirical fitting parameter, The third empirical fitting parameter, The fourth empirical fitting parameter, The fifth empirical fitting parameter, The sixth empirical fitting parameter is used to regulate the decay behavior of the bond order as a function of distance. For a given key level threshold Key-level functions With geometric distance The value increases and then strictly decreases monotonically, thus there exists a unique critical distance. Satisfy the equation ; The critical distance is calculated by iteratively solving the problem within a set distance search interval using a binary search algorithm. The specific values ​​of are such that the following relationship holds: ； in, This represents a logically necessary and sufficient condition, if and only if the real-time geometric distance between the atomic pairs is... Less than the critical distance of this type of atom pair At that time, its calculated real-time key level Only then can it exceed the set key level threshold. ; For each pair of atom types, a binary search algorithm is used to search within the interval The inner solution satisfies The critical distance is used as the adaptive cutoff radius for the current atom type pair; where, This represents the set lower limit distance for the binary search, used to avoid non-physical minimum distance intervals where atoms are extremely overlapping, and to prevent overflow of the exponential term calculation; This represents the predefined geometric cutoff radius of the unified space of the ReaxFF reactive force field, which serves as the upper limit distance for the binary search. The molecular dynamics simulation module is configured to: perform accelerated reactive molecular dynamics simulation on the initial structural model of the lithium-silicon alloy using the atomic pair adaptive cutoff radius lookup table under preset temperature and simulation conditions, so as to perform relaxation and evolution of atomic positions and charge states, and obtain lithium-silicon alloy structural data after the simulation is completed. The output module is configured to output the atomic topological evolution process, bond order characteristics, or material stability of the lithium-silicon alloy under test, based on the lithium-silicon alloy structure data.

Images