Method for solving abnormal expansion of calculation box in molecular dynamics simulation process

Abstract

The application discloses a method for solving abnormal expansion of a calculation box in a molecular dynamics simulation process, and comprises the following steps: crystal configuration construction, determining the atomic arrangement structure characteristics and model size range of a calculation system, and constructing an initial crystal configuration; model setting, performing relevant model setting on the initial crystal configuration before MD calculation operation, and the model setting comprises the setting of boundary conditions, potential function selection, energy minimization, relaxation balance and output parameters of the calculation system; calculation operation, the MD method is mainly based on classical Newtonian mechanics to simulate the motion and distribution of atoms in the calculation system, that is, the distance between atoms determines the interaction force between atoms, which leads to the motion and redistribution of atoms, and the cycle is repeated until the atomic position and velocity information under certain conditions are obtained; and data analysis, according to the atomic velocity and position information obtained by calculation, the thermodynamic parameters and other macroscopic properties of the object system can be analyzed and obtained.

Description

Technical Field

[0001] This invention relates to the field of molecular dynamics simulation technology, and in particular to a method for solving the problem of abnormal enlargement of the calculation box during molecular dynamics simulation. Background Technology

[0002] Since the beginning of the 21st century, the widespread development and application of ultra-high-strength steel has effectively improved the lightweighting level and safety of equipment. However, with the increase in steel strength, ultra-high-strength steels such as hot-formed steel for automobiles, steel for bridge cables, and steel for hydrogen storage and transportation are prone to hydrogen embrittlement (i.e., delayed fracture caused by hydrogen in the material under stress), which seriously endangers the safety of ultra-high-strength steel components. Although scholars have been studying hydrogen embrittlement for over 100 years, there is still a lack of unified conclusions regarding the mechanism of hydrogen embrittlement in steel. Therefore, in-depth research on the diffusion and segregation of hydrogen atoms in the steel matrix under stress and their influence on crack formation remains essential, as this will help clarify the mechanism of hydrogen embrittlement in ultra-high-strength steel.

[0003] Given the small radius of hydrogen atoms and the very low hydrogen content in steel (typically less than 2 ppm), conventional experimental methods are often insufficient to effectively characterize hydrogen atoms in steel. Currently, the most advanced method for characterizing hydrogen atoms is cryogenic atom probe tomography (APT), but this technique can only statically analyze the local spatial distribution of hydrogen atoms at the nanoscale, and cannot dynamically characterize the diffusion and aggregation of hydrogen atoms under external stress. Furthermore, cryogenic APT is complex and costly, and very few institutions, both domestically and internationally, are proficient in this technique. Therefore, it is essential to explore other methods for studying hydrogen embrittlement.

[0004] The rapid development of computer technology has led to the widespread application of molecular dynamics (MD) methods in scientific exploration. MD methods rely on Newtonian mechanics to simulate the motion of atomic / molecular systems, enabling dynamic study at the atomic scale of the evolution of atomic distribution and structure in specific regions during material deformation and fracture (such as atomic diffusion, dislocation slip, and crack propagation). This gives it unparalleled advantages over experimental methods in simulating hydrogen atom diffusion and aggregation, dislocation slip accumulation, and crack initiation and propagation in ultra-high-strength steel under external stress. However, due to the small radius of hydrogen atoms, they exist in a free state in the steel matrix. During MD simulations, they easily diffuse and overflow from the steel matrix, causing an abnormal expansion of the "computation box" (i.e., the size range of the computational model constructed in MD), leading to calculation failures. This patent aims to solve this problem. Summary of the Invention

[0005] The purpose of this invention is to provide a method to address the above-mentioned shortcomings by providing a solution to the problem of abnormal expansion of the computation box during molecular dynamics simulation. This method solves the problem in the prior art where, due to the small radius of hydrogen atoms, they exist in a free state in the steel matrix and easily diffuse and overflow from the steel matrix during MD simulation, leading to abnormal expansion of the "computation box" (i.e., the size range of the computational model constructed in MD) and thus causing calculation failure.

[0006] This invention is achieved through the following scheme:

[0007] A method for resolving the abnormal scaling of the computational box during molecular dynamics simulations includes the following steps:

[0008] Step 1: Crystal configuration construction, determine the atomic arrangement structure characteristics and model size range of the computational system, and construct the initial crystal configuration;

[0009] Step 2: Model setup. Before running the MD calculation, the relevant model settings are performed on the initial crystal configuration. The model settings include the boundary conditions of the calculation system, selection of potential functions, energy minimization, relaxation equilibrium, shear deformation, and setting of output parameters.

[0010] Step 3: Calculation and execution. The MD method is mainly based on classical Newtonian mechanics to simulate and calculate the motion and distribution of atoms within the system. That is, the distance between atoms determines the interatomic force, which leads to the motion and redistribution of atoms. This process is repeated until the atomic position and velocity information under certain conditions is obtained.

[0011] Step 4: Data analysis. Based on the calculated atomic velocities and positions, the thermodynamic parameters and other macroscopic properties of the system can be analyzed.

[0012] In step one, the crystal configuration of the substrate is constructed so that the final state of the crystal configuration in the computer is a cuboid structure.

[0013] In step two, when adding the system to be studied to the matrix, it is only added to the interior of the matrix, while the surface area of ​​the matrix is ​​not added. This surface area serves as a "buffer zone" for the diffusion of the system to be studied.

[0014] In step three, during the relaxation equilibrium process of the crystal configuration, the atoms in the computational system will undergo diffusion motion, and some of the relevant systems to be studied will diffuse to the "buffer zone". After the relaxation equilibrium reaches a predetermined time period, the relevant systems to be studied in the "buffer zone" will be directly deleted.

[0015] In step three, the number of related systems to be studied that are directly deleted shall not exceed 5% of the total number of related systems to be studied initially added.

[0016] In step three, after deleting the relevant system to be studied, the continued relaxation equilibrium time of the system is calculated to be no less than 20% of the total relaxation equilibrium time; the “thickness” of the surface region should be no less than 5.

[0017] This scheme uses the LAMMPS software platform to perform MD calculations to simulate the shearing behavior of the bcc Fe-H system; and uses the conjugate gradient method to minimize the energy of the constructed initial crystal configuration under 0K conditions.

[0018] When setting up the model on the LAMMPS software platform, the boundary conditions in the X, Y and Z axis directions are all set as shrinking boundary conditions (s-boundaries). The model size (or model boundary position) of the s-boundary will expand or shrink as the movement range of all atoms in the computational system increases, ensuring that all atoms are within the model boundary range.

[0019] After setting the system boundary conditions, this scheme allows for the energy minimization and relaxation equilibrium of the initial crystal configuration by selecting an appropriate potential function based on the research object; the MEAM type Fe-H potential function created by Lee-Jang is selected.

[0020] The calculation system in this scheme uses the NPT ensemble for relaxation equilibrium at a temperature of 300K and a relaxation equilibrium time of 100ps; shearing is performed in the NVT ensemble at a temperature of 300K and a shear strain rate of 10. 8 / s. In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are:

[0021] (1) The proposed methods of “not adding the relevant system to be studied to the surface part of the matrix” and “directly deleting the relevant system to be studied in the “buffer zone”” can ensure the normal diffusion distribution evolution of the relevant system to be studied in the final crystal configuration during the shearing process. Until the shearing ends (800ps), the relevant system to be studied will not diffuse outside the matrix, and there is no abnormal expansion of the calculation box. Finally, the calculation simulation is successfully completed. Attached Figure Description

[0022] Figure 1 It has a bcc Fe matrix crystal configuration;

[0023] Figure 2 This is a diagram showing the initial distribution of H atoms in the bcc Fe matrix;

[0024] Figure 3 This is a diagram showing the evolution of H atom diffusion distribution in the bcc Fe-H crystal configuration during relaxation equilibrium.

[0025] Figure 4 Add a region map to the improved H atom;

[0026] Figure 5The initial H atom distribution diagram in the improved bcc Fe matrix;

[0027] Figure 6 The evolution of H atom diffusion distribution during relaxation equilibrium in the improved bcc Fe-H crystal configuration;

[0028] Figure 7 The evolution of H atom diffusion distribution during shearing in the improved bcc Fe-H crystal configuration; Detailed Implementation

[0029] All features disclosed in this specification, or all steps in all disclosed methods or processes, may be combined in any way, except for mutually exclusive features and / or steps.

[0030] Any feature disclosed in this specification (including any appended claims and abstract) may be replaced by other equivalent or similar features, unless specifically stated otherwise. That is, unless specifically stated otherwise, each feature is merely one example of a series of equivalent or similar features.

[0031] In the description of this invention, it should be understood that the terms "upper", "lower", "left", "right", etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this invention.

[0032] Furthermore, the terms "first," "second," etc., are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined with "first," "second," etc., may explicitly or implicitly include one or more of that feature.

[0033] Example 1

[0034] This embodiment uses the crystal configuration of bcc Fe matrix as an example for explanation and illustration. The relevant system to be studied uses H atom as an example for explanation and illustration. It should be noted that the relevant system to be studied is determined by the influence of a single element or a combination of multiple elements on the target crystal configuration. It can include any type of element to be studied.

[0035] Meanwhile, the crystal configuration of the matrix can be the crystal structure of common pure metals, including: face-centered cubic (FCC) structure, common metals with FCC structure include β-iron, gold, silver, copper, etc.; or body-centered cubic (BCC) structure, common metals with BCC structure include α-iron, molybdenum, tungsten, vanadium, etc.; or hexagonal close-packed (HCP) structure, common metals with HCP structure include cobalt, zinc, magnesium, etc.; or simple cubic structure: polonium.

[0036] This approach primarily targets multi-element computational systems, rather than single-element systems; that is, the computational system must contain at least two types of atoms of different element types. Taking the bcc Fe-H system as an example, when constructing the initial crystal model, the system is... Figure 1 After establishing the bcc Fe matrix crystal structure in (a), H atoms need to be added to the central region (bcc Fe) (the two side regions are shear-fixed layers and do not require the addition of H atoms) to study the motion behavior of H atoms during shearing. In similar related studies, non-matrix atoms (in this patent, the matrix atoms are Fe atoms, and the non-matrix atoms are H atoms) are generally randomly and uniformly added to the entire central region according to a certain proportion, i.e., as shown in the figure. Figure 2 As shown.

[0037] After minimizing the energy of the constructed bcc Fe-H crystal configuration, relaxation equilibrium was achieved. However, at 40 ps relaxation, it was found that some H atoms had diffused out of the bcc Fe matrix, causing an abnormal expansion of the computational model size, i.e., an abnormal expansion of the computational box (in MD, the region containing all computational atoms is called the computational box), such as... Figure 3 As shown in (c), as the relaxation equilibrium continues, H atoms continue to diffuse away from the bcc Fe matrix region, causing the computational box to expand abruptly further, as... Figure 3 As shown in (d), an abnormal expansion of the computation box can cause MD calculations to become difficult to continue normally, resulting in the failure of the computational simulation.

[0038] in Figure 1 (a) is a full diagram of the crystal structure of the bcc Fe matrix (the inset is a magnified view of a part to show the arrangement and distribution of Fe atoms); Figure 1 (b) Figure 1 (c) and Figure 1 (d) are all partial cross-sectional views to fully show the crystal configuration; in the figures, the filling parts are all Fe atoms; the central region is the shear region, and the two side regions are the fixed layers;

[0039] exist Figure 2 Note: In the image above, the spheres represent H atoms; the areas where H atoms are added are... Figure 1(a) Central region; all Fe atoms have been hidden in the above figure to better illustrate the distribution of H atoms;

[0040] Figure 3 This describes the evolution of H atom diffusion distribution during relaxation equilibrium. Figure 3 (c) and Figure 3 (d) The black circles represent H atoms that diffused out of the bcc Fe matrix;

[0041] exist Figure 4 The bcc Fe matrix crystal configuration and Figure 1 Completely identical, except that no H atoms are added to the surface layer of the central region;

[0042] First, let's introduce the relevant terminology:

[0043] Minimizing energy refers to the process of rearranging some atoms in a local region within a computational system. Since the initial crystal configuration may have some disordered atomic arrangement in local regions, the entire computational system is not in the minimum energy state, that is, the computational system is not in the most stable state. Therefore, it is necessary to minimize the energy of the initial crystal configuration first.

[0044] like Figure 1 As shown, the present invention provides a technical solution:

[0045] A method for resolving the abnormal scaling of the computational box during molecular dynamics simulations, comprising the following steps:

[0046] Step 1: Crystal Configuration Construction. When using the MD method for related research, the first step is to determine the atomic arrangement structure characteristics and model size range of the computational system based on the characteristics of the research object and research objectives, thereby constructing the initial crystal configuration. Specifically, constructing the crystal configuration of the bcc Fe matrix requires that the arrangement of Fe atoms in the crystal configuration satisfies a body-centered cubic (bcc) structure, such as... Figure 1 As shown; the initial dimensions of this crystal configuration in the X, Y and Z axis directions are 21.1 nm, 20.1 nm and 12.4 nm, respectively, with a total of 450,000 Fe atoms;

[0047] In the aforementioned method, H atoms are randomly and uniformly added to Figure 1 (a) The entire bcc Fe matrix (central region) is included. This approach may result in H atoms being added to the surface edges of the bcc Fe matrix (i.e., the edges of the calculation box), such as... Figure 2As shown, the initial bcc Fe-H crystal configuration constructed in this way is highly likely to cause H atoms to diffuse out of the bcc Fe matrix during subsequent relaxation equilibrium processes due to the random diffusion of H atoms, ultimately leading to an abnormal expansion of the computational box. To solve this problem, this patent restricts and improves the range of H atom addition regions;

[0048] Specifically, once the construction is complete, as follows: Figure 1 (a) shows the bcc Fe crystal configuration. When adding H atoms to the bcc Fe matrix (central region), the addition is only made inside the bcc Fe matrix (central region), while the surface portion of the bcc Fe matrix (central region) is not added. This surface portion serves as a "buffer zone" for H atom diffusion, i.e., as shown in the diagram. Figure 4 As shown, no H atoms are added to the surface layer; Figure 5 The initial crystal configuration of bcc Fe-H constructed using this method is compared with... Figure 5 (d) and Figure 2 (d) It can be found that, Figure 5 (d) No H atoms are present in the surface region. This method ensures that even if H atoms are randomly added to the initial bcc Fe-H crystal configuration, the H atoms will remain intact. Figure 4 At the outer edge of the bcc Fe matrix (central region) (i.e., the junction between the bcc Fe matrix (central region) and the surface region), during the subsequent energy minimization and relaxation equilibrium process, when H atoms diffuse outward at this edge, they will remain within the surface region and will not diffuse out of the bcc Fe matrix. In short, the surface region acts as a "buffer zone" for H atom diffusion, preventing hydrogen atoms from diffusing out of the bcc Fe matrix and thus avoiding abnormal expansion of the computational box. In actual simulation calculations, the "thickness" of the surface region needs to be determined comprehensively based on factors such as the atom type in the computational system, simulation duration, and model size; this patent suggests that the "thickness" of the surface region should not be less than 5.

[0049] Step 2: Model Setup. For molecular dynamics simulations, calculations can be performed either by self-programming or by utilizing numerous existing mature commercial software platforms. Before running MD calculations, the initial crystal configuration needs to be set up, which involves writing scripts that meet the calculation rules of the relevant software platforms to achieve specific research objectives. Model setup mainly includes setting the boundary conditions of the calculation system, potential function selection, energy minimization, relaxation equilibrium, shear deformation, output parameters, etc.

[0050] During the relaxation equilibrium process, atoms in the computational system undergo diffusion motion. Figure 3In this case, H atoms diffuse within the bcc Fe matrix during relaxation equilibrium until they overflow outside the matrix, causing the computational box to expand abnormally and the simulation to fail. Similarly, although... Figure 4 The middle layer constructs a "buffer zone" (surface region) for the diffusion of H atoms, but the H atoms initially added near the surface region may diffuse into the surface region during the relaxation equilibrium process, or even eventually diffuse out of the surface region.

[0051] Figure 6 To observe the diffusion distribution of H atoms in the improved bcc Fe-H crystal configuration during relaxation equilibrium, it can be found that at 80 ps of relaxation equilibrium, a large number of H atoms have diffused into the surface region ("buffer zone"). Figure 6 (b) H atoms outside the dashed box (indicated by small arrows; see enlarged view). Figure 6 (e) H atoms in the "buffer zone" may continue to diffuse outwards beyond the bcc Fe matrix during subsequent relaxation equilibrium (depending on factors such as relaxation equilibrium time, atom type, and "buffer zone" thickness). Furthermore, even if H atoms do not diffuse outwards beyond the bcc Fe matrix during relaxation equilibrium, they may diffuse outwards during subsequent shearing (or other research needs, such as MD simulations of stretching, impact, and nanoindentation), ultimately leading to abnormal expansion of the computational box and simulation failure. Therefore, H atoms diffused into the "buffer zone" require certain processing. This patent proposes to directly delete H atoms in the "buffer zone" after a certain relaxation equilibrium period, i.e., direct deletion. Figure 6 (b) All H atoms outside the dashed box (indicated by the small black arrows); after deleting the H atoms in the "buffer zone", the distribution of H atoms in the calculated system is as follows. Figure 6 As shown in (c), this prevents H atoms from diffusing out of the bcc Fe matrix during subsequent relaxation equilibrium and shearing processes, thus avoiding abnormal expansion of the computation box.

[0052] It should be noted that in MD computational simulations, directly deleting some H atoms will affect the dynamic equilibrium of local regions of the computational system. Therefore, after deleting some H atoms, the computational system needs to continue to relax and reach equilibrium. In addition, in order to ensure the stability of the computational system and the accuracy of the calculation results, this patent suggests that the number of H atoms directly deleted should not exceed 5% of the total number of H atoms initially added to the crystal configuration, and the relaxation and equilibrium time of the computational system after deleting H atoms should not be less than 20% of the total relaxation and equilibrium time.

[0053] The two improvement methods mentioned above can be used alone or in combination, depending on the actual research object and research objectives. Figure 7To combine the two improved methods mentioned above, the diffusion distribution evolution of H atoms in the bcc Fe-H crystal configuration obtained during the shearing process was observed. Until the shearing ended (800 ps), H atoms were all distributed within the bcc Fe matrix, and there was no phenomenon of H atoms diffusing out of the bcc Fe matrix. In other words, there was no abnormal expansion of the computational box, and the computational simulation was successfully completed.

[0054] Furthermore, this patent only uses the bcc Fe-H system as an example; it is still applicable to other multi-atom type calculation systems.

[0055] Specifically, this scheme uses the LAMMPS (Large-scale Atomic / Molecular Massively Parallel Simulator) software platform to perform MD calculations to simulate the shear behavior of the bcc Fe-H system;

[0056] When setting up the model on the LAMMPS software platform, the boundary conditions in the X, Y and Z axis directions are all set as shrinking boundary conditions (s boundary). The model size (or model boundary position) of the s boundary will expand or shrink with the range of motion of all atoms in the computation system to ensure that all atoms are within the model boundary range.

[0057] After setting the boundary conditions of the computational system, select an appropriate potential function according to the research object (the potential function is an expression function used to calculate the interaction force between atoms. The potential function is different for different research objects or atomic systems; select the MEAM type Fe-H potential function created by Lee-Jang), and then perform energy minimization and relaxation equilibrium on the initial crystal configuration.

[0058] This scheme utilizes the conjugate gradient method (CG method) to minimize the energy of the initial crystal configuration under 0K conditions.

[0059] During the energy minimization process, the atoms in the initial crystal configuration can only undergo a small rearrangement within a limited range at 0 K, and this is only considered from the perspective of minimizing the energy of the computational system. In reality, atoms experience thermal perturbations at certain temperatures, meaning they can undergo long-range or short-range diffusion within a certain range. The positions of atoms are not fixed, and the positional relationships between atoms do not necessarily strictly conform to the atomic arrangement rules of a specific crystal structure (such as the bcc structure). Therefore, after minimizing the energy, the initial crystal model needs to undergo relaxation equilibrium. The relaxation equilibrium process involves all atoms in the computational system undergoing sufficient diffusion under certain ensemble conditions, causing the computational system to reach a quasi-equilibrium steady state that approximates the actual state. In this scheme, the relaxation equilibrium of the computational system uses an isothermal-pressure, constant-temperature (NPT) ensemble at 300 K and a relaxation equilibrium time of 100 ps.

[0060] After the system has reached relaxation equilibrium, the atomic arrangement in the crystal configuration is in a near-steady state, and shearing can then be applied to the crystal configuration. In this scheme, shearing is performed under a canonical ensemble (NVT) at a temperature of 300 K and a shear strain rate of 10⁻⁶. 8 / s.

[0061] In the LAMMPS software platform, crystal configuration shearing can be achieved through the following settings: Figure 1 In (a), both the left and right sides of the Y-axis are fixed layers; during shearing, the right fixed layer (i.e., Figure 1 (a) The region indicated by the arrow moves along the +X axis at a certain speed (0.0171 ps in this patent), while the left fixed layer remains stationary with a movement speed of 0. The movement speed of the middle region increases linearly from 0 to 0.0171 ps along the +Y axis. In this way, an approximate shearing effect can be obtained. It should be noted that during the energy minimization, relaxation equilibrium, and shearing process, the relative positions of all Fe atoms in the regions on both sides of the Y axis (fixed layer) remain unchanged. That is, there is no diffusion motion of atoms in the regions on both sides, and only the right fixed layer undergoes overall movement.

[0062] Step 3: Computational Execution. The MD method primarily uses classical Newtonian mechanics to simulate the motion and distribution of atoms within a system. The distance between atoms determines the interatomic forces, leading to atomic motion and redistribution. This process repeats until the atomic positions and velocities under specific conditions are obtained. In short, the MD computational execution involves using a computer to solve for the forces and motions between atoms in the system based on the crystal configuration and model settings, until the atomic velocities and positions under the target conditions are finally obtained; this computational process is entirely computer-driven.

[0063] Step 4: Data analysis. Based on the calculated atomic velocities and positions, the thermodynamic parameters and other macroscopic properties of the target system can be analyzed and obtained.

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

Claims

1. A method for resolving the abnormal enlargement of the computational box during molecular dynamics simulations, characterized in that: The steps are as follows: Step 1: Crystal configuration construction, determine the atomic arrangement structure characteristics and model size range of the computational system, and construct the initial crystal configuration; In step one, the crystal configuration of the matrix is ​​constructed so that the final state of the crystal configuration in the computer is a cuboid structure; Step 2: Model setup. Before running the MD calculation, the relevant model settings are performed for the initial crystal configuration. The model settings include the boundary conditions of the calculation system, the selection of potential functions, energy minimization, relaxation equilibrium, and the setting of output parameters. In step two, when adding the system to be studied to the matrix, it is only added to the interior of the matrix, while the surface layer of the matrix is ​​not added. This surface layer serves as a "buffer zone" for the diffusion of the system to be studied. Step 3: Calculation and execution. The MD method is based on classical Newtonian mechanics to simulate and calculate the motion and distribution of atoms within the system. That is, the distance between atoms determines the force between atoms, which leads to the motion and redistribution of atoms. This process is repeated until the atomic position and velocity information under certain conditions is obtained. Step 4: Data analysis. Based on the calculated atomic velocities and positions, the thermodynamic parameters and other macroscopic properties of the system can be analyzed.

2. The method for solving the problem of abnormal expansion of the computational box during molecular dynamics simulation as described in claim 1, characterized in that: In step three, during the relaxation equilibrium process of the crystal configuration, the atoms in the computational system will undergo diffusion motion, and some of the relevant systems to be studied will diffuse to the "buffer zone". After the relaxation equilibrium reaches a predetermined time period, the relevant systems to be studied in the "buffer zone" will be directly deleted.

3. The method for solving the problem of abnormal scaling of the computational box during molecular dynamics simulation as described in claim 2, characterized in that: In step three, the number of related systems to be studied that are directly deleted shall not exceed 5% of the total number of related systems to be studied initially added.

4. The method for solving the problem of abnormal scaling of the computational box during molecular dynamics simulation as described in claim 3, characterized in that: In step three, after deleting the relevant system to be studied, the calculated continued relaxation equilibrium time of the system should not be less than 20% of the total relaxation equilibrium time; the "thickness" of the surface region should not be less than 5 Å.

5. The method for solving the problem of abnormal scaling of the computational box during molecular dynamics simulation as described in claim 1, characterized in that: The LAMMPS software platform was used to realize the shear behavior of the system in MD computational simulation; the energy of the initial crystal configuration was minimized under 0K conditions using the conjugate gradient method.

6. The method for solving the problem of abnormal scaling of the computational box during molecular dynamics simulation as described in claim 5, characterized in that: When setting up the model on the LAMMPS software platform, the boundary conditions in the X, Y and Z axis directions are all set as shrinking boundary conditions, specifically s-boundaries. The model size or model boundary position of the s-boundary will expand or shrink with the range of motion of all atoms in the computational system, ensuring that all atoms are within the model boundary range.

7. The method for solving the problem of abnormal scaling of the computational box during molecular dynamics simulation as described in claim 6, characterized in that: After setting the system boundary conditions, and selecting an appropriate potential function according to the research object, the initial crystal configuration can be minimized and relaxed to achieve equilibrium. The MEAM type Fe-H potential function created by Lee-Jang is selected.

8. The method for solving the problem of abnormal expansion of the computational box during molecular dynamics simulation as described in claim 4, characterized in that: The system was brought to relaxation equilibrium using an NPT ensemble at 300 K and a relaxation time of 100 ps. Shearing was performed in an NVT ensemble at 300 K and a shear strain rate of 10. 8 / s.

Images