Improved method and apparatus for computer simulation of substances and materials

Abstract

The application provides an improved method and device for computer simulation of substances and materials, which firstly constructs a multi-particle model and system of the substances and materials to be simulated; characterized in that: a multi-particle system dynamics equation set containing a Nose-Hoover⊕Langevin equation set is used to calculate the change of the coordinates and velocities of all particles in the multi-particle model of the substances and materials with time, the Nose-Hoover⊕Langevin equation set gives each particle in the multi-particle system three additional dynamic quantities, and the additional dynamic quantities of each particle are set as the thermodynamic and statistical mechanical quantities corresponding to the eccentricity of the particle; and thermodynamic and dynamic analysis is performed on the corresponding substances and materials according to the coordinates and velocities of all particles in the multi-particle model of the substances and materials at each discrete time step.

Description

Technical Field

[0001] This invention belongs to the research field of condensed matter and materials science and computer simulation engineering, and specifically relates to an improved technical solution for computer numerical simulation based on the accelerated molecular dynamics method of condensed matter multi-particle systems. Background Technology

[0002] Molecular dynamics numerical simulation is currently one of the most important simulation methods in the field of materials research. Because it can directly simulate the real-time dynamic processes of matter at the atomic and molecular scale, it is widely used in materials science, physics, chemistry, life sciences, and other related technical fields. The basic principles of commonly used traditional molecular dynamics simulation methods can be described as follows.

[0003] For a multi-particle system of any substance or material containing N particles (atoms or molecules), molecular dynamics simulation methods generally use the coordinates of all particles in three-dimensional Euclidean space X = (x1, x2, ..., xn). N ) and the corresponding particle momentum P=(p1,p2,…,p N (p) i =m i v i m i Let v be the mass of particle i. i Let the velocity of particle i (i∈{1,2,…,N}) be used as a dynamic quantity and constitute a phase space to describe any microstate of the multi-particle system. Then, a specific set of dynamic equations is used to describe the dynamic process of the evolution of the microstate of the multi-particle system over time. In traditional molecular dynamics simulation methods, commonly used sets of dynamic equations include (see reference [1]):

[0004] A. Newtonian dynamics (1687)

[0005] F i =m i a i

[0006] B. Langevin's Dynamics (1908)

[0007]

[0008] C. Nose-Hoover extended dynamics (1984)

[0009]

[0010] In the above formula, the physical quantities corresponding to each symbol are as follows:

[0011] F—The force acting on the particle (generally a conservative force);

[0012] m—mass of the particle;

[0013] v—the velocity of the particle ( x is the coordinate of the particle in three-dimensional Euclidean space;

[0014] a—the acceleration of the particle ( x is the coordinate of the particle in three-dimensional Euclidean space;

[0015] γ—the viscosity coefficient of particles moving in the Langevin kinetic hypothetical solution;

[0016] T—The temperature of the heat bath for the multi-particle system;

[0017] R(t) — White noise function;

[0018] k B —Boltzmann constant;

[0019] ζ—Nose-Hoover dynamics extension variable for multi-particle systems;

[0020] Q—the mass corresponding to the Nose-Hoover dynamic extension variable ζ;

[0021] N—The number of all particles (atoms or molecules) in the system.

[0022] In all physical quantity symbols, the subscript i indicates that the physical quantity belongs to any particle labeled i in a multi-particle system. The single dot and double dot on the physical quantity symbol represent the first derivative and the second derivative with respect to time, respectively.

[0023] Within the aforementioned model and methodological framework, to realistically and reliably reflect the motion of atoms and molecules in matter, the time step Δt after numerical discretization of the dynamic equations is generally no more than a few femtoseconds (approximately 5% of the vibration period of atoms in a typical solid lattice). This limits the simulation capabilities of current computers to simulating physical processes on a nanosecond timescale. However, the microstructural evolution processes involved in commonly studied physical phenomena, such as thin film growth, solid-state phase transitions, material deformation, and heat treatment, involve time spans exceeding milliseconds, or even days or years (e.g., material creep). For example, in iron-carbon alloys like steel, the time required for a carbon atom to hop once between iron atoms at 500°C is approximately one microsecond. Such a vast timescale limitation (difference of approximately 3 to 13 orders of magnitude) makes it difficult to employ traditional molecular dynamics simulation methods to study the physical processes and their laws occurring in matter and materials under normal conditions.

[0024] Although some open-source software (see references [2-4]) and commercial software (see references [5-8]) have incorporated molecular dynamics simulation methods into their functionality, these methods are still largely confined to academic laboratories and used to a limited extent by researchers. The aforementioned timescale limitation is the key bottleneck hindering its wider application in industrial and engineering practice. If this bottleneck can be resolved and overcome, it is foreseeable that it will bring a new look to the research and development methods of related industrial fields, especially the materials industry. At that time, molecular dynamics computer numerical simulation methods will become an essential tool for materials research, similar to electron microscopy. Furthermore, due to the significant cost advantage of high-performance computers compared to expensive experimental equipment, the time cycle and economic cost of developing new materials and materials processing technologies will be greatly reduced, thereby effectively promoting the development of various industrial and engineering technologies.

[0025] To address this bottleneck on the timescale, numerous research teams worldwide have conducted continuous and tireless research, developing and proposing several solutions for so-called accelerated molecular dynamics simulations. These solutions can be broadly categorized into three types based on their technical approaches: ① Simulation methods for raising the system's potential energy surface, including Hyperdynamics (see references [9,10]), Metadynamics (see references [11,12]), Strain-boost molecular dynamics (see reference

[13] ), and Adaptive boost molecular dynamics (see reference

[14] ). The drawback of this type of method is that, since the potential energy surface of a multi-particle system is a hypersurface in 3N-dimensional space, the morphology of this hypersurface becomes extremely complex as the number of particles N increases. Therefore, careful attention must be paid to whether the superposition of bias potentials violates the fundamental laws of chemical reaction dynamics. At the same time, there are no fixed rules for constructing bias potentials, and they are generally quite complex, limiting their widespread application to specific dynamic processes such as atomic diffusion, dislocation nucleation and motion. ② Multi-replica simulation methods, mainly including Parallel ReplicaDynamics (see references [10,15]) and Parallel Trajectory Splicing (see reference

[16] ). When using these methods, because the same molecular dynamics simulation is performed for each replica, they are only suitable for models of material systems with a small number of particles, thus limiting their widespread application. ③ High-temperature-low-temperature mapping simulation methods, mainly including Temperature Accelerated Dynamics methods and models (see references [10,17,18]). The main drawback of this method is that, since the microscopic state changes of the same system at higher and lower temperatures can differ significantly, molecular dynamics simulations at higher temperatures may not effectively and reasonably represent the microscopic state changes at lower temperatures, potentially leading to kinetic distortion.

[0026] In addition to the three technical solutions mentioned above, by drawing on the idea of ​​dynamical extended variables proposed by Shuichi Nose in developing Nose-Hoover dynamics (see references [19,20]), and based on the concept of "collective variables" in multi-particle dynamics systems introduced to describe the dynamic processes of specific microstructure transition events in condensed matter multi-particle systems (see references [11,21]), Luca Maragliano and Eric Vanden-Eijnden published an article in 2006 (see reference

[22] ) proposing that the thermodynamic and statistical mechanical quantities corresponding to the collective variables of multi-particle dynamics systems can be set as dynamical extended variables, and then a set of extended Langevin dynamic equations can be used to describe the dynamic processes of the evolution of the system's microstates, thereby achieving an accelerated effect in molecular dynamics simulation of specific microstructure transition events corresponding to the definition expressions of the collective variables. The effectiveness and reliability of this method strongly depend on the rationality of the construction method of the collective variables (see references [21,23]). Although many scholars and researchers have tried to improve the construction of set variables within this framework (such as the use of artificial intelligence methods, see references [24,25]) to accelerate molecular dynamics simulation, so far there is no reasonable, reliable, universal, simple and effective technical solution to support the achievement of good time scale crossing effect in the computer numerical simulation of molecular dynamics of multi-particle condensed matter systems.

[0027] The applicant's research team filed a patent application in August 2020 entitled "A Computer Simulation Method and Apparatus for Matter and Materials," which has been granted (patent number 202010884514.3). This patent includes: constructing a multi-particle model of the matter and materials system to be simulated in a computer simulation; describing the microstate of the system using the coordinates and momentum of all particles (atoms or molecules); setting a particle stepping motion metric D for each particle in the multi-particle model to describe the reaction coordinates of any microstructural transformation process that can occur in the model; setting the initial reference state of the multi-particle model; calculating the evolution of the microstate of the matter and materials system over time using the four-dimensional extended Langevin kinetic equations of the multi-particle system, where each particle is assigned an additional kinetic quantity s, and s is coupled to D in a harmonic oscillator manner; and performing thermodynamic and kinetic analysis on the corresponding matter and materials system based on the obtained coordinates and velocities of all particles in the model at each discrete time step.

[0028] In May 2021, the applicant's research team filed a patent application entitled "A Novel Computer Simulation Method and Apparatus for Matter and Materials," patent application number 202110528679.1. This patent application improves upon the method in the aforementioned patent (patent number 202010884514.3) by introducing a new particle nearest neighbor eccentricity metric D to describe the degree of non-cooperative motion between an arbitrary particle and its nearest neighbor particles, replacing the particle shuffling motion set in the method of the aforementioned patent (patent number 202010884514.3). After replacing the metric D with this setting, it is not necessary to set the reference state of the multi-particle system. Therefore, it is more convenient to implement and use than the method given in the patent with patent number 202010884514.3. At the same time, since the calculation formula of the particle nearest neighbor eccentricity metric D is simpler than the calculation formula of the particle shuffling motion metric D set in the method of patent number 202010884514.3, it is easier to program in actual computer simulation and has better computational efficiency.

[0029] In July 2021, the applicant's research team filed a patent application entitled "Computer Simulation Method and Apparatus for Matter and Materials Across Time Scales," patent application number 202110778355.3. This patent application further improves upon the method in the aforementioned patent application number 202110528679.1 by introducing a new particle nearest neighbor eccentricity metric D to describe the degree of non-cooperative motion between any particle and its nearest neighbor particles. This new metric replaces the particle nearest neighbor eccentricity metric D in the method of patent number 202110528679.1. Since the newly introduced particle nearest neighbor eccentricity metric D is a linear function of particle displacement, while the particle nearest neighbor eccentricity metric D in the method of patent number 202110528679.1 is a non-linear function of particle displacement, the new patent method can more accurately describe and reflect the reaction coordinates of the microstructure transformation process in the multi-particle model. Therefore, when performing computer simulations on actual matter and materials, higher fidelity computer simulation data results can be obtained.

[0030] Further testing and research by the applicant's research team revealed some shortcomings in the practical application of the patented method (patent number 202110778355.3). For example, because it uses the four-dimensional extended Langevin kinetic equations of a multi-particle system to calculate the evolution of the microscopic state of a matter or material system over time, it is difficult to accurately control the temperature of the multi-particle model of the simulated matter or material when using this patented method for computer simulation. In practical use, a series of tests and runs are required beforehand to continuously adjust the value of the main heat bath temperature parameter in the four-dimensional extended Langevin kinetic equations of the multi-particle system to determine a suitable value for the parameter in order to accurately control the temperature of the simulated multi-particle model of the matter or material. This brings inconvenience to the implementation and use of the patented method.

[0031] Related literature for this invention:

[0032] [1] Chen Minbo. Computational Chemistry: From Theoretical Chemistry to Molecular Simulation. Beijing: Science Press, 2009.

[0033] [2]LAMMPS Molecular Dynamics Simulator.https: / / lammps.sandia.gov / .

[0034] [3]NAMD-Scalable Molecular Dynamics.http: / / www.ks.uiuc.edu / Research / namd / .

[0035] [4]GROMACS.http: / / www.gromacs.org / .

[0036] [5]BIOVIA Materials Studio.https: / / 3dsbiovia.com / products / collaborative-science / biovia-materials-studio / .

[0037] [6]QuantumATK.https: / / www.synopsys.com / silicon / quantumatk.html.

[0038] [7]DL_POLY Molecular Simulation Package.https: / / www.scd.stfc.ac.uk / Pages / DL_POLY.aspx.

[0039] [8]The Vienna Ab initio Simulation Package.https: / / www.vasp.at / .

[0040] [9]VOTER A F.Hyperdynamics:Accelerated Molecular Dynamics ofInfrequent Events.Phys Rev Lett,1997,78(20):3908-11.

[0041]

[10] VOTER AF,MONTALENTI F,GERMANN T C.Extending the timescale inatomistic simulation of materials.Annual Review of Materials Research,2002,32(1):321-46.

[0042]

[11] LAIO A,PARRINELLO M.Escaping free-energy minima.Proceedings ofthe National Academy of Sciences of the United States of America,2002,99(20):12562-6.

[0043]

[12] LAIO A,GERVASIO F L.Metadynamics:a method to simulate rare eventsand reconstruct the free energy in biophysics,chemistry and materialscience.Rep Prog Phys,2008,71(12):126601.

[0044]

[13] HARA S,LI J.Adaptive strain-boost hyperdynamics simulations ofstress-driven atomic processes.Physical Review B,2010,82(18):184114.

[0045]

[14] ISHII A,OGATA S,KIMIZUKA H,et al.Adaptive-boost moleculardynamics simulation of carbon diffusion in iron.Physical Review B,2012,85(6):064303.

[0046]

[15] VOTER A F.Parallel replica method for dynamics of infrequentevents.Physical review B,Condensed matter,1998,57(22):R13985-8.

[0047]

[16] PEREZ D,CUBUK E D,WATERLAND A,et al.Long-time dynamics throughparallel trajectory splicing.Journal of Chemical Theory&Computation,2015,acs.jctc.5b00916.

[0048]

[17] SORENSEN M R,VOTER A F J.Temperature-accelerated dynamics forsimulation of infrequent events.Journal of Chemical Physics,2000,112(21):9599-606.

[0049]

[18] ZAMORA R J,UBERUAGA B P,PEREZ D,et al.The Modern Temperature-Accelerated Dynamics Approach.Annu Rev Chem Biomol Eng,2016,7(1):87-110.

[0050]

[19] NOSE S.A Unified Formulation of the Constant TemperatureMolecular-Dynamics Methods.Journal of Chemical Physics,1984,81(1):511-9.

[0051]

[20] NOSE S.A molecular dynamics method for simulations in thecanonical ensemble.Molecular Physics,1984,52(2):255-68.

[0052]

[21] FIORIN G,KLEIN M L.Using collective variables to drive moleculardynamics simulations.Molecular Physics,2013,111(22-23):3345-62.

[0053]

[22] MARAGLIANO L,VANDEN-EIJNDEN E.A temperature accelerated methodfor sampling free energy and determining reaction pathways in rare eventssimulations.Chemical Physics Letters,2006,426(1-3):168-75.

[0054]

[23] VALSSON O,TIWARY P,PARRINELLO M.Enhancing Important Fluctuations:Rare Events and Metadynamics from a Conceptual Viewpoint.Annual review ofphysical chemistry,2016,67(1):annurev-physchem-040215-112229.

[24] ZHANG J,CHENM.Unfolding Hidden Barriers by Active Enhanced Sampling.Physical ReviewLetters,2017,121(1):010601.

[0055]

[25] CHEN W,TAN AR,FERGUSON A L.Collective variable discovery and enhanced sampling using autoencoders:Innovations in network architecture and error function design.Journal of Chemical Physics,2018,149(7):072312.

[0056]

[26] D. FRENKEL et al. Molecular Simulation: From Algorithms to Applications. Beijing: Chemical Industry Press, 2004.

[0057]

[27] Interatomic Potentials Repository.https: / / www.ctcms.nist.gov / potentials / .

[0058]

[28] BECQUART CS,RAULOT JM,BENCTEUX G,et al.Atomistic modeling of anFe system with a small concentration of C.Computational Materials Science,2007,40(1):119-29.

[0059]

[29] VEIGA RG, PEREZ M, BECQUART CS, et al. Atomistic modeling of carbonCottrell atmospheres in bcc iron. Journal of physics Condensed matter, 2013, 25(2):025401.

[0060]

[30] APOSTOL F, MISHIN Y. Angular-dependent interatomic potential for the aluminum-hydrogen system. Physical Review B, 2010, 82(14): 144115. Summary of the Invention

[0061] The purpose of this invention is to address the problems existing in the prior art by proposing a new technical solution to simulate the evolution of the microscopic state of multi-particle (atomic or molecular) systems corresponding to matter and materials over time on a computer, so as to achieve highly reliable, efficient, versatile and more convenient computer simulation analysis of matter and materials.

[0062] The technical solution of this invention provides an improved method for computer simulation of substances and materials, which first involves the following steps: 1.

[0063] Step 1: For the substance and material to be simulated in computer simulation, construct a multi-particle model and system for that substance and material, as follows:

[0064] For any substance or material containing N particles, the set of all particles is denoted as {1,2,…,N}. All particles are represented in three-dimensional Euclidean space. The coordinates in the figure are X = (x1, x2, ..., x). N ) and the corresponding particle momentum P=(p1,p2,…,p N As a dynamic quantity and constituting a phase space, it describes any microscopic state of the substance or material, where the coordinates of any particle i are x. i The momentum of particle i is p. i =m i v i m i Let v be the mass of particle i. i Let V be the velocity of particle i, i∈{1,2,…,N}, and the velocities of all particles are V=(v1,v2,…,v…). N ); Assume a three-dimensional Euclidean space The coordinate axes of the Cartesian coordinate system are x i With v i and p i All are 3-dimensional vectors, X, V, and P are all 3N-dimensional vectors; these particles interact with each other, and U(X) represents the potential energy function corresponding to the interaction potential field between all particles. Then, the force on any particle i in the α-axis direction is expressed as: The particle exists in three-dimensional Euclidean space. The coordinates along the α-axis are The particles are atoms or molecules;

[0065] Then proceed to steps 2 and 3.

[0066] Step 2, using a method containing The system of equations for multi-particle system dynamics is used to calculate the changes in the coordinates X and velocity V of all particles in the multi-particle model of the substance and material over time.

[0067] The The equations assign three additional dynamic quantities to each particle in the multi-particle system. The three additional dynamic quantities of each particle form a 3-dimensional vector. The additional kinetic quantities of all particles constitute a 3N-dimensional vector S = (s1, s2, s3, ..., s N ); Set additional dynamic quantities for each particle. The nearest neighbor eccentricity of the particle Corresponding thermodynamic and statistical mechanical quantities;

[0068] The nearest neighbor eccentricity of the particle The calculation method is as follows: first, the nearest neighbor eccentric displacement R(X) of the particle is obtained by calculating the three-dimensional Euclidean space vector from the geometric center of all nearest neighbor particles to the particle's position; then, the absolute values ​​of the three coordinate axis components of the nearest neighbor eccentric displacement R(X) are used. The values ​​of the three components of the nearest neighbor eccentricity d of the particle are respectively used. The nearest-neighbor eccentricities of all particles constitute a 3N-dimensional vector D = (d1, d2, d3, ..., dn). N The nearest neighbor eccentricity d of the particle is a function of the coordinates X of all particles and is used to describe the amplitude of the non-cooperative motion between the particle and its nearest neighbor.

[0069] Due to the nearest neighbor eccentricity d of the particle, the particle's additional kinetic quantity s, and the particle's position in three-dimensional Euclidean space... Position x has the same physical dimensions, and the multi-particle system is considered to be in a 6-dimensional space. Motion in this 6-dimensional space; In this process, each particle is placed in the original three-dimensional Euclidean space. The nearest neighbor eccentricity of the particle as defined in [the document / reference] With the particle in extra space Additional dynamic quantities The coupling is achieved using a harmonic oscillator, and then the potential energy function corresponding to the coupling interaction potential field is superimposed onto the multi-particle system in the original three-dimensional Euclidean space. From the potential energy function U(X) corresponding to the interparticle interaction potential field, we obtain the state of all particles in 6-dimensional space. The coupled superposition potential energy function U corresponding to the interparticle interaction potential field κ (X,S); with the coupled superposition potential function Uκ Based on (X,S), using... The system of equations for multi-particle system dynamics is used to calculate the time-varying coordinates (X) and velocities (V) of all particles in a multi-particle model of the substance or material. The system of equations is constructed such that the particle exists in the original three-dimensional Euclidean space. The dynamic quantity X in the particle changes with time according to the Nose-Hoover dynamic equations, while the particle in the extra space... The variation of the additional kinetic quantity S with time follows the Langevin kinetic equation; the viscosity coefficients of the Nose-Hoover kinetic equation and the Langevin kinetic equation are independent of each other, and their respective bath temperatures are also independent of each other;

[0070] For the contents The dynamic equations of the multi-particle system are subjected to time discretization numerical processing. The coordinates X and velocity V of all particles in the multi-particle model are automatically calculated step by step in the order of time steps to obtain the values ​​of coordinates X and velocity V of all particles in the multi-particle model of matter and materials as a function of time.

[0071] Step 3: Based on the calculated values ​​of the coordinates X and velocity V of all particles corresponding to each discrete time step in the multi-particle model of the substance and material, perform thermodynamic and kinetic analysis on the corresponding substance and material.

[0072] Furthermore, in step 2, the nearest neighbor eccentricity D of the particle is (d1, d2, d3, ..., d N The definition and calculation method of ) are as follows:

[0073]

[0074] in,

[0075] N—The total number of particles in a multi-particle system;

[0076] —respectively, three-dimensional Euclidean space The labels of the three coordinate axes;

[0077] X—All particles in three-dimensional Euclidean space The coordinates in the graph are X = (x1, x2, ..., x...). N ),

[0078] d i —The nearest neighbor eccentricity of any particle i in a multi-particle system is a 3-dimensional vector.

[0079] — These are the nearest neighbor eccentricities d of particle i. i In three-dimensional Euclidean space The three coordinate axes The three components of direction;

[0080] R i (X)—The nearest neighbor eccentric displacement of particle i is a function of the coordinates X of all particles;

[0081] — These represent the nearest neighbor eccentricity displacements R of particle i. i (X) in three-dimensional Euclidean space The three coordinate axes The directional component is a function of the coordinates X of all particles; They represent The corresponding absolute mathematical value;

[0082] N i —i particles with R d The number of all nearest-neighbor particles within the radius;

[0083] —i particles with R d It is the set of all nearest-neighbor particles within the radius;

[0084] x i —At the current moment, particle i is in three-dimensional Euclidean space. Coordinates in;

[0085] x j —At the current moment, particle j is in three-dimensional Euclidean space. The coordinates in the diagram.

[0086] Furthermore, in step 2, a multi-particle system is used in 6-dimensional space. The coupled superposition potential energy function U corresponding to the interparticle interaction potential field κ Based on (X,S), using... The multi-particle system dynamics equations are used to calculate the time-varying coordinates X and velocities V of all particles in a multi-particle model of matter and materials, and the coupled superposition potential energy function U... κ (X,S) and the aforementioned The system of equations has the following expression:

[0087]

[0088] in,

[0089] N—The total number of particles in a multi-particle system;

[0090] —respectively, three-dimensional Euclidean space The labels of the three coordinate axes;

[0091] X—All particles in three-dimensional Euclidean space The coordinates in the graph are X = (x1, x2, ..., x...). N ),

[0092] S—All particles in extra space Additional kinetic quantities, S=(s1,s2,s3,…,s N ),

[0093] —representing the i-particle in three-dimensional Euclidean space The coordinates, velocity, and acceleration along the α-axis in the middle coordinate system.

[0094] —representing the i-particle in the extra space Additional dynamic quantity s i exist The components along the coordinate axes, their first derivative with respect to time, and their second derivative with respect to time.

[0095] ζ—Nose-Hoover dynamics extension variable for multi-particle systems;

[0096] —The first time derivative of the extended variable ζ in the Nose-Hoover dynamics of a multi-particle system;

[0097] m i —i particles in three-dimensional Euclidean space The quality within;

[0098] μ i —i particles in extra space The virtual mass in;

[0099] γ x —Multi-particle systems in three-dimensional Euclidean space The viscosity coefficient corresponding to the Nose-Hoover dynamics;

[0100] T x —Multi-particle systems in three-dimensional Euclidean space The main heat bath temperature corresponding to the Nose-Hoover dynamics;

[0101] k B —Boltzmann constant;

[0102] γ s —Multi-particle systems in extra space The viscosity coefficient corresponding to Langevin dynamics;

[0103] T s —Multi-particle systems in extra space The auxiliary heat bath temperature corresponding to the Langevin dynamics;

[0104] R s (t) — noise function;

[0105] U(X) — A multi-particle system in three-dimensional Euclidean space The potential energy function corresponding to the inter-particle interaction potential field is a function of the coordinates X of all particles;

[0106] —The additional kinetic quantity of the j-particle in Components along the coordinate axes

[0107] —The nearest-neighbor eccentricity of particle j in three-dimensional Euclidean space The components along the α-axis are functions of the coordinates X of all particles.

[0108] κ—The force coupling parameter between the additional dynamic quantity S of a multi-particle system and the nearest neighbor eccentricity D of a particle;

[0109] U κ (X,S)—A multi-particle system in 6-dimensional space The coupled superposition potential energy function corresponding to the inter-particle interaction potential field is a function of the coordinates X of all particles and the additional dynamic quantity S.

[0110] Moreover, in step 2, the content containing The dynamic equations of the multi-particle system are subjected to time discretization numerical processing. The discrete time step is set to Δt, which can be a fixed time step or a non-fixed time step. The automated stepwise calculation solves for the coordinates X and velocity V of all particles in the multi-particle model of matter and materials at each discrete time step, and gives the evolution of the microstate of the matter and materials over time.

[0111] Furthermore, it is used to simulate and analyze the diffusion and migration dynamics of interstitial carbon atoms and iron atom vacancies in a body-centered cubic iron lattice in iron-carbon alloys.

[0112] Alternatively, it can be used to simulate and analyze the kinetics of hydrogen atom diffusion and segregation at grain boundaries in aluminum bicrystals containing a certain concentration of dissolved hydrogen atoms.

[0113] On the other hand, the present invention provides an apparatus for computer simulation of matter and materials, for implementing the improved method for computer simulation of matter and materials as described above.

[0114] The advantage of this invention is that, since the nearest neighbor eccentricity d of the particle proposed and defined in this invention is a 3-dimensional vector... Essentially, it reflects the amplitude of the non-cooperative motion (shuffling, also known as "shuffling motion") between each particle and its nearest neighbor in a multi-particle system. This amplitude of non-cooperative motion typically corresponds to the atomic-molecular scale reaction coordinates of arbitrary internal microstructural transformation processes in condensed matter and materials. Based on this, the present invention cleverly utilizes the nearest-neighbor eccentricity D = (d1, d2, d3, ..., d...) of all particles. N The set of variables (3N-dimensional) for the multi-particle dynamics system is constructed using a 3N-dimensional vector S = (s1, s2, s3, ..., s...). N As an additional dynamic quantity, it is substituted into the extended Langevin dynamic equations proposed by Luca Maragliano and Eric Vanden-Eijnden (see reference

[22] ), and the Nose-Hoover dynamic equations are used as a heat bath method for controlling the temperature of the multi-particle system, thus obtaining the temperature of the multi-particle system. The system of equations describes the evolution of the coordinates X and velocity V of all particles corresponding to their microstate over time. According to the analysis given by Luca Maragliano and Eric Vanden-Eijnden (see reference

[22] ), as long as the parameter μ is set reasonably... i κ, γ x With γ s The value of is then adjusted by regulating the auxiliary heat bath temperature T of the multi-particle system, which corresponds to the additional kinetic quantity S. sBy controlling and enhancing the fluctuations in the non-cooperative motion between each particle (atom or molecule) and its nearest neighbor in condensed matter and materials, accelerated molecular dynamics computer simulations of arbitrary internal microstructural transformation processes in condensed matter and materials can be achieved. Under conditions of a numerical simulation time step Δt comparable to traditional molecular dynamics simulation methods, and within a limited number of computational simulation time steps and computational load comparable to traditional methods, atomic-scale computer numerical simulations can be used to study the physical processes of microstructural evolution in condensed matter and materials with real physical time spans ranging from 1 picosecond to several days. This significantly expands the time scale range of molecular dynamics simulations (by approximately 1 to 13 orders of magnitude) compared to traditional methods. Based on this, the processes and dynamic characteristics of internal microstructural transformations in condensed matter and materials under certain macroscopic thermodynamic conditions can be effectively analyzed. Furthermore, various physical and chemical properties of the studied condensed matter and materials can be analyzed, demonstrating significant application value and broad commercial prospects. Based on the above characteristics, the improved method for general accelerated molecular dynamics computer simulation proposed in this invention is named "Shuffling Accelerated Molecular Dynamics (SAMD)".

[0115] Compared with the patent application "Computer Simulation Method and Apparatus for Matter and Materials Across Time Scales" submitted by the applicant's research team in July 2021, the present invention, in step 2, calculates the evolution of the coordinates X and velocity V of all particles in the multi-particle model of matter and materials over time using a method that includes... The system of equations for the dynamics of a multi-particle system, where for particles in the original three-dimensional Euclidean space... The dynamic quantity X in the previous example was calculated using the Nose-Hoover dynamic equations, while the method proposed in the application titled "Computer Simulation Method and Device for Matter and Materials Across Time Scales" uses the four-dimensional Langevin dynamic equations for multi-particle systems, where the particle's behavior in the original three-dimensional Euclidean space... The change of the kinetic quantity X over time is calculated using the Langevin kinetic equation. Since the Nose-Hoover kinetic equation is closer to real Newtonian kinetics in describing classical multi-particle dynamics than the Langevin equation, it has higher kinetic reliability and fidelity. Furthermore, the Nose-Hoover kinetic equation can control the temperature of the multi-particle system more accurately and effectively than the Langevin equation. Therefore, in practical implementation, the method proposed in the application titled "Method and Device for Computer Simulation of Matter and Materials Across Time Scales" requires a series of tests and runs to continuously adjust the value of the main heat bath temperature parameter in the four-dimensional extended Langevin kinetic equation set of the multi-particle system to determine a suitable value for the parameter in order to accurately control the temperature of the simulated multi-particle model of matter and materials. The method of this invention eliminates this cumbersome process, directly setting the value of the main heat bath temperature parameter in the Nose-Hoover kinetic equation set to the target temperature of the simulated multi-particle model of matter and materials, greatly improving the convenience and ease of use in implementation and application.

[0116] Furthermore, the nearest neighbor eccentricity d set in step 2 of this invention to describe the amplitude of non-cooperative motion between any particle and its nearest neighbor particles has the same function as the particle nearest neighbor eccentricity metric D set in the method of the patent application entitled "Method and Apparatus for Computer Simulation of Matter and Materials Across Time Scales". However, the two have significantly different calculation methods and data dimension characteristics. In the method of the patent application entitled "Method and Apparatus for Computer Simulation of Matter and Materials Across Time Scales", the particle nearest neighbor eccentricity metric D is calculated by averaging the sum of the three coordinate axis components of the nearest neighbor eccentricity displacement R(X) obtained by taking a certain particle as the starting point and then the sum of the three-dimensional Euclidean space vectors to all its nearest neighbor particles. The particle nearest neighbor eccentricity metric D obtained by this calculation method is a scalar (1-dimensional) and uses a vector S = (s1, s2, s3, ..., s2) with dimension N. N The method of coupling the nearest neighbor eccentricity measure of all particles with an additional dynamic quantity is a relatively general averaging treatment of the geometric characteristics of the particle's nearest neighbor eccentricity. However, this invention uses the absolute values ​​of the three coordinate axis components of the nearest neighbor eccentricity displacement R(X) to calculate the values ​​of the three components of the particle's nearest neighbor eccentricity d. The nearest neighbor eccentricity d of the particle defined and calculated in this way is a 3-dimensional vector. Therefore, it provides a more comprehensive description of the geometric characteristics of the particle's nearest neighbor eccentricity, and correspondingly adopts a 3N-dimensional vector S = (s1, s2, s3, ..., s N ) as an additional kinetic quantity (where The model is coupled with the nearest neighbor eccentricity of all particles, so that the absolute values ​​of the three coordinate axis components of the nearest neighbor eccentricity displacement R(X) of each particle can be extended and dynamically coupled. Therefore, it can more accurately and effectively describe and reflect the reaction coordinates of the microstructure transformation process in the multi-particle model, and obtain computer simulation data results with higher fidelity when performing computer simulation on actual materials.

[0117] It should be noted that the new calculation method proposed in this invention is a new improvement approach obtained by the applicant's research team from another perspective and with a different mindset, after extensive exploration, research and practice.

[0118] The present invention is simple and convenient to implement, highly practical, and highly reliable. It solves the problems of low practicality, low fidelity, and inconvenience in practical application of related technologies, improves user experience, and has significant market value. Attached Figure Description

[0119] Figure 1 This is a schematic diagram of the main execution steps of the SAMD simulation method proposed in the embodiments of the present invention.

[0120] Figure 2 This is a schematic diagram of the atomic model of an iron crystal containing one iron atom vacancy and three interstitial carbon atoms, as discussed in Example 1 of the specific application of the embodiment of the present invention. Sub-figure (A) shows the lattice of the iron crystal with a body-centered cubic structure, and sub-figure (B) shows a schematic diagram of the iron atom vacancy and interstitial carbon atoms in the iron crystal.

[0121] Figure 3 The display is based on Figure 2 (B) Model, given the actual thermodynamic temperature of the system as 300K, shows the change of the root mean square displacement of carbon atoms with the number of simulation time steps, as given by computer numerical simulation using the Conventional Molecular Dynamics (MD) method and the SAMD simulation method proposed in this invention.

[0122] Figure 4 The display is based on Figure 2 Model (B) employs the SAMD simulation method proposed in this invention, with the auxiliary heat bath temperature of the system set to T. s Computer numerical simulations were performed at 50000K, and the resulting model was... Figure 3 The middle arrows indicate different time steps (n). steps The results of the positions and migration of interstitial carbon and iron vacancies at the corresponding time: (A)n steps =0; (B)n steps =5e+05; (C)n steps=1e+06; (D)n steps =1.5e+06; (E)n steps =2e+06.

[0123] Figure 5 This is a schematic diagram of the atomic model of a metallic aluminum bicrystal with a certain concentration of dissolved hydrogen atoms, as discussed in Example 2 of the specific application of this invention. Sub-figure (A) shows a projection of the aluminum bicrystal model with a certain concentration of dissolved hydrogen atoms (the large sphere represents aluminum atoms, the small sphere represents hydrogen atoms, and the distribution of hydrogen atoms has not reached equilibrium). Sub-figure (B) shows the case where all aluminum atoms are hidden and only the initial distribution of hydrogen atoms is shown.

[0124] Figure 6 The display is based on Figure 5 The model, given the actual thermodynamic temperature of the system as 100K, shows the change of the root mean square displacement of hydrogen atoms with the number of simulation time steps, obtained by computer numerical simulation using both the Conventional Molecular Dynamics (MD) method and the SAMD simulation method proposed in this invention.

[0125] Figure 7 The display is based on Figure 5 The model, using the SAMD simulation method proposed in this invention, sets the auxiliary heat bath temperature of the system to T. s Computer numerical simulations were performed at 500K, and the resulting model was... Figure 6 The middle arrows indicate different time steps (n). steps The results of the distribution and migration of hydrogen atoms at the corresponding time point: (A)n steps =0; (B)n steps =5e+05; (C)n steps =1e+06; (D)n steps =1.5e+06; (E)n steps =2e+06.

[0126] Figure 8 The display is based on Figure 5 The model, using the SAMD simulation method proposed in this invention, sets the auxiliary heat bath temperature T of the system. s Computer numerical simulations were performed at 100K (subgraph (A)), 300K (subgraph (B)), and 500K (subgraph (C)) respectively, with a total time step of 2.0 × 10⁻⁶. 6 The simulation results show the final hydrogen atom distribution under the given simulation time length; subplot (D) shows the hydrogen atom distribution obtained by using the Monte Carlo (MC) computer numerical simulation method, also at the actual thermodynamic temperature of the model system of 100K, reflecting the equilibrium distribution state of hydrogen atoms at this thermodynamic temperature.

[0127] Figure 9 This is a flowchart of a computer simulation method for materials in the prior art. Detailed Implementation

[0128] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings and embodiments.

[0129] This invention employs a physical quantity, which can be called the "nearest neighbor eccentricity" of a particle, as the "collective variable" of multi-particle system dynamics. It defines the corresponding thermodynamic and statistical mechanical quantities as extended dynamic variables of the multi-particle system. Based on this, it proposes an improved method for general accelerated molecular dynamics computational simulation, which can be named "Shuffling Accelerated Molecular Dynamics (SAMD)" (hereinafter referred to as the SAMD simulation method). By numerically solving the Nose-Hoover⊕Langevin dynamic equations of the multi-particle system (to be presented below) on a computer, it simulates the evolution of the microscopic state of condensed matter and materials over time, and analyzes the process and dynamic characteristics and laws of internal microstructural transformations of condensed matter and materials under certain macroscopic thermodynamic conditions.

[0130] It is important to note that all uses of the term "particle" in this article refer to the small, physical, and chemically small parts of an object with volume and mass. Depending on their size, this range can encompass subatomic particles such as electrons, protons, and neutrons, microscopic particles such as atoms and molecules, mesoscopic particles such as biomolecules and colloidal particles, and even macroscopic particles such as powders and other granular substances. The "particles" used in the description of the invention's techniques and methods in this article primarily refer to the microscopic particles that constitute matter and materials, namely atoms or molecules. One of their main characteristics is that they can be described using the concept of point masses in physical theory. This does not preclude the application of the invention's techniques and methods to various other particles, such as biomolecules, colloidal particles, and other granular substances, as research objects; therefore, their use in these contexts is also protected by this invention patent.

[0131] See Figure 1 The SAMD simulation method proposed in this embodiment of the invention mainly includes the following steps:

[0132] Step 1: For the substances and materials under study and requiring computer simulation, construct a multi-particle model and system for those substances and materials. The implementation method for this step can be found in the patent applications previously filed by the inventors' research team. For ease of implementation and reference, the implementation instructions are provided below:

[0133] For any condensed matter or material containing N particles (atoms or molecules), the set of all its particles is denoted as {1,2,…,N}. The arrangement of all particles in three-dimensional Euclidean space... (using Cartesian coordinate system) express, For three-dimensional Euclidean space The coordinates X = (x1, x2, ..., x) in the three coordinate axes) N ) and the corresponding particle momentum P=(p1,p2,…,p N As a dynamic quantity and constituting a phase space, it describes any microscopic state of a multi-particle system. The coordinates of any particle i are... The momentum of particle i is p i =m i v i m i Let v be the mass of particle i. i Let be the velocity of particle i, i∈{1,2,…,N}. The velocities of all particles are V = (v1, v2, ..., v...). N ), x i With v i and p i All are 3-dimensional vectors, and X, V, and P are all 3N-dimensional vectors. These particles interact with each other. Let U(X) represent the potential energy function of the interactions between all particles in the system. Then, the force on any particle i along the α-axis can be expressed as... The particle exists in three-dimensional Euclidean space. The coordinates along the α-axis are

[0134] The construction of the multi-particle model of this matter and material requires the values ​​of N and M (where M = (m1, m2, ..., m...). N The data for the U(X) function and the initial values ​​of X and V are given. In practice, a multi-particle model of matter and materials is constructed by providing the values ​​of N and M (N and M are generally constants), the data for the U(X) function, and the initial values ​​of X and V. The next task is to determine how the values ​​of X and V of this multi-particle model change with time under the influence of the inter-particle interaction potential field with potential energy function U(X).

[0135] The technical features that distinguish this invention from the prior art are described in subsequent steps 2 and 3.

[0136] Step 2: Studies of condensed matter and material systems typically focus on the dynamics of various internal microstructural transformations, rather than the random thermal vibrations of each particle near its equilibrium position. This invention proposes using the nearest-neighbor eccentricity *d*, defined by the following formula, to describe the amplitude of the non-cooperative motion (shuffling, also known as "walking motion") between each particle and its nearest neighbor in a multi-particle model:

[0137]

[0138] in,

[0139] N—The total number of particles in a multi-particle system;

[0140] —respectively, three-dimensional Euclidean space The labels of the three coordinate axes;

[0141] X—All particles in three-dimensional Euclidean space The coordinates in the graph are X = (x1, x2, ..., x...). N ),

[0142] d i —The nearest neighbor eccentricity of any particle i is a 3-dimensional vector;

[0143] — These are the nearest neighbor eccentricities d of particle i. i In three-dimensional Euclidean space The three coordinate axes The three components of direction;

[0144] R i (X)—The nearest neighbor eccentric displacement of particle i is a function of the coordinates X of all particles;

[0145] — These represent the nearest neighbor eccentricity displacements R of particle i. i (X) in three-dimensional Euclidean space The three coordinate axes The directional component is a function of the coordinates X of all particles; They represent The corresponding absolute mathematical value;

[0146] N i —i particles with R d The number of all nearest-neighbor particles within the radius;

[0147] —i particles with R d It is the set of all nearest-neighbor particles within the radius;

[0148] x i —At the current moment, particle i is in three-dimensional Euclidean space. Coordinates in;

[0149] x j —At the current moment, particle j is in three-dimensional Euclidean space. The coordinates in the diagram.

[0150] The nearest-neighbor eccentricities of all particles form a 3N-dimensional vector D = (d1, d2, d3, ..., dn). N The amplitude of the non-cooperative motion of particles, given by this formula, usually corresponds to the reaction coordinates of any internal microstructural transformation process in condensed matter and materials. Therefore, the nearest neighbor eccentricity D = (d1, d2, d3, ..., d) of all particles is defined by this formula. N ) can be used to describe the reaction coordinates of arbitrary internal microstructural transformation processes in condensed matter and materials.

[0151] From the definition of the nearest neighbor eccentricity d of the particle above, it can be seen that as long as the parameter R is given... d The value of the nearest neighbor eccentricity d of each particle in a multi-particle system at any given time can be directly calculated based on the coordinates X of all particles at the current time using the above expression.

[0152] Based on the definition of the nearest neighbor eccentricity d of the aforementioned particles, this invention proposes that the values ​​of the coordinates X and velocity V of all particles in a multi-particle model change with time under the influence of the inter-particle interaction potential field with potential energy function U(X). The proposed approach is to use the following form... The calculations are performed using the system of dynamic equations:

[0153]

[0154] in,

[0155] N—The total number of particles in a multi-particle system;

[0156] —respectively, three-dimensional Euclidean space The labels of the three coordinate axes;

[0157] X—All particles in three-dimensional Euclidean space The coordinates in the graph are X = (x1, x2, ..., x...). N ),

[0158] S—All particles in extra space Additional kinetic quantities, S=(s1,s2,s3,…,s N ),

[0159] —representing the i-particle in three-dimensional Euclidean space The coordinates, velocity, and acceleration along the α-axis in the middle coordinate system.

[0160] —representing the i-particle in the extra space Additional dynamic quantity s i exist The components along the coordinate axes, their first derivative with respect to time, and their second derivative with respect to time.

[0161] ζ—Nose-Hoover dynamics extension variable for multi-particle systems;

[0162] —The first time derivative of the extended variable ζ in the Nose-Hoover dynamics of a multi-particle system;

[0163] m i —i particles in three-dimensional Euclidean space The quality within;

[0164] μ i —i particles in extra space The virtual mass in;

[0165] γ x —Multi-particle systems in three-dimensional Euclidean space The viscosity coefficient corresponding to the Nose-Hoover dynamics;

[0166] T x —Multi-particle systems in three-dimensional Euclidean space The main heat bath temperature corresponding to the Nose-Hoover kinetics; k B —Boltzmann constant;

[0167] γ s —Multi-particle systems in extra space The viscosity coefficient corresponding to Langevin dynamics;

[0168] T s —Multi-particle systems in extra space The auxiliary heat bath temperature corresponding to the Langevin dynamics;

[0169] R s (t) — Noise function, such as white noise function;

[0170] U(X) — A multi-particle system in three-dimensional Euclidean space The potential energy function corresponding to the inter-particle interaction potential field is a function of the coordinates X of all particles;

[0171] —The additional kinetic quantity of the j-particle in Components along the coordinate axes

[0172] —The nearest-neighbor eccentricity of particle j in three-dimensional Euclidean space The components along the α-axis are functions of the coordinates X of all particles.

[0173] κ—The force coupling parameter between the additional dynamic quantity S of a multi-particle system and the nearest neighbor eccentricity D of a particle;

[0174] U κ (X,S)—A multi-particle system in 6-dimensional space The coupled superposition potential energy function corresponding to the inter-particle interaction potential field is a function of the coordinates X of all particles and the additional dynamic quantity S.

[0175] Here, the present invention endows each particle in a multi-particle system with three additional kinetic quantities. These three additional dynamic quantities constitute a vector. The additional kinetic quantities s of all particles constitute a 3N-dimensional vector S = (s1, s2, s3, ..., s N The nearest-neighbor eccentricity D of S and all particles in the multi-particle system is (d1, d2, d3, ..., d...). N The physical meaning of this correspondence is that, if we take the nearest neighbor eccentricity D = (d1, d2, d3, ..., d...) of all particles as the basis, then... N As a set of variables (3N-dimensional) in a multi-particle dynamics system, the 3N-dimensional vector S = (s1, s2, s3, ..., s...) N ) are the thermodynamic and statistical mechanical quantities corresponding to the variables in this set.

[0176] The physical meaning of this system of equations can be understood as follows: the particle's additional dynamic quantity s corresponds to the particle's nearest neighbor eccentricity d, and according to the definition of the particle's nearest neighbor eccentricity d, d is related to the particle in three-dimensional Euclidean space. The coordinates have the same physical dimensions, therefore the dynamic quantity s is also the same as that of the particle in three-dimensional Euclidean space. The coordinates have the same physical dimensions, while each particle is given an additional dynamic quantity. This is equivalent to adding 3 extra dynamic degrees of freedom to each particle, so the multi-particle system can be considered to exist in a 6-dimensional space. Motion in this 6-dimensional space In the middle, all particles are in the original three-dimensional Euclidean space The particle's nearest neighbor eccentricity D and the particle's additional kinetic quantity S are coupled in a harmonic oscillator manner, and the potential energy function corresponding to this coupled interaction potential field is then superimposed onto the multi-particle system in the original three-dimensional Euclidean space. From the potential energy function U(X) corresponding to the interparticle interaction potential field, we obtain the state of all particles in 6-dimensional space. The coupled superposition potential energy function U corresponding to the interparticle interaction potential field κ (X,S), under the coupled superposition potential function U κ Under the influence of the interparticle interaction potential field corresponding to (X,S), the dynamic quantities X and S of the multi-particle system evolve over time according to the Nose-Hoover dynamic equation and the Langevin dynamic equation, respectively, with independent viscosity coefficients and bath temperature. Therefore, this invention refers to the above set of equations as the multi-particle system in 6-dimensional space. The system of dynamic equations.

[0177] The main task of the SAMD simulation method is to simulate the above-mentioned... The dynamic equations were solved numerically on a computer. It can be seen that this system of equations is a coupled second-order ordinary differential equation system with time as the independent variable. This can be solved through simple variable substitution (such as...). This can be transformed into a system of coupled first-order ordinary differential equations, then the time can be numerically discretized, with the discrete time step set to a fixed Δt (or a non-fixed time step), and the differential can be replaced with a difference scheme (e.g., This is transformed into a system of ordinary algebraic equations, based on the initial values ​​of (X,V), (N,M), and the data of the U(X) function given in step 1, as well as the calculation formula for the particle's nearest neighbor eccentricity d and the parameter R. d The value is then given. The parameters μ in the system of equations i κ, γ x γ s T x T s The value of can be obtained step by step to solve for the coordinates X and velocity V of all particles corresponding to each discrete time step (Δt, 2Δt, 3Δt, ..., nΔt), giving the evolution of the microscopic state of condensed matter and materials over time.

[0178] The specific computational solution to this system of equations can be achieved through different implementations and numerical processing details, such as the Verlet algorithm, the Leap frog algorithm, the Gear prediction-correction algorithm, etc. (see references [1,26]). These specific numerical computational processing techniques are already quite mature and will not be elaborated further. It should be noted that, although the above... Different numerical solutions and algorithms can be used to solve the kinetic equations, but as long as a calculation method incorporating these equations is employed and the computer numerical simulation of matter and materials is performed according to the steps outlined in this specification, its use is protected by this invention patent. Step 3: Based on the values ​​of (X,V) corresponding to each discrete time step (Δt, 2Δt, 3Δt, ..., nΔt) of the multi-particle model of condensed matter and materials calculated in Step 2, various thermodynamic and kinetic analyses are performed on the corresponding condensed matter and materials, including the spatial configuration and temporal evolution analysis of the multi-particle system, the process of internal microstructural transformation of the system under certain macroscopic thermodynamic conditions, and its dynamic characteristics and laws, etc. Since these analytical methods, means, and techniques can refer to the methods and means used in traditional ordinary molecular dynamics simulation methods, this invention will not elaborate further.

[0179] Figure 9 The method flowchart from patent application No. 202110778355.3, entitled "Method and Apparatus for Computer Simulation of Matter and Materials Across Time Scales," is provided. The differences in the flowchart between the method of this invention and the method of the cited patent are clearly visible. Specifically, the method of this invention employs a multi-particle system in step 2. The dynamic equations are used to calculate the changes in the coordinates X and velocity V of all particles in a multi-particle model of matter and materials over time. Figure 1 The method in the patent with patent number 202110778355.3 uses the four-dimensional Langevin dynamic equations of a multi-particle system for related calculations. The new calculation method proposed in this invention is a new improvement approach obtained by the applicant's research team from another perspective and idea after a lot of exploration, research and practice. In specific implementation, it eliminates the tedious process of the patent method with patent number 202110778355.3, which requires a series of tests and runs in advance to continuously adjust the value of the main heat bath temperature parameter in the four-dimensional extended Langevin dynamic equations of the multi-particle system to determine the appropriate value of the parameter.

[0180] For ease of implementation and reference, two specific application examples are provided below. First, follow the basic framework of steps 1-2-3 listed in the above implementation method (see appendix). Figure 1 The SAMD simulation method (mainly in step 2) The dynamic equations are numerically discretized and processed using a computer-automated computation process. Then, for any substance or material of interest to the user, simply prepare a multi-particle model of that substance or material according to step 1, including the values ​​of (N,M), the data of the interparticle (atom or molecule) interaction potential energy function U(X), and the initial values ​​of (X,V) corresponding to the microstate of the system, using these as input data. Simultaneously, the numerical values ​​of various parameters required for the SAMD simulation method are provided, mainly including:

[0181] 1) The time step Δt used for time discretization of the dynamic equations (the recommended range is 0.1 to 2.0 fs);

[0182] 2) The cutoff radius R required to construct the nearest neighbor set for each particle in order to calculate the nearest neighbor eccentricity D of the particle. d (Generally, it can be set to the position corresponding to the first valley of the radial distribution function curve of the multi-particle system. Its value should be such that the number of particles in the nearest neighbor set of each particle in the multi-particle system is between 8 and 20).

[0183] 3) The parameter μ of the dynamic equation system i (Generally, m can be taken as the actual mass of the corresponding particle) i (value);

[0184] 4) The parameters κ and γ of the dynamic equations x γ s (The recommended value ranges are as follows:) 10-200ps -1 1~10ps -1 );

[0185] 5) The parameter T of the dynamic equations x (The main heat bath temperature of the multi-particle system is set according to the simulation requirements and actual conditions.)

[0186] 6) The parameter T of the dynamic equations s (The auxiliary heat bath temperature of the multi-particle system, which corresponds to the additional dynamic quantity S of the particles, is set according to the simulation requirements and actual conditions. It can generally be taken as a value in the range of 100 to 100,000 K.)

[0187] 7) Total running time steps N Run (The value is set according to the simulation requirements and actual situation, and is generally taken as 100 to 1×10.) 8 (The range varies).

[0188] Then, a computer-automated computation process is used to solve for each discrete time step (Δt, 2Δt, 3Δt, ..., N). Run The value of (X,V) corresponding to Δt can give the evolution of the microscopic state of the substance and material over time, enabling computer simulation of the dynamic process of the multi-particle system of the substance and material. Finally, based on this, the microstructural transformation events and corresponding dynamic characteristics and laws occurring inside the substance and material system are analyzed.

[0189] In the process of implementing the specific application example, the most popular and widely used molecular dynamics simulation software LAMMPS (see reference [2]) was used. With the open source nature, robust code quality, rich functions, and its excellent modular programming framework and convenient extension interface of LAMMPS, the computer automated operation process of the SAMD simulation method can be easily compiled into computer program code and integrated into the LAMMPS software. In addition, two new LAMMPS input script commands (compute nnoad / atom and fix samd / nnoad / atom) specifically designed for the SAMD simulation method are provided for the LAMMPS software. When using this modified and expanded LAMMPS software, the source code is compiled on the computer to obtain an executable file. Then, the simulation task is performed on the computer according to the usual steps of using LAMMPS software for ordinary molecular dynamics simulation (see LAMMPS User Manual, Reference [2]): setting up and compiling the LAMMPS input command line script - providing various necessary input data files - submitting the LAMMPS executable program to perform calculations - analyzing the output result data files. It is only necessary to add and set the two new command lines (compute nnoad / atom and fix samd / nnoad / atom) specifically for the SAMD simulation method in the LAMMPS input command line script file to perform SAMD simulation.

[0190] It should be emphasized that although the LAMMPS open-source software was used in the specific application examples below, the computer-automated calculation process of the SAMD simulation method was compiled into computer program code and integrated into the LAMMPS software, and then the SAMD simulation method proposed in this invention was implemented using this modified and functionally expanded LAMMPS software, the implementation of this invention is not limited by the LAMMPS software in any way. Any other feasible computer programming scheme and corresponding implementation method can be used to specifically implement the technical solution of this invention, and all implementations and uses of this invention using these computer programming schemes are protected by this invention patent.

[0191] Specific application example 1: The diffusion and migration dynamics of interstitial carbon and iron vacancies in a body-centered cubic iron lattice in iron-carbon alloys.

[0192] The most basic component of commonly used steel materials is an iron-carbon alloy, consisting of iron atoms arranged periodically in a body-centered cubic lattice, while a small number of carbon atoms reside as interstitial atoms in the octahedral interstices of the iron lattice, as shown in the attached diagram. Figure 2 As shown in (A). At normal temperatures (such as room temperature), iron atoms undergo thermal vibrations at their lattice sites, and have a certain probability of detaching from their sites to become interstitial atoms, while leaving an iron atom vacancy at one lattice site, as shown in the attached diagram. Figure 2 As shown in (B), the diffusion and migration of carbon atoms and the diffusion and migration of iron atoms through vacancies are important kinetic processes in iron-carbon alloys. They play a crucial role in material processing such as carburizing of steel surfaces and heat treatment of steel that cause changes in microstructure and mechanical properties. They are key physicochemical processes that affect these processes and material properties.

[0193] Step 1: For the attached Figure 2 The atomic model of iron crystal (B) (containing 1 iron atom vacancy and 3 interstitial carbon atoms) was constructed and the data preparation was carried out as follows:

[0194] 1) The interatomic interaction potential energy function U(X) of the Fe-C binary alloy system is provided by a special potential function file (which can be obtained from the publicly available potential function file library on the Internet, see reference

[27] , and here the file is named "potential.file"). The potential function file "potential.file" contains the mass data of Fe atoms and C atoms respectively.

[0195] 2) Set up and write a LAMMPS input command script file. In this file, add commands such as lattice, region, create_box, create_atoms, and delete_atoms to fill the model space with iron atoms (a total of 5400 iron atoms). Also, set up to add 3 carbon atoms at different octahedral interstices in the iron crystal model. At the same time, specify to delete 1 iron atom in the model so that the position of the atom becomes a vacancy. Use the pair_style and pair_coeff commands to set the interatomic interaction potential function file "potential.file" of the Fe-C binary alloy system in the model. Then use commands such as minimize, timestep, velocity, fix npt, and run to specify the relaxation of the initial model. Use the write_restart command to set the final output results to the file "model.file".

[0196] 3) Once the input command script file is generated, it can be passed to the LAMMPS executable program and submitted to the computing server for computation to obtain the corresponding output file "model.file". This file contains the value of the total number of atoms N in the model and the initial values ​​of the coordinates X and velocity V of all atoms.

[0197] Step 2: Based on the model state corresponding to the output file "model.file" obtained in Step 1, perform SAMD simulation on the model using LAMMPS software. The specific implementation process is as follows:

[0198] 1) Set up and write a new LAMMPS input command script file. In this file, add the `read_restart` command to read the output file "model.file" obtained in step 1 to import the model. Use the `pair_style` and `pair_coeff` commands to specify the interatomic interaction potential function file "potential.file" (same as the file in step 1) for the Fe-C binary alloy system. Use commands such as `timestep`, `compute nnoad / atom`, `fix nvt`, and `fix samd / nnoad / atom` to specify the values ​​of various parameters for SAMD simulation of the model, as shown in the table below:

[0199]

[0200] Here, the cutoff radius R d value This is exactly the position of the first valley of the radial distribution function of the body-centered cubic iron atom lattice, meaning each iron atom normally has 14 nearest-neighbor atoms. The `compute msd` command is used to specify the calculation of the root mean square displacement of carbon atoms during the simulation. The `fix nph` command is used to specify the use of Nose-Hoover dynamics to control the model's pressure to a zero-stress state. The `thermo` command is used to specify outputting the model's thermodynamic state information (including the aforementioned root mean square displacement of carbon atoms) to the file "thermo.file" every 5 time steps. The `restart` command is used to specify outputting various dynamic state information of the model at the current time step (including the coordinates X and velocity V of all atoms at the current moment) to the file "kinetics.file" every 10000 time steps for subsequent analysis of the model's spatial configuration evolution over time. The `run` command specifies the total number of running time steps N for the SAMD simulation. Run It is 2,000,000.

[0201] 2) Once the input command script file is generated, it can be passed to the LAMMPS executable program and submitted to the computing server for calculation, resulting in the corresponding SAMD simulation output data files "thermo.file" and "kinetics.file".

[0202] Step 3: Analyze the output data (thermodynamic state information and kinetic state information) files "thermo.file" and "kinetics.file" obtained in Step 2. Extract the root mean square displacement of carbon atoms and the corresponding time step number from the thermodynamic state information output file "thermo.file". Use plotting software to show the change of the root mean square displacement of carbon atoms with the increase of time step number (see appendix). Figure 3 Based on the spatial configuration data of the model given in the kinetics.file, the OVITO atomic configuration visualization program was used to perform computer visualization image processing and analysis on the model's spatial configuration, and to derive the model at different time steps (n). steps Images of the atomic spatial configurations (including carbon and iron vacancies) under (1)n steps =0 (attached) Figure 4 (A)); (2)n steps =5e+05 (attached) Figure 4 (B));(3)n steps =1e+06 (Appendix) Figure 4 (C)); (4)n steps =1.5e+06 (Appendix) Figure 4 (D));(5)n steps =2e+06 (Appendix) Figure 4 (E)).

[0203] In steps 1 and 2 of this specific application example, apart from the newly added compute nnoad / atom and fix samd / nnoad / atom commands specifically designed for the SAMD simulation method, the use and settings of all LAMMPS commands involved are performed in the conventional manner according to the instructions in the LAMMPS user manual (see reference [2]). Some related details will not be elaborated here. Detailed instructions for these commands can be found in the LAMMPS user manual (see reference [2]), and they will not be further explained here.

[0204] It should be noted that the interatomic interaction potential function for the Fe-C binary alloy system in this model is the embedded atom method-type interatomic interaction potential developed by CSBecquart et al. for the Fe-C binary system (see references [28,29]). Verification shows that this potential function can realistically reflect the basic physical properties of the Fe-C alloy and has good reliability and rationality. Calculations show that the activation energy barriers for vacancy diffusion of interstitial carbon atoms and iron atoms given by this potential function are 0.81 eV / atom and 0.63 eV / atom, respectively, which are close to the experimental results.

[0205] From the appendix Figure 3 As can be seen from the above, the SAMD simulation method proposed in this invention is used to simulate the attached... Figure 2 (B) The model was simulated using computer numerical simulation. When the auxiliary heat bath temperature Ts < 10000K, the total time step was 2.0 × 10⁻⁶. 6 Within the simulation timeframe, no migration of carbon or iron atom vacancies was observed; from the attached Figure 3 and attached Figure 4 As can be seen from (A)-(E), as the secondary heat bath temperature Ts of the system increases from 20000K to 50000K, the migration frequency and migration distance of carbon atoms gradually increase. This indicates that by adjusting the value of the secondary heat bath temperature Ts parameter, the extent to which the SAMD simulation method proposed in this invention expands and spans across time scales can be effectively controlled.

[0206] Simple estimations show that for SAMD simulations of the system's auxiliary heat bath temperatures Ts of 20000K, 30000K, and 50000K, with a total time step of 2.0 × 10⁻⁶, the results are satisfactory. 6 Within the simulation timeframe, the average number of migration jumps for carbon atoms in the model was approximately 4.81, 62.11, and 75.32, respectively. Experimental calculations and analysis show that at a thermodynamic temperature of 300 K, carbon atoms jump an average of approximately 2 times per second between adjacent octahedral interstitial positions. Comparative analysis reveals that in this specific application example, for SAMD simulations with secondary heat bath temperatures Ts of 20000 K, 30000 K, and 50000 K, the corresponding time speedup ratios are approximately 1.2 × 10⁻⁶. 9 1.5×10 10 2.1×10 10 .

[0207] Specific application example 2: The kinetics of hydrogen atom diffusion and segregation at grain boundaries in a metallic aluminum bicrystal containing a certain concentration of dissolved hydrogen atoms.

[0208] Various hydrogen-containing substances exist in nature and engineering environments, such as widely distributed water, as well as various acid, alkali, and salt solutions (like seawater). When metal structural components come into contact with these hydrogen-containing substances during service, hydrogen atoms can penetrate into the metal due to surface chemical or electrochemical reactions, weakening the interatomic bonding strength and causing what is known as "hydrogen embrittlement." This poses a serious safety hazard to the service life of metal structural components and the use of engineering equipment. Because hydrogen atoms are the smallest atoms, their diffusion in metallic materials is generally relatively easy. Furthermore, experimental studies show that hydrogen atoms dissolved in metals tend to segregate at crystal defects such as dislocations and grain boundaries. Segregation at grain boundaries, in particular, often leads to significant weakening of the grain boundaries, becoming the source of crack initiation and the preferred path for crack propagation. Therefore, a deep understanding of the diffusion and segregation behavior of hydrogen atoms in aluminum crystals is of significant reference value for studying the hydrogen embrittlement behavior of metallic aluminum.

[0209] Step 1: For the attached Figure 5 The image shows a metallic aluminum bicrystal containing a certain concentration of dissolved hydrogen atoms (with grain boundaries of Σ5). <110> The atomic model of the {310} symmetric tilted grain boundary was first constructed and the data was prepared using LAMMPS software as follows:

[0210] 1) The interatomic interaction potential energy function U(X) of the Al-H binary alloy system is provided by a special potential function file (which can be obtained from the publicly available potential function file library on the Internet, see reference

[27] , and here the file is named "potential.file"). The potential function file "potential.file" contains the mass data of Al atoms and H atoms respectively.

[0211] 2) Set up and write a LAMMPS input command script file. In this file, add commands such as lattice, region, create_box, create_atoms, delete_atoms to specify that aluminum atoms and hydrogen atoms (a total of 5664 aluminum atoms and 200 hydrogen atoms) are filled in the model space according to the lattice type and grain orientation. Also, specify that one atom in a pair of atoms that are too close to each other should be deleted. Use the pair_style and pair_coeff commands to specify the interatomic interaction potential function file "potential.file" for the Al-H binary alloy system in the model. Then use commands such as minimize, timestep, velocity, fix npt, run to specify the relaxation of the initial model. Use the write_restart command to specify the final output results to the file "model.file".

[0212] 3) Once the input command script file is generated, it can be passed to the LAMMPS executable program and submitted to the computing server for computation, resulting in the corresponding output file "model.file". This file contains the value of the total number of atoms N in the model and the initial values ​​of the coordinates X and velocity V of all atoms. Finally, the model contains a total of 5846 atoms (5664 aluminum atoms and 182 hydrogen atoms).

[0213] Step 2: Based on the model state corresponding to the output file "model.file" obtained in Step 1, perform SAMD simulation on the model using LAMMPS software. The specific implementation process is as follows:

[0214] 1) Set up and write a new LAMMPS input command script file. In this file, add the `read_restart` command to read the output file "model.file" obtained in step 1 to import the model. Use the `pair_style` and `pair_coeff` commands to specify the interatomic interaction potential function file "potential.file" (same as the file in step 1) for the Al-H binary alloy system. Use commands such as `timestep`, `compute nnoad / atom`, `fix nvt`, and `fix samd / nnoad / atom` to specify the values ​​of various parameters for SAMD simulation of the model, as shown in the table below:

[0215]

[0216] Here, the cutoff radius R d value To determine the location of the first valley of the radial distribution function of the face-centered cubic aluminum atomic lattice, the `compute msd` command is used to specify the calculation of the root mean square displacement of hydrogen atoms during the simulation. The `fix nph` command is used to specify the use of Nose-Hoover dynamics to control the model's pressure to a zero-stress state. The `thermo` command is used to specify outputting the model's thermodynamic state information (including the aforementioned root mean square displacement of hydrogen atoms) to the file "thermo.file" every 5 time steps. The `restart` command is used to specify outputting various dynamic state information of the model at the current time step (including the coordinates X and velocities V of all atoms at the current moment) to the file "kinetics.file" every 10,000 time steps for subsequent analysis of the model's spatial configuration evolution over time. The `run` command is used to specify the total number of running time steps N for the SAMD simulation. Run It is 2,000,000.

[0217] 2) Once the input command script file is generated, it can be passed to the LAMMPS executable program and submitted to the computing server for calculation, resulting in the corresponding SAMD simulation output data files "thermo.file" and "kinetics.file".

[0218] Step 3: Analyze the output data (thermodynamic state information and kinetic state information) files "thermo.file" and "kinetics.file" obtained in Step 2. Extract the root mean square displacement of hydrogen atoms and the corresponding time step number from the thermodynamic state information output file "thermo.file". Use plotting software to show the change of the root mean square displacement of hydrogen atoms with the increase of time step number (see appendix). Figure 6 Based on the spatial configuration data of the model given in the kinetics.file, the OVITO atomic configuration visualization program was used to perform computer visualization image processing and analysis on the model's spatial configuration, and to derive the model at different time steps (n). steps Image of the atomic spatial configuration (hidden aluminum atoms) under (1)n steps =0 (attached) Figure 7 (A)); (2)n steps =5e+05 (attached) Figure 7 (B));(3)n steps =1e+06 (Appendix) Figure 7 (C)); (4)n steps =1.5e+06 (Appendix) Figure 7 (D));(5)n steps =2e+06 (Appendix) Figure 7 (E), Appendix Figure 8 (A)-(C)).

[0219] In steps 1 and 2 of this specific application example, apart from the newly added compute nnoad / atom and fix samd / nnoad / atom commands specifically designed for the SAMD simulation method, the use and settings of all LAMMPS commands involved are performed in the conventional manner according to the instructions in the LAMMPS user manual (see reference [2]). Some related details will not be elaborated here. Detailed instructions for these commands can be found in the LAMMPS user manual (see reference [2]), and they will not be further explained here.

[0220] It should be noted that the interatomic interaction potential function of the Al-H binary alloy system in this model is the bond angle-dependent multibody interatomic interaction potential developed by F. Apostol and Y. Mishin for the Al-H binary system (see reference

[30] ). Verification shows that this potential function can reflect the basic physical properties of Al-H alloys more realistically and has good reliability and rationality. Calculations show that this potential function can correctly give the lower energy of the tetrahedral interstitial position of hydrogen atoms in the aluminum lattice compared to the octahedral interstitial position. The given diffusion activation energy barrier of interstitial hydrogen atoms is 0.19 eV / atom, which is very close to the results of first-principles calculations and experimental measurements. It also shows that hydrogen atoms diffuse more easily in metallic aluminum.

[0221] From the appendix Figure 6 As can be seen from the above, the SAMD simulation method proposed in this invention is used to simulate the attached... Figure 5 The model shown was subjected to computer numerical simulation. As the temperature of the secondary heat bath, Ts, increased from 100K to 500K, the migration frequency of hydrogen atoms gradually increased, resulting in a considerable migration and diffusion distance (as shown in the attached figure). Figure 7 As shown in (A)-(E), it demonstrates that adjusting the value of the auxiliary heat bath temperature Ts parameter of the system can effectively control the degree of amplification and spanning on the time scale of the SAMD simulation method proposed in this invention. (See attached...) Figure 7 As can be seen from (A)-(E), as the simulation time increases, hydrogen atoms diffuse and migrate, changing from their initial relatively dispersed distribution state (see appendix). Figure 7 (A)) results in the gradual segregation of hydrogen atoms at the grain boundaries (see appendix). Figure 7 (B)-(E)). Appendix Figure 8 The results given in (A)-(C) also show that as long as the migration frequency of hydrogen atoms is high enough (e.g., by increasing the number of simulation time steps or appropriately increasing the temperature of the auxiliary heat bath in the system), the SAMD simulation method proposed in this invention can give results of hydrogen atom segregation at grain boundaries. This is consistent with the results given by the Monte Carlo (MC) method (see appendix). Figure 8 (D) is relatively consistent. The main reason for the segregation of hydrogen atoms at grain boundaries is that the activation energy barrier for the diffusion and migration of hydrogen atoms in the aluminum lattice is different from that for their back-and-forth diffusion and migration at and near grain boundaries. The SAMD simulation method proposed in this invention can give a good result of the segregation of hydrogen atoms at grain boundaries, indicating that the SAMD simulation method proposed in this invention can well realize the assignment of corresponding activation frequencies with appropriate magnitudes to structural transformation events and processes with different activation energy barriers and reaction coefficients, thereby better ensuring the physical reality of the structural transformation dynamics in matter and materials.

[0222] Simple estimations show that for SAMD simulations of the system's auxiliary heat bath temperatures Ts of 100K, 300K, and 500K, with a total time step of 2.0 × 10⁻⁶, the results are as follows: 6 Within the simulation timeframe, the average number of migration jumps for hydrogen atoms in the model were approximately 12.0, 57.6, and 132.5, respectively, while the experimentally calculated diffusion coefficient of hydrogen atoms in aluminum crystal at a thermodynamic temperature of 100 K was approximately 5.98 × 10⁻⁶. -17 m 2 / s, which is equivalent to approximately 8750 jumps per second for hydrogen atoms in the tetrahedral interstices of an aluminum crystal. Therefore, in this specific application example, the corresponding SAMD simulation speedup can be estimated to be approximately 1.4 × 10⁻⁶. 6 6.3×10 6 1.5×10 7 .

[0223] In specific implementation, the method proposed in the technical solution of this invention can be automatically executed by those skilled in the art using computer software technology. System devices for implementing the method, such as computer-readable storage media storing the corresponding computer program of the technical solution of this invention and computer equipment including the computer program running the corresponding computer program, should also be within the protection scope of this invention.

[0224] In some possible embodiments, a computer simulation apparatus for materials is provided, including a processor and a memory, wherein the memory is used to store program instructions and related input and output data, and the processor is used to invoke the instructions stored in the processor to execute an improved method for computer simulation of materials as described above.

[0225] In some possible embodiments, a computer simulation apparatus for matter and materials is provided, including a readable storage medium on which a computer program is stored, which, when executed, implements an improved method for computer simulation of matter and materials as described above.

[0226] The above description is merely a specific embodiment of the present invention and is not intended to limit the present invention. Any modifications or improvements made within the guiding principles, spirit, and rules of the present invention should be included within the scope of protection of this patent.

Claims

1. An improved method for computer simulation of substances and materials, first performing the following step 1, Step 1, for the substances and materials to be simulated by computer, constructing a multi-particle model and system of the substances and materials, achieving the following, For any substance or material containing N particles, the set of all particles is denoted as {1,2,…,N}. All particles are represented in three-dimensional Euclidean space. The coordinates in the figure are X = (x1, x2, ..., x). N ) and the corresponding particle momentum P=(p1,p2,…,p N As a dynamic quantity and constituting a phase space, it describes any microscopic state of the substance or material, where the coordinates of any particle i are x. i The momentum of particle i is p. i =m i v i m i Let v be the mass of particle i. i Let V be the velocity of particle i, i∈{1,2,…,N}, and the velocities of all particles are V=(v1,v2,…,v…). N ); Assume a three-dimensional Euclidean space The coordinate axes of the Cartesian coordinate system are x i With v i and p i All are 3-dimensional vectors, X, V, and P are all 3N-dimensional vectors; these particles interact with each other, and U(X) represents the potential energy function corresponding to the interaction potential field between all particles. Then, the force on any particle i in the α-axis direction is expressed as: The particle exists in three-dimensional Euclidean space. The coordinates along the α-axis are The particles are atoms or molecules; characterized in that Then performing the following steps 2 and 3, Step 2, using the equations of motion of the system of particles The system of equations of motion of the system of particles is used to calculate the time evolution of the coordinates X and velocities V of all particles in the multi-particle model of the substance and material. The The system of equations gives each particle in the many-particle system three additional dynamical quantities The three additional dynamical quantities of each particle form a three-dimensional vector The additional dynamical quantities of all particles form a 3N-dimensional vector S = (s1, s2, s3,..., s N ); the additional dynamical quantities of each particle are set the thermodynamic and statistical mechanical quantities corresponding to the eccentricity of the nearest neighbors of this particle ; The near-neighbor eccentricity of the particle The near-neighbor eccentricity of the particle is calculated by first calculating the three-dimensional Euclidean space vector from the geometric center of all the nearest neighbor particles to the position of the particle, and then taking the absolute value of the three coordinate axis components of the near-neighbor eccentric displacement R(X) of the particle respectively as the values of the three components of the near-neighbor eccentricity d of the particle The near-neighbor eccentricities of all the particles form a 3N-dimensional vector D=(d1, d2, d3, …, dN) N ); The near-neighbor eccentricity d of the particle is a function of the coordinates X of all the particles, and is used to describe the amplitude of the non-coordinated motion between the particle and its nearest neighbor particles; Due to the nearest neighbor eccentricity d of the particle, the particle's additional kinetic quantity s, and the particle's position in three-dimensional Euclidean space... Position x has the same physical dimensions, and the multi-particle system is considered to be in a 6-dimensional space. Motion in this 6-dimensional space; In this process, each particle is placed in the original three-dimensional Euclidean space. The nearest neighbor eccentricity of the particle as defined in [the document / reference] With the particle in extra space Additional dynamic quantities The coupling is achieved using a harmonic oscillator, and then the potential energy function corresponding to the coupling interaction potential field is superimposed onto the multi-particle system in the original three-dimensional Euclidean space. From the potential energy function U(X) corresponding to the interparticle interaction potential field, we obtain the state of all particles in 6-dimensional space. The coupled superposition potential energy function U corresponding to the interparticle interaction potential field κ (X,S); with the coupled superposition potential function U κ Based on (X,S), using... The system of equations for multi-particle system dynamics is used to calculate the time-varying coordinates (X) and velocities (V) of all particles in a multi-particle model of the substance or material. The system of equations is constructed such that the particle exists in the original three-dimensional Euclidean space. The dynamic quantity X in the particle changes with time according to the Nose-Hoover dynamic equations, while the particle in the extra space... The variation of the additional kinetic quantity S with time follows the Langevin kinetic equation; the viscosity coefficients of the Nose-Hoover kinetic equation and the Langevin kinetic equation are independent of each other, and their respective bath temperatures are also independent of each other; comprising The multi-particle system dynamics equation set is time-discretized and numerically processed, and the coordinates X and velocities V of all particles of the multi-particle model at each discrete time step are automatically and step by step calculated in time step order to obtain the changes and information of the values of the coordinates X and velocities V of all particles of the multi-particle model of the substance and material with time. Step 3, according to the values of the coordinates X and velocities V of all particles corresponding to each discrete time step of the multi-particle model of the substances and materials obtained in step 2, performing thermodynamic and dynamic analysis on the corresponding substances and materials.

2. The improved method for computer simulation of substances and materials of claim 1, wherein: In step 2, the definition and calculation expression of the near-neighbor eccentricity D = (d1, d2, d3,..., dn) of the particles are as follows, N ) Wherein, N - the total number of particles in the multi-particle system; - the indices of the three coordinate axes of the three-dimensional Euclidean space, respectively; - the indices of the three coordinate axes of the three-dimensional Euclidean space, respectively; X - coordinates of all particles in three-dimensional Euclidean space, X = (xl, x2,..., xn) , N , d i The eccentricity of the near neighbors of any particle i in the system is a 3D vector. - the near neighbor eccentricity d of i particles, respectively i In three-dimensional Euclidean space The three coordinate axes The 3 components of a direction; R i (X)—i the near neighbor eccentric displacement of particle i is a function of the coordinates X of all particles; - the respective nearest neighbor eccentric displacement R of the i-particle i (X) in three-dimensional Euclidean space of the three coordinate axes is a function of the coordinates X of all particles; respectively the respective mathematical absolute value; N i —i the number of all near neighbor particles within a radius of R d from the particle i; - the set of all near neighbors of particle i within a radius of R d as the radius of the set of all near neighbors of particle i x i — the coordinates of the i particle in the three-dimensional Euclidean space at the current moment of time; x j —At the current moment, particle j is in three-dimensional Euclidean space. The coordinates in the diagram.

3. The improved method for computer simulation of substances and materials of claim 1, wherein: In step 2, the inter-particle potential field in the 6-dimensional space corresponding to the coupling superposition potential function U κ (X,S) is used as a basis, and a multi-particle system dynamics equation set including a coupling superposition potential function U κ (X,S) and the equation set has the following expression: Wherein, N - the total number of particles in the multi-particle system; - the indices of the three coordinate axes of the three-dimensional Euclidean space, respectively; - the indices of the three coordinate axes of the three-dimensional Euclidean space, respectively; X - coordinates of all particles in three-dimensional Euclidean space, X = (xl, x2,..., x ), N ), S - all particles in extra space N ​​ - the coordinates, velocities, accelerations of the i particles in the direction of the alpha coordinate axis in three-dimensional Euclidean space, respectively, ​ - the i-particle in the additional space the additional dynamical quantity s i In components of the coordinate axes, their first derivatives with respect to time, their second derivatives with respect to time, ζ - the Nose-Hoover dynamic extension variable of the multi-particle system; - the first derivative of the Nose-Hoover kinetic extension variable ζ of the multi-particle system with respect to time; m i —i particle in three-dimensional Euclidean space of mass; μ i —i particle in extra space of virtual mass; gamma x - A multi-particle system in three-dimensional Euclidean space corresponding to Nose-Hoover dynamics; T x —Multi-particle system in three-dimensional Euclidean space corresponding to Nose-Hoover dynamics; k B Boltzmann constant; gamma s — A multi-particle system in additional space corresponding to Langevin dynamics in the presence of a friction coefficient; T s — Multi-particle system in additional space corresponding to the temperature of the secondary heat bath in Langevin dynamics; R s (t)— noise function; U(X) - the potential energy function corresponding to the field of interparticle interaction potentials of the many-particle system in three-dimensional Euclidean space is a function of the coordinates X of all particles; - the additional kinetic quantity of the particle in the component of the coordinate axis direction, the component of the near-neighbor eccentricity of a particle in the direction of the alpha coordinate axis in three-dimensional Euclidean space is a function of the coordinates X of all particles, ​ κ - the force coupling parameter between the additional dynamic quantity S of the multi-particle system and the near neighbor eccentricity D of the particles; U κ (X,S)—Multi-particle systems in 6-dimensional space The coupling potential function corresponding to the inter-particle potential field in the center is a function of the coordinates X of all particles and the additional dynamic quantity S.

4. The improved method for computer simulation of substances and materials of claim 1, wherein: In step 2, the containing The multi-particle system dynamics equation set is time-discretized and numerically processed, a discrete time step is set as Δt, the discrete time step Δt adopts a fixed time step or a non-fixed time step, and values of coordinates X and velocities V of all particles corresponding to each discrete time step of the multi-particle model of the substance and material are solved out through automatic step-by-step operation, and evolution of the microstate of the substance and material with time is given.

5. The improved method for in silico modeling of substances and materials according to any of claims 1 or 2 or 3 or 4, characterized by: For simulating and analyzing the diffusion and migration dynamics of interstitial carbon atoms and iron atom vacancies in a body-centered cubic iron lattice in iron-carbon alloys.

6. The improved method for in silico modeling of substances and materials according to any of claims 1 or 2 or 3 or 4, characterized by: For simulating and analyzing the diffusion and migration dynamics of hydrogen atoms in aluminum bicrystals and their segregation on grain boundaries in aluminum bicrystals dissolved with a certain concentration of hydrogen atoms.

7. A computer simulation modeling apparatus for substances and materials, characterized by: For implementing the improved method for computer simulation of substances and materials as claimed in any one of claims 1 to 6.

Images