Long-range many-body effect and short-range angle-dependent potential model and construction method

Abstract

The application provides a long-range many-body effect and short-range angle-dependent interaction potential model and a construction method, and belongs to the field of material simulation, and specifically comprises the following steps: firstly, fitting various potential functions, fitting target data is obtained from external input or through first-principle calculation, and the fitting target data is used as a fitting target set; fitting each potential function, preparing a predictive property calculation structure according to the fitting target set, and calculating a corresponding predictive property; the calculation result of the predictive property is consistent with each fitting target; then, using the least square method as a cost function, evaluating the potential function by calculating the energy value of each correlation model and the square sum of the deviation of each fitting target, and setting a weight factor for each fitting target; finally, setting an optimization algorithm, taking the calculation result of the potential function as a fitting target, and adopting step-by-step fitting to construct a long-range many-body effect and short-range angle-dependent interaction potential hybrid model; and the application effectively guarantees the hcp phase stability of high c / a ratio hcp metals.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of material simulation, and particularly relates to a long-range many-body effect and short-range angle-dependent interaction potential model and a construction method. BACKGROUND

[0002] The interatomic potential is a function representing the bonding or non-bonding interaction between atoms or molecules, which fundamentally determines all results of atomic-scale simulation of materials, and the development of related models and the fitting of potential functions are one of the most core tasks in molecular simulation.

[0003] At present, the long-range pair-wise functional potentials such as EAM and MEAM are generally applicable to metal systems, while the short-range cluster potentials such as Tersoff potential and SW potential are more suitable for covalent systems; however, these models also face a series of problems in their respective application scenarios, for example, the lattice stability of high c / a ratio hcp metals and the conflict between stacking fault energy order and splitting work of covalent systems, due to the fact that their respective physical basis is often limited to a single description of metal bonds or covalent bonds.

[0004] The interatomic interaction in actual materials is often complex, and sometimes exhibits both metal or covalent bond components, and there is an urgent need for an interatomic potential model that has good compatibility for both metal and covalent systems and can effectively ensure the hcp phase stability of high c / a ratio hcp metals. SUMMARY

[0005] Considering that the hybrid model is more complex than the traditional single physical model, the application provides a long-range many-body effect and short-range angle-dependent interaction potential model and a construction method, which can reproduce the correct order of

[111] face dislocation and

[111] ,

[110] and

[100] face full dislocation stacking fault energy through step-by-step fitting, while ensuring reasonable splitting work.

[0006] The long-range many-body effect and short-range angle-dependent interaction potential model has an expression as follows:

[0007]

[0008] Wherein i, j and k represent atoms in the atomic coordinate and model box structure file inputted by a user from the outside;

[0009] F i (ρ i ) is a long-range embedding energy term containing many-body effects, and adopts an EAM embedding energy expression:

[0010] F i (ρ i )=F α [1-ηlnρ i ]ρi η

[0011] where ρ i is the background charge density of central atom i, which is expressed as

[0012]

[0013] where r ij is the distance between atom i and its nearest neighbor atom j; is a smooth transition function, F α , η, f e , and β e are fitting parameters, and r e is the first nearest neighbor distance parameter of the ground state structure input by the user.

[0014] Φ ij(k) is a short-range correction term containing angular dependence, which is expressed as

[0015]

[0016] where

[0017]

[0018]

[0019] r ik is the distance between atom i and atom k in the structure input by the user, θ ijk is the included angle between coordination bonds ij and ik, A m , λ m , λ2, β, and n are fitting parameters, is a smooth transition function; g(θ ijk ) is an angular-dependent term, which is expressed as

[0020]

[0021] c, d, and h are fitting parameters.

[0022] The method for constructing the long-range many-body effect and short-range angular-dependent potential model comprises the following specific steps:

[0023] Step 1: fitting various potential functions, fitting target data is input externally or calculated by first principles, as a fitting target set;

[0024] 1.1 For fcc metal elemental potential fitting, the equilibrium lattice constant, cohesive energy, vacancy formation energy, elastic constants, unstable stacking fault energy, stable stacking fault energy,

[111] surface formation energy, fcc phase energy difference and bcc phase energy difference are inputted from outside as the fitting target set.

[0025] 1.2 For bcc metal elemental potential fitting, the equilibrium lattice constant, cohesive energy, vacancy formation energy, elastic constants,

[110] unstable stacking fault energy,

[112] unstable stacking fault energy,

[111] surface formation energy, fcc phase energy difference and hcp phase energy difference are inputted from outside as the fitting target set.

[0026] 1.3 For hcp metal elemental potential fitting, the equilibrium lattice constant, equilibrium c / a ratio, cohesive energy, vacancy formation energy, elastic constants, unstable stacking fault energy, stable stacking fault energy,

[0001] surface formation energy, fcc phase energy difference and bcc phase energy difference are inputted from outside as the fitting target set.

[0027] 1.4 For diamond structure elemental potential fitting, the equilibrium lattice constant, cohesive energy, elastic constants, non-full dislocation unstable stacking fault energy,

[111] surface formation energy, hexagonal diamond structure phase energy difference are inputted from outside as the fitting training value array.

[0028] 1.5 For cross potential fitting, the cohesive energy, equilibrium lattice constant and elastic constants of L12 (A3B), L12 (AB3), B1 (AB) and B2 (AB) structures are inputted from outside as the fitting target set.

[0029] Step two, fitting each potential function, preparing the respective predictive property calculation structure according to the fitting target set, and calculating the corresponding predicted properties; making the calculation results of the predicted properties consistent with each fitting target;

[0030] Specifically: for each metal elemental potential fitting, the predictive property calculation structure includes the fcc (face-centered cubic), bcc (body-centered cubic) and hcp (hexagonal close-packed) structures inputted from outside;

[0031] For diamond structure elemental potential fitting, the predictive property calculation structure includes the cubic diamond and hexagonal diamond structures inputted from outside by the user;

[0032] For cross potential fitting, the predictive property calculation structure includes the alloy L12 (A3B), L12 (AB3), B1 (AB) and B2 (AB) structures inputted from outside by the user.

[0033] The process of calculating the predicted properties is:

[0034] Step a, input all the predicted property calculation structures into the molecular dynamics calculation software respectively and relax to obtain the respective lattice constants, energies and stresses, and calculate the polymerization energy of each structure using the following formula:

[0035]

[0036] wherein is the energy obtained after relaxation of each structure, N i is the number of atoms contained in the structure.

[0037] Step b, based on the relaxed base structure model, calculate the elastic constant;

[0038] The formula is as follows:

[0039]

[0040] wherein ∈ = 10 -5 , and are the j components of the box stress obtained by the molecular dynamics calculation software after applying + ∈ and - ∈ strain along the i direction, wherein i and j traverse the six directions of x, y, z, yz, xz and xy.

[0041] Step c, for the metal and diamond structure element potential, fit the structure with the lowest polymerization energy, and construct the correlation model by adding various lattice and atomic transformations on the basis of its cell;

[0042] The correlation model specifically includes the following:

[0043] a) Vacancy formation energy model: 3 times expansion of the cell in x, y, z direction, and after expansion, the first atom in the original cell is deleted.

[0044] b) Stacking fault energy model:

[0045] For the fcc structure as the ground state structure, the unit cell with the orientation of [11-2], [-110] and

[111] direction is constructed based on the ground state structure, and the z direction is expanded by 6 periods; then based on the z coordinate, the part with z coordinate lower than half period of z direction is fixed, and the other part is respectively moved along x axis by 1 / 6 and 1 / 3 lattice units, to obtain the unstable stacking fault and stable stacking fault energy calculation model.

[0046] For the bcc structure as the ground state structure, including

[110] and

[112] unstable stacking fault calculation model; specifically:

[0047] Based on the ground state structure, construct the unit cell with orientation of

[111] , [11-2], [-110] direction, and expand the cell in z direction for 6 periods. Then based on the z coordinate, fix the part with z coordinate lower than the half period of z direction, and move the other part along the x axis by 1 / 6 and 1 / 3 lattice units, to obtain the

[112] unstable stacking fault calculation model.

[0048] Based on the ground state structure, construct the unit cell with orientation of

[111] , [11-2], [-110] direction, and expand the cell in z direction for 6 periods. Then based on the z coordinate, fix the part with z coordinate lower than the half period of z direction, and move the other part along the x axis by 1 / 6 and 1 / 3 lattice units, to obtain the

[112] unstable stacking fault calculation model.

[0049] Based on the ground state structure, construct the unit cell with orientation of

[111] , [11-2], [-110] direction, and expand the cell in z direction for 6 periods. Then based on the z coordinate, fix the part with z coordinate lower than the half period of z direction, and move the other part along the x axis by 1 / 6 and 1 / 3 lattice units, to obtain the

[112] unstable stacking fault calculation model.

[0050] Based on the ground state structure, construct the unit cell with orientation of

[111] , [11-2], [-110] direction, and expand the cell in z direction for 6 periods. Then based on the z coordinate, fix the part with z coordinate lower than the half period of z direction, and move the other part along the x axis by 1 / 6 and 1 / 3 lattice units, to obtain the

[112] unstable stacking fault calculation model.

[0051] c) Surface formation model:

[0052] Based on the relaxed ground state structure, construct a model with z direction orientation as close-packed direction, and expand the cell in z direction for 10 periods; then add a vacuum layer with thickness of 10 A on the top of the model as the surface energy calculation model.

[0053] Step d, for each element potential fitting, use molecular dynamics to calculate the energy of each associated model;

[0054] Specifically, the following calculations are included:

[0055] a) After obtaining the relaxation energy by fixing the lattice and relaxing the atoms, the vacancy formation energy is calculated using the following formula:

[0056]

[0057] Where E tot is the energy of the vacancy formation energy model, N tot is the number of atoms included in the vacancy formation energy model, ​is the energy of the ground structure, and N0 is the number of atoms in the ground structure.

[0058] b) Fixing the lattice, relaxing the atoms along the z direction, and obtaining the relaxation energy, the stable or unstable stacking fault energy is calculated using the following formula:

[0059]

[0060] where E sf is the energy of the stable or unstable stacking fault model after relaxing along the z direction, N sf is the number of atoms contained in the corresponding model.

[0061] c) Fixing the lattice, only relaxing the atoms, and obtaining the relaxation energy, the surface formation energy is calculated using the following formula:

[0062]

[0063] where E surf is the energy of the surface atom of the surface formation model after relaxing, N surf is the number of atoms contained in the corresponding model.

[0064] d) The energy difference between the remaining base structure and the ground state structure except the ground state structure is calculated using the following formula:

[0065]

[0066] where E is the energy of the i-th structure after relaxing, is the number of atoms contained in the corresponding model.

[0067] Step three, using the least square method as the cost function, by calculating the predicted energy value of each associated model and the square sum of the deviation of each fitting target to evaluate the potential function, and setting the weight factor for each fitting target;

[0068] The potential function takes the following form:

[0069]

[0070] where N target is the number of fitting training values, is the i-th fitting target value in the fitting target set, is the set predictive property value corresponding to the i-th fitting target, w i is the weight factor corresponding to the i-th fitting target, T i is the proportionality coefficient.

[0071] The higher weight is adopted for the polymerization energy and the lattice constant, and the smaller weight is adopted for other fitting targets, such as the vacancy formation energy, the elastic constant, the stacking fault energy and the surface energy, and the weight values of different fitting targets can be fine-tuned according to the fitting effect.

[0072] Step four, setting an optimization algorithm, taking the calculation result of the potential function as a fitting target, and adopting a step-by-step fitting to construct a hybrid model;

[0073] Specifically includes the following steps:

[0074] 4.1 First, globally fitting the undetermined fitting parameters of the angle-dependent correction term in the hybrid model:

[0075] The fitting variable is the short-range correction term Φ ij(k) defined in the undetermined fitting parameters:

[0076]

[0077] The initial value of the fitting variable is initially defined, the global particle swarm optimization algorithm is adopted, the fitting termination tolerance threshold is artificially set, and the output quantity is the fitting result of the current fitting variable.

[0078] The undetermined fitting parameters refer to A0, A1, λ0, λ1, λ2, β, n, c, d and h in the hybrid model.

[0079] 4.2 Then, locally fitting the undetermined fitting parameters of the angle-dependent correction term in the hybrid model:

[0080] The fitting variable remains unchanged, the fitting result output by 4.1 is taken as the initial fitting variable of this round, the fitting target remains unchanged, the local Nelder-Mead simplex algorithm is adopted, 5% of each component is first added to the initial fitting parameter to generate a simplex around the initial fitting parameter, and then the simplex is repeatedly modified according to the Nelder-Mead simplex optimization algorithm until the set fitting termination tolerance is reached, and the output quantity is the fitting result of the current fitting variable.

[0081] 4.3 Then, globally fitting all the undetermined parameters in the hybrid model:

[0082] The fitting result output by 4.2 is taken as the fitting variable initial value, the global particle swarm optimization algorithm is adopted, and the output quantity is the fitting result of the current fitting variable; all the undetermined parameters include F α , η, f e , β e , A0, A1, λ0, λ1, λ2, β, n, c, d and h.

[0083] 4.4 Then, locally fitting all the undetermined parameters in the hybrid model:

[0084] The fitting variable is not changed, and the fitting result output by 4.3 is taken as the initial value of the fitting variable in this round of fitting, and the target of fitting is not changed. The local Nelder-Mead simplex algorithm is adopted, 5% of each component is added to the initial fitting parameter to generate a simplex around the initial fitting parameter, and then the simplex is repeatedly modified according to the Nelder-Mead simplex optimization algorithm until the set fitting termination tolerance is reached, and the output quantity is the final fitting result of the current fitting variable.

[0085] 4.5 The final fitting result of the fitting variable is output according to the list of molecular dynamics software and saved as a potential function.

[0086] The beneficial effects of the present application are as follows:

[0087] The present application discloses a long-range many-body effect and short-range angle-dependent interaction potential model and a construction method, which breaks through the limitation that traditional potential function models can only describe metal or covalent single systems, has good compatibility for metal and covalent systems, effectively guarantees the hcp phase stability of high c / a ratio hcp metal, reliably reproduces the correct order of

[111] surface dislocation and

[111] ,

[110] and

[100] surface full dislocation stacking fault energy of diamond structure, and ensures reasonable splitting work. BRIEF DESCRIPTION OF DRAWINGS

[0088] Figure 1 The flowchart of the long-range many-body effect and short-range angle-dependent interaction potential model and the construction method of the present application.

[0089] Figure 2 The generalized stacking fault energy curve of the diamond and silicon hybrid potential constructed by the present application, and the comparison with the DFT calculation value and the previously published potential function.

[0090] Figure 3 The lattice constant, c / a ratio, cohesive energy, phase energy difference, vacancy formation energy, unstable stacking fault energy, stable stacking fault energy, elastic constant and surface formation energy of the Zn hybrid potential constructed by the present application, and the comparison with the DFT calculation value, experimental data and previously published potential function. DETAILED DESCRIPTION

[0091] In order to facilitate those skilled in the art to understand and implement the present application, the present application will be further described in detail below in combination with the drawings and examples.

[0092] The application provides a long-range many-body effect and short-range angle-dependent potential model and a construction method, the model is a hybrid interatomic potential model composed of a long-range many-body embedding energy item based on a density functional theory and a short-range angle-dependent correction item based on a valence electron pair repulsion theory, and has good compatibility for metal and covalent systems; considering that the hybrid model is more complex than a traditional single physical model, the application also provides a step-by-step fitting technology combining a global optimization algorithm and a local optimization algorithm, fitting the short-range angle-dependent item first and then constructing the hybrid potential based on the short-range angle-dependent item. The hybrid interatomic potential is not only suitable for metal and covalent material systems, but also solves a series of bottleneck problems of traditional long-range many-body potential and short-range angle-dependent potential functions, and lays a foundation for further improving the accuracy of interatomic potential and improving the reliability of related simulation.

[0093] The long-range many-body effect and short-range angle-dependent potential model has an expression as follows:

[0094]

[0095] Wherein i, j and k represent atoms in a box structure file inputted by a user from outside and containing atomic coordinates and model;

[0096] F i (ρ i ) is a long-range embedding energy item containing a many-body effect, and an embedding energy expression in an EAM (embedded-atom method) embedding energy model and a MEAM (modified embedded-atom method) embedding energy model is adopted:

[0097]

[0098] Wherein α and β are initially taken as 1.0 and 0.0 respectively, and ρ i is a background charge density of the central atom i, and an expression thereof is

[0099]

[0100] Wherein r ij is a distance between the atom i and a near-neighbor atom j; is a smooth transition function, in the application, when r ij >r c2 , f(r ij ) is 0, when r c1 <r ij , f(r c1 ) is 1, and when r c2 is between r α and r e , the transition is as follows:

[0101]

[0102] F α , η, f e β e and r e The first nearest neighbor distance parameter is the ground state structure input from the user; the initial values ​​are 2.5, 3.0, 1.0, 1.5, and 5.0, r c2 =0.98r c1 ;

[0103] After the fitting is complete, the user can set β to 1.0 as needed to improve the fitting accuracy. In this case, the angular dependence contribution of the electron is taken into account. The expression is

[0104]

[0105] and The fitting parameters are initially set to 1.0 and the polymerization energy of the ground-state structure as an external input.

[0106] The expressions for the electron densities in the summation term are:

[0107]

[0108]

[0109]

[0110]

[0111] Among them, α, β and γ are ergodic coordinate bonds r ij The x, y, and z components of Cartesian coordinates, and ρ i (k) The expression is:

[0112]

[0113] Where β0, β1, β2, β3, t0, t1, t2 and t3 are fitting parameters, with initial values ​​of 4.0, 4.0, 5.0, 3.0, 1.0, 2.0, 2.0 and -1.0.

[0114] Φ ij(k) It is a short-range correction term that includes angle dependency effects, expressed as follows:

[0115]

[0116] in

[0117]

[0118]

[0119] r ik is the distance between atom i and atom k in the user input structure, θ ijk is the angle between the bonds ij and ik, A m , λ m , λ2 and β are fitting parameters with initial values of 150, 250, 1.4, 2.0, 0.7 and 10 -7 . is a smooth transition function; in the present invention it is 0 when r ij > R + D , r ij < R - D , r ij is between R - D and R + D, it is transitioned as follows:

[0120]

[0121] where R and D depend on the user-supplied input ground state structure, let the first and second neighbor distances of the ground state structure be r e and r2, then g(θ ijk ) is an angular dependence term, which is expressed as:

[0122]

[0123] c, d and h are fitting parameters; initial values are 10 4 , 5.0, and -0.5.

[0124] After the fitting is done, the user can replace the Stillinger-Weber model with:

[0125]

[0126] where

[0127]

[0128]

[0129] where A, B, p, q, λ, ω and σ are fitting parameters with initial values of 15.0, 1.0, 4.0, 0.0, 20.0, -0.3 and 2.0.

[0130] Since the hybrid model contains both long-range many-body effect and short-range angle-dependent effect, it is more complex than the traditional single model, which leads to the fact that the traditional potential function fitting scheme is not applicable to the construction of the hybrid model. Figure 1 As shown in the following steps:

[0131] Step one, fitting various potential functions, fitting target data is inputted from outside or calculated by first principle calculation as fitting target set;

[0132] 1.1 For fcc metal element potential fitting, input the equilibrium lattice constant, cohesive energy, vacancy formation energy, elastic constant, unstable stacking fault energy, stable stacking fault energy,

[111] surface formation energy, hcp phase energy difference and bcc phase energy difference from outside as fitting target set.

[0133] 1.2 For bcc metal element potential fitting, input the equilibrium lattice constant, cohesive energy, vacancy formation energy, elastic constant,

[110] surface unstable stacking fault energy,

[112] surface unstable stacking fault energy,

[111] surface formation energy, fcc phase energy difference and hcp phase energy difference from outside as fitting target set.

[0134] 1.3 For hcp metal element potential fitting, input the equilibrium lattice constant, equilibrium c / a ratio, cohesive energy, vacancy formation energy, elastic constant, unstable stacking fault energy, stable stacking fault energy,

[0001] surface formation energy, fcc phase energy difference and bcc phase energy difference from outside as fitting target set.

[0135] 1.4 For diamond structure element potential fitting, input the equilibrium lattice constant, cohesive energy, elastic constant, non-full dislocation unstable stacking fault energy,

[111] surface formation energy, hexagonal diamond structure phase energy difference from outside as fitting training value array.

[0136] 1.5 For cross potential fitting, input the aggregation energy, equilibrium lattice constant and elastic constant of L12(A3B), L12(AB3), B1(AB) and B2(AB) structures from outside as fitting target set.

[0137] Step two, fitting each potential function, according to the fitting target set, preparing the corresponding predicted property calculation structure and calculating the corresponding predicted property, so that the calculation result of the predicted property is consistent with each fitting target;

[0138] Specifically, for each metal element potential fitting, the predicted property calculation structure includes fcc (face-centered cubic), bcc (body-centered cubic) and hcp (hexagonal close-packed) structures inputted from outside;

[0139] For diamond structure element potential fitting, the predictive property calculation structure includes that the user inputs cubic diamond and hexagonal diamond structure from outside;

[0140] For cross potential fitting, the predictive property calculation structure includes that the user inputs alloy L12 (A3B), L12 (AB3), B1 (AB) and B2 (AB) structure from outside.

[0141] The calculation of the predicted property process is:

[0142] Step a, input all the predictive property calculation structures into the molecular dynamics calculation software (such as LAMMPS) respectively and relax to obtain the respective lattice constants, energies and stresses, and calculate the aggregation energy of each structure using the following formula:

[0143]

[0144] Wherein is the energy obtained after relaxation of each structure, N i is the number of atoms contained in the structure.

[0145] Step b, based on the relaxed basic structure model, calculate the elastic constant;

[0146] The formula is as follows:

[0147]

[0148] Wherein ∈ = 10 -5 , and are the j components of the box stress obtained by the molecular dynamics calculation software after applying + ∈ and - ∈ strain along the i direction, respectively, wherein i and j traverse the six directions of x, y, z, yz, xz and xy.

[0149] Step c, for metal and diamond structure element potential fitting, the structure with the lowest aggregation energy (the most stable in the basic structure) is taken as the ground state structure, and various lattice and atomic transformations are added on the basis of its cell to construct a correlation model; to cover various possible potential intermediate structures related to phase transition or crystal plastic deformation;

[0150] The correlation model specifically includes the following:

[0151] a) Vacancy formation energy model: 3 times of cell expansion in x, y, z direction, and after expansion, the first atom in the initial cell is deleted.

[0152] b) Stacking fault energy model:

[0153] For the fcc structure, the unit cells with orientations of [11-2], [-110] and

[111] are constructed based on the ground structure, and the z direction is expanded by 6 periods. Then, based on the z coordinate, the part with a z coordinate lower than half of the z direction is fixed, and the other part is moved by 1 / 6 and 1 / 3 lattice units along the x axis, respectively, to obtain the unstable stacking fault and stable stacking fault energy calculation model.

[0154] For the bcc structure, the

[110] and

[112] unstable stacking fault calculation models are included. Specifically:

[0155] The unit cells with orientations of

[111] , [11-2] and [-110] are constructed based on the ground structure, and the z direction is expanded by 6 periods. Then, based on the z coordinate, the part with a z coordinate lower than half of the z direction is fixed, and the other part is moved by 1 / 2 lattice unit along the x axis, to obtain the

[110] unstable stacking fault calculation model.

[0156] The unit cells with orientations of

[111] , [1-10] and [11-2] are constructed based on the ground structure, and the z direction is expanded by 6 periods. Then, based on the z coordinate, the part with a z coordinate lower than half of the z direction is fixed, and the other part is moved by 1 / 6 and 1 / 3 lattice units along the x axis, to obtain the

[112] unstable stacking fault calculation model.

[0157] For the hcp structure, the unit cells with orientations of [-1-120], [1-100] and

[0001] are constructed based on the ground structure, and the z direction is expanded by 6 periods. Then, based on the z coordinate, the part with a z coordinate lower than half of the z direction is fixed, and the other part is moved by 1 / 6 and 1 / 3 lattice units along the x axis, respectively, to obtain the unstable stacking fault and stable stacking fault energy calculation model.

[0158] For the diamond structure, the unit cells with orientations of [11-2], [-110] and

[111] are constructed based on the ground structure, and the z direction is expanded by 6 periods. Then, based on the z coordinate, the part with a z coordinate lower than half of the z direction is fixed, and the other part is moved by 1 / 6 lattice unit along the x axis, to obtain the non-full dislocation unstable stacking fault energy calculation model.

[0159] c) Surface formation model:

[0160] Based on the relaxed ground structure, a model with z direction orientation as close-packed direction is constructed, and the z direction is expanded by 10 periods; then a vacuum layer with a thickness of 10 is added on the top of the model as the surface energy calculation model.

[0161] Step d, for each element potential fitting, all the constructed correlation models are calculated by molecular dynamics software (such as LAMMPS) to obtain their respective energies;

[0162] Specifically, the following calculation is included:

[0163] a) Fixing the lattice, only relaxing the atoms, after obtaining the relaxation energy, the following formula is used to calculate the vacancy formation energy:

[0164]

[0165] Where E tot is the energy of the vacancy formation model, N tot is the number of atoms included in the vacancy formation model, is the energy of the ground state structure, and N0 is the number of atoms in the ground state structure.

[0166] b) Fixing the lattice, relaxing the atoms along the z direction, after obtaining the relaxation energy, the following formula is used to calculate the stable or unstable stacking fault energy:

[0167]

[0168] Where E sf is the energy of the stable or unstable stacking fault model after relaxing along the z direction, N sf is the number of atoms included in the corresponding model.

[0169] c) Fixing the lattice, only relaxing the atoms, after obtaining the relaxation energy, the following formula is used to calculate the surface formation energy:

[0170]

[0171] Where E surf is the energy of the surface formation model after relaxing the surface atoms, N surf is the number of atoms included in the corresponding model.

[0172] d) The energy difference between the remaining base structure and the ground state structure is calculated using the following formula:

[0173]

[0174] Where is the energy of the i-th structure after relaxation, excluding the ground state structure, is the number of atoms included in the corresponding model.

[0175] Step three, using the least squares method as the cost function, by calculating the predicted energy value of each correlation model and the square sum of the deviation of each fitting target to evaluate the potential function, and setting the weight factor for each fitting target;

[0176] The fitting process is carried out around the minimization of the difference between the predictive properties of the potential function and the fitting target;

[0177] The potential function takes the following form:

[0178]

[0179] where N target is the number of training values fitted in step one, is the ith fitted target value in the fitting target set, is the predicted property value corresponding to the ith fitting target in step two, w i is the weight factor corresponding to the ith fitting target, T i is the proportionality coefficient, generally equal to If T i takes the value of 0, and T ij(k) takes the value of 1.

[0180] A higher weight of 100 is given to the aggregation energy and lattice constant, and a smaller weight of 10 is given to other fitting targets such as vacancy formation energy, elastic constant, stacking fault energy and surface energy. The weight values of different fitting targets can be fine-tuned according to the fitting effect, for example, the weight is doubled when the fitting deviation of a certain property is too large (higher than 20%) to achieve better fitting effect.

[0181] Step four, set the optimization algorithm, take the calculation result of the potential function as the fitting target, and use step-by-step fitting to construct the hybrid model;

[0182] Specifically includes the following steps:

[0183] 4.1 First, globally fit the undetermined fitting parameters of the angle-dependent correction term in the hybrid model:

[0184] The fitting variable is the undetermined fitting parameter defined in the short-range correction term Φ ij(k) containing angle-dependent effect:

[0185]

[0186] The initial value of the fitting variable is initially defined, the fitting target is the calculated value of the potential function, and the global particle swarm optimization algorithm is used, wherein the particle number is set to 50, the fitting termination tolerance threshold is artificially set to 10 -4 , and the output is the fitting result of the current fitting variable.

[0187] The undetermined fitting parameters refer to A0, A1, λ0, λ1, λ2, β, c, d and h in the hybrid model.

[0188] 4.2 Then locally fit the undetermined fitting parameters of the angle-dependent correction term in the hybrid model:

[0189] The fitting variable is unchanged, the fitting result output by 4.1 is used as the initial value of the fitting variable, the fitting target is unchanged, the local Nelder-Mead simplex algorithm is used, 5% of each component is added to the initial fitting parameter to generate a simplex around the initial fitting parameter, then the simplex is repeatedly modified according to the Nelder-Mead simplex optimization algorithm until the set fitting termination tolerance is reached, and the output quantity is the fitting result of the current fitting variable.

[0190] 4.3 Then globally fit all undetermined parameters in the hybrid model:

[0191] The fitting variable is unchanged, the initial value is defined as the initial value of the fitting variable, the fitting result output by 4.2 is used as the initial value of the fitting variable, the fitting target is unchanged, the global particle swarm optimization algorithm is used, the number of particles is set to 50, and the fitting termination tolerance is set to 10 -4 The fitting is performed, and the output quantity is the fitting result of the current fitting variable;

[0192] All undetermined parameters include: F α , η, f e , β e , A0, A1, λ0, λ1, λ2, β, c, d and h.

[0193] 4.4 Then locally fit all undetermined parameters in the hybrid model:

[0194] The fitting variable is unchanged, and the fitting result output by 4.3 is used as the initial value of the fitting variable, the fitting target is unchanged. The local Nelder-Mead simplex algorithm is used, 5% of each component is added to the initial fitting parameter to generate a simplex around the initial fitting parameter, then the simplex is repeatedly modified according to the Nelder-Mead simplex optimization algorithm until the set fitting termination tolerance 10 -4 is reached, and the output quantity is the final fitting result of the current fitting variable.

[0195] 4.5 The final fitting result of the fitting variable is output according to the list of molecular dynamics software, and is saved as a potential function.

[0196] For example Figure 2As shown, the generalized stacking fault energy curve of the diamond and silicon hybrid potential constructed by the present application is consistent with the first-principles calculation. In contrast, the diamond potential constructed by Erhart and Albe et al. has a non-physical concave, which may cause a non-physical phase transition in the simulation of mechanical properties, and is inconsistent with the actual results. Although the shape of the generalized stacking fault energy curve of the Tersoff potential published in 1994 is reasonable, the unstable stacking fault energy of the incompleted dislocation on the

[112] plane is higher than that of the

[111] full dislocation, which is inconsistent with the first-principles calculation. In comparison, the diamond hybrid potential constructed by the present application does not have these problems and is suitable for the simulation of mechanical properties. In addition, for the silicon potential function, the two previously published potential functions have a non-physical concave, while the hybrid potential constructed by the present application can give a reasonable stacking fault energy curve, ensuring the reliability of the simulation results.

[0197] In addition, as Figure 3 As shown, for the metal Zn, due to its excessively high c / a ratio and stacking fault energy, it has been a difficult point for potential function fitting for a long time. The previously published potential functions either only meet the high c / a ratio but cannot guarantee that the hcp phase is the most stable and give reasonable stacking fault energy values, such as the EAM potential constructed by Sheng et al., the FS potential constructed by Igarashi et al., and the MEAM potential constructed by Dickel et al.; or although the hcp phase can be guaranteed to be the most stable, it cannot meet the high c / a ratio, for example, the MEAM potential constructed by Jang et al., the Tersoff potential constructed by Erhart et al., and the BOP potential constructed by Ward et al., which may cause a non-physical twin. In contrast, the hybrid potential constructed by the present application can reflect a reasonable c / a ratio while guaranteeing the stability of the hcp phase and giving a reasonable stacking fault energy value compared with the DFT calculation results.

Claims

1. A long-range many-body effect and short-range angle-dependent potential model, characterized in that, The potential model expression is: wherein and respectively represent atoms contained in the atomic coordinates and the model box structure file inputted from the outside by the user; is the long-range embedding energy term including many-body effects, expressed using the EAM embedding energy expression: wherein the central atom background charge density, expressed as wherein is an atom and its nearest neighbor atoms distance between; is a smooth transition function; , , and are fitting parameters, r e is a first-neighbor distance parameter of the ground state structure input externally by a user is a short-range correction term that includes angular dependence, expressed as follows: Wherein the distance between atoms in the user input structure and atoms , the angle between coordination bonds and , , , , and are fitting parameters, and are smoothing transition functions; is an angular dependency term, which is expressed as follows: wherein, , and are fitting parameters; the model is simultaneously applicable to metallic and covalent material systems, whose construction proceeds as follows: First, for various metal elements, diamond structure elements, and cross potential function fitting, the target data is fitted from external input or by first principle calculation, as the fitting target set; Then, for each potential function fitting, the corresponding predicted property calculation structure is prepared according to the fitting target set, and the corresponding predicted property is calculated; so that the calculation results of the predicted property are consistent with each fitting target; The calculation process of the predicted property is as follows: Step a, all the predicted property calculation structures were input into the molecular dynamics calculation software respectively and relaxed to obtain the respective lattice constants, energies and stresses, and the polymerization energy of each structure was calculated using the following formula: ; the energy obtained after relaxing each structure, the number of atoms contained in the structure; Step b, based on the relaxed basic structure model, calculate the elastic constant; The formula is as follows: wherein, and are applied along the directions and after straining, the component of the box stress, =10 -5 wherein and traverse the six directions x, y, z, yz, xz, and xy; Step c, for metal and diamond structure element potential fitting, the lowest energy structure is aggregated, and various lattice and atomic transformations are added on the basis of the unit cell to construct a correlation model; The correlation model specifically includes the following: a) Vacancy formation energy model: 3 times expansion of the unit cell in x, y, z direction, and deletion of the first atom in the initial unit cell after expansion; b) Layer fault energy model: For the fcc structure as the ground state structure, the unit cell with the orientation of [11-2], [-110], [111] direction is constructed based on the ground state structure, and the z direction is expanded by 6 periods; then based on the z coordinate, the part below the z direction half period is fixed, and the other part is moved along the x axis by 1 / 6 and 1 / 3 lattice units respectively, to obtain the unstable layer fault and stable layer fault energy calculation model; For the bcc structure as the ground state structure, including [110] and [112] unstable layer fault calculation model; Specifically: Based on the ground state structure, the unit cell with the orientation of [111], [11-2], [-110] direction is constructed, and the z direction is expanded by 6 periods; then based on the z coordinate, the part below the z direction half period is fixed, and the other part is moved along the x axis by 1 / 2 lattice unit, to obtain the [110] unstable layer fault calculation model; Based on the ground state structure, the unit cell with the orientation of [111], [1-10], [11-2] direction is constructed, and the z direction is expanded by 6 periods; then based on the z coordinate, the part below the z direction half period is fixed, and the other part is moved along the x axis by 1 / 6 and 1 / 3 lattice units respectively, to obtain the [112] unstable layer fault calculation model; For the hcp structure as the ground state structure, based on the ground state structure, the unit cell with the orientation of [-1-120], [1-100], [0001] direction is constructed, and the z direction is expanded by 6 periods; then based on the z coordinate, the part below the z direction half period is fixed, and the other part is moved along the x axis by 1 / 6 and 1 / 3 lattice units respectively, to obtain the unstable layer fault and stable layer fault energy calculation model; For the diamond structure as the ground state structure, based on the ground state structure, the unit cell with the orientation of [11-2], [-110], [111] direction is constructed, and the z direction is expanded by 6 periods; then based on the z coordinate, the part below the z direction half period is fixed, and the other part is moved along the x axis by 1 / 6 lattice unit, to obtain the not full dislocation unstable layer fault energy calculation model; c) Surface formation model: Based on the relaxed ground state structure, a model with z direction orientation as close-packed direction is constructed, and the z direction is expanded by 10 periods; then a vacuum layer with thickness of 10 A is added on the top of the model as the surface energy calculation model; ​ Step d, for each element potential fitting, all the constructed correlation model is calculated by molecular dynamics energy of each; Specifically includes the following calculations: a) fixed lattice, only relax the atoms, get relaxation energy, and use the following formula to calculate the vacancy formation energy: wherein is the energy of the vacancy formation energy model, is the number of atoms included in the vacancy formation energy model, is the energy of the ground state structure, is the number of atoms of the ground state structure; b) fixed lattice, atoms along the z direction relaxation, get relaxation energy, and use the following formula to calculate the stable or unstable stacking fault energy: wherein is the energy of the stable or unstable stacking fault energy model after relaxation along the z direction, is the number of atoms contained in the respective model; c) fixed lattice, only relax the atoms, get relaxation energy, and use the following formula to calculate the surface formation energy: wherein, to calculate the energy of the surface atoms of the model surface after relaxation, the number of atoms comprised by the respective model; d) use the following formula to calculate the energy difference between the remaining basic structure and the ground state structure: wherein is the energy of the structure after relaxation of the ground structure, is the number of atoms comprised by the respective model; Then, using the least square method as the cost function, by calculating the predicted energy value of each correlation model and the square sum of the deviation of each fitting target, the potential function is evaluated, and each fitting target is set with its own weight factor; Finally, set the optimization algorithm, take the calculation result of the potential function as the fitting target, and use step-by-step fitting to construct the hybrid model; Specifically includes the following steps: 4.1 first global fitting hybrid model in the angle dependent correction term pending fitting parameters: Fitted variables are short-range correction terms including angular dependence Pending fitted parameters defined in the middle: The initial value of the fitting variable is defined initially, the global particle swarm optimization algorithm is used, the fitting termination tolerance threshold is artificially set, and the output quantity is the fitting result of the current fitting variable; The pending fitting parameters are the ones in the hybrid model: , , , , , , , c, d and h; 4.2 then local fitting hybrid model in the angle dependent correction term pending fitting parameters: The fitting variable remains unchanged, the fitting result output by 4.1 is taken as the initial value of the fitting variable, the fitting target remains unchanged, the local Nelder-Mead simplex algorithm is used, first add 5% of each component to the initial fitting parameter to generate a simplex around the initial fitting parameter, then repeatedly modify the simplex according to the Nelder-Mead simplex optimization algorithm until the set fitting termination tolerance is reached, and the output quantity is the fitting result of the current fitting variable; 4.3 then global fitting hybrid model in all pending parameters: The fitting result of the output of 4.2 is taken as the initial value of the fitting variable, a global particle swarm optimization algorithm is used, and the output quantity is the fitting result of the current fitting variable; all the undetermined parameters include: 、 、 、 、 、 、 、 、 、 、 、c、d and h; 4.4 then local fitting hybrid model in all pending parameters: The fitting variable remains unchanged, and the fitting result output by 4.3 is taken as the initial value of the fitting variable, and the fitting target remains unchanged; the local Nelder-Mead simplex algorithm is used, first add 5% of each component to the initial fitting parameter to generate a simplex around the initial fitting parameter, then repeatedly modify the simplex according to the Nelder-Mead simplex optimization algorithm until the set fitting termination tolerance is reached, and the output quantity is the final fitting result of the current fitting variable; 4.5 the final fitting result of the fitting variable is output according to the list of molecular dynamics software and saved as the potential function.

2. A long-range many-body effect and short-range angle-dependent interaction potential model as claimed in claim 1, wherein, The specific formation process of the fitting target is: 1.1 for fcc metal element potential fitting, the equilibrium lattice constant, cohesive energy, vacancy formation energy, elastic constant, unstable stacking fault energy, stable stacking fault energy, [111] surface formation energy, HCP phase energy difference and BCC phase energy difference are input from the outside as the fitting target set; 1.2 For bcc metal elemental potential fitting, the equilibrium lattice constant, cohesive energy, vacancy formation energy, elastic constants, [110] unstable stacking fault energy, [112] unstable stacking fault energy, [111] surface formation energy, fcc phase energy difference and hcp phase energy difference are inputted from outside as the fitting target set; 1.3 For hcp metal elemental potential fitting, the equilibrium lattice constant, equilibrium c / a ratio, cohesive energy, vacancy formation energy, elastic constants, unstable stacking fault energy, stable stacking fault energy, [0001] surface formation energy, fcc phase energy difference and bcc phase energy difference are inputted from outside as the fitting target set; 1.4 For diamond structure elemental potential fitting, the equilibrium lattice constant, cohesive energy, elastic constants, non-full dislocation unstable stacking fault energy, [111] surface formation energy, hexagonal diamond structure phase energy difference are inputted from outside as the fitting training value array; 1.5 For cross potential fitting, the cohesive energy, equilibrium lattice constant and elastic constant of L12 (A3B), L12 (AB3), B1 (AB) and B2 (AB) structure are inputted from outside as the fitting target set.

3. A long-range many-body effect and short-range angle-dependent interaction potential model as in claim 1, wherein, The predicted property calculation structure of each potential function fitting is specifically: For each metal element potential fitting, the predicted property calculation structure includes the fcc (face-centered cubic), bcc (body-centered cubic) and hcp (hexagonal close-packed) structures inputted from outside; For diamond structure elemental potential fitting, the predicted property calculation structure includes the cubic diamond and hexagonal diamond structures inputted from outside by the user; For cross potential fitting, the predicted property calculation structure includes the alloy L12 (A3B), L12 (AB3), B1 (AB) and B2 (AB) structures inputted from outside by the user; wherein, is the number of fitted training values, is the ith fitted target value in the set of fitted targets, is the set predictive property value corresponding to the ith fitted target, w i is the weight factor corresponding to the ith fitted target, is the proportionality coefficient; The aggregation energy and lattice constant are given a higher weight, and other fitting targets, such as vacancy formation energy, elastic constant, stacking fault energy and surface energy, are given a smaller weight. The weight values of different fitting targets are fine-tuned according to the fitting effect.

Images