A Material Force Field Fitting Method Based on Potential Energy Surface Matching

Abstract

This invention discloses a material force field fitting method based on potential energy surface matching, comprising: S1 obtaining the corresponding chemical structural formula of the target material; S2 constructing a chemical structural model of the connecting bonds in the chemical structural formula based on the atomic hybridization state of the chemical structural formula; S3 performing preliminary structural optimization on the chemical structural model to obtain the energy-minimizing structure; S4 calculating the chemical structural parameters of the target material using the DFT method based on the energy-minimizing structure; S5 obtaining the Morse bond stretching potential parameter using the Morse algorithm based on the chemical structural parameters and saving it in the molecular simulation file; S6 repeating S2-S5 until the calculation of all chemical structural models is completed, and integrating all molecular simulation files into the force field file; S7 performing molecular simulation based on the force field file obtained in S6 to calculate the mechanical properties of the target material. This method, by introducing the Morse bond stretching potential parameter, makes the simulation results of fiber materials more consistent with reality.

Description

Technical Field

[0001] This invention relates to the field of materials informatics, and in particular to a material force field fitting method based on potential energy surface matching. Background Technology

[0002] Based on quantum chemistry and molecular mechanics, micro- and nanoscale integrated computational methods are gradually replacing the high-cost trial-and-error methods in traditional laboratories. When applied to atomic-level molecular dynamics models under the classical Newtonian mechanics framework such as LAMMPS, these methods can accurately describe microstructure changes, visualize interaction mechanisms, and predict the thermodynamic, transport, and constitutive properties of materials. This helps guide the macroscopic performance regulation and high-throughput screening of materials at the atomic level.

[0003] Fibers are solid structural materials with microscopic transverse scales and macroscopic axial scales. They are characterized by flexibility, fineness, and high aspect ratio. They are widely used in textile structures and fiber-reinforced composite materials. In engineering, they have designability and huge application potential. For common reinforcing fibers such as carbon fiber, glass fiber, and aramid fiber, their yield strength, breaking strength, and breaking elongation are important indicators for evaluating the reinforcement effect.

[0004] Molecular dynamics treats covalent bonds as simple harmonic oscillators with quantum effects and uses force fields to empirically fit the potential energy surface of the system's three-dimensional configuration. Therefore, the quality of the force field determines the accuracy and guiding effect of the simulation results. Common classical force fields include CFF91, AMBER, CHARMm, and CVFF; thermal programming force fields include UFF and VALBOND; and second-generation force fields include PCFF, MMFF, and COMPASS. These force fields have all achieved good results in research, but many parameters are not publicly available or lack effective calibration, making it difficult to accurately assess the reliability of the parameters.

[0005] The tensile changes in materials can be decomposed into bond interactions (bond potentials) and non-bond interactions (pair potentials). Bond interactions include bond stretching potentials, bond angle bending energies, dihedral angle torsional energies, out-of-plane angle torsional energies, and more complex cross-terms of coupling effects. Non-bond interactions include van der Waals interactions, Coulomb interactions, hydrogen bonds, etc. Bond stretching potentials often play a dominant role in the system. Traditional harmonic bond potentials cannot describe the bond breaking process, and simply using the sum of atomic covalent bond radii as a bond breaking criterion is too inaccurate.

[0006] Patent document CN109994158A discloses a system and method for constructing a molecular reaction force field based on reinforcement learning, including an input / output module, a parameter and configuration module, a molecular dynamics interface module, an environment setting module, and an optimization module; the molecular force field is constructed by encapsulating a pre-trained deep learning model.

[0007] The system is highly versatile, but it relies heavily on a pre-trained model. However, this model is based on sampling and iteratively selects the force field that best meets the requirements. This results in excessive computational load and a lack of targeted parameters, leading to poor simulation results in the output force field files.

[0008] Patent document CN110473596A discloses a force field fitting method for an aluminum electrolytic molten salt system based on potential energy surface scanning:

[0009] 1) Calculate the potential energy surface of any two ionic bonds in the molten salt system;

[0010] 2) Calculate the charge and electrostatic interaction potential energy of ions in molten salt; decompose the electrostatic interaction potential energy between ions to obtain the Buckingham potential energy;

[0011] 3) The Buckingham potential function was fitted using a nonlinear fitting method to obtain potential parameters suitable for the molecular dynamics of the molten salt system;

[0012] 4) Perform mechanical tests based on the obtained potential parameters, and calculate the ionic structure and transport properties of the molten salt system.

[0013] This method is only suitable for force field fitting of non-bonded interactions between ions and is not suitable for high-precision mechanical property analysis regarding strength. Summary of the Invention

[0014] To address the aforementioned issues, this invention provides a material force field fitting method based on potential energy surface matching. This method introduces the Morse bond stretching potential parameter with a finite potential well depth at infinity, making the final simulation results of yield strength, fracture strength, and fracture elongation more consistent with reality. This allows for a better understanding of the physical and mechanical properties and structure-property relationships of fiber materials.

[0015] A material force field fitting method based on potential energy surface matching includes:

[0016] S1 obtains the corresponding chemical structural formula based on the target material;

[0017] Based on the atomic hybridization state of the chemical structure in S1, S2 constructs a chemical structure model of the connecting bonds in the chemical structure.

[0018] S3 performs preliminary structural optimization on the chemical structure model obtained in S2 to obtain the energy-minimizing structure.

[0019] S4 is based on the energy-minimizing structure obtained in S3, and the chemical structure parameters of the target material are calculated by the DFT method.

[0020] Based on the chemical structure parameters obtained in S4, S5 uses the Morse algorithm to obtain the Morse bond stretching potential parameters and saves them in the molecular simulation file.

[0021] S6 repeats S2-S5 until all chemical structure models have been calculated and all molecular simulation files have been integrated into the force field file;

[0022] Based on the force field file obtained in S6, S7 performs molecular simulation to calculate the mechanical properties of the target material.

[0023] The method provided by this invention is for materials for which tensile yield strength, fracture strength and elongation at break are the main evaluation indicators.

[0024] Preferably, the target material is a fiber material.

[0025] Specifically, the chemical structure model in S2 includes a cluster structure model and a free radical fragment model after bond breaking.

[0026] Specifically, the initial structural optimization in S3 is based on the molecular force field of MMFFs, and the chemical structure model is optimized by combining the fastest descent method and the conjugate gradient method to facilitate subsequent quantitative calculations.

[0027] Specifically, the chemical structural parameters in S4 include equilibrium bond length, electronic energy, equilibrium vibrational frequency, and corresponding vibrational wavenumber.

[0028] Specifically, the process of S4 is as follows:

[0029] S4.1 The vibration frequency of the energy-minimizing structure is analyzed by secular equation to obtain the equilibrium vibration frequency;

[0030] S4.2 combined with S4.1 to obtain the mode shape diagram of the equilibrium vibration frequency, and the corresponding vibration wavenumber is obtained by infrared spectroscopy analysis;

[0031] S4.3 By setting the functional, basis set, static charge number, and spin multiplicity parameters, wavefunction analysis of the energy-minimizing structure is performed using density functional theory to obtain the corresponding electronic energy and equilibrium bond length.

[0032] Specifically, the Morse bond stretching potential parameters in S5 include the bond dissociation potential well depth, the equilibrium bond length, and the curvature of the potential energy curve controlling the equilibrium position.

[0033] Specifically, the process of S5 is as follows:

[0034] S5.1 calculates the bond dissociation potential well depth using the electronic energy difference between the cluster structure model and the post-bond-break radical fragment model:

[0035]

[0036] Where h is Planck's constant, c is the speed of light, and D ν Bond dissociation potential well expressed in wavelength units. Let i be the molar enthalpy of the i-th free radical fragment. For the molar enthalpy of the cluster model, ΔE c This is the temperature correction term, where R is the ideal gas constant, T is the temperature, and E is the temperature correction term. el,fr,i E represents the electron energy of the free fragment after bond breakage. el,cl The electronic energy of the cluster structure;

[0037] S5.2 Calculation of the curvature of the potential energy curve at the control equilibrium position:

[0038] α=(8π 2 cμν e χ / h) 1 / 2 =0.2454(Mν) e χ) 1 / 2 =0.1227ν e (M / D ν ) 1 / 2

[0039] Where M = M1M2 / (M1+M2), M is the reduced mass of the diatomic system, M1 and M2 are the molar masses of each atom in the diatomic system, χ is the anharmonicity constant, and v e The wavenumber is the frequency of vibration in equilibrium.

[0040] The mechanical properties of the target material include yield strength, fracture strength, and elongation at break.

[0041] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0042] (1) In the process of force field fitting, the Morse bond stretching potential parameters with more physical meaning are obtained by fitting the nonlinear Morse exponential function. The model is constructed by the force field file containing the Morse bond stretching potential parameters, and the model can better reflect the physical and mechanical properties and structure-property relationship of the material.

[0043] (2) Compared with the traditional method of directly using the sum of covalent bond radii as the bond breaking judgment, the present invention adopts a nonlinear Morse-type exponential function fitting method to avoid the bond-level reaction force field involving electron density integral calculation, thereby reducing the calculation time.

[0044] (3) The method of the present invention requires less experimental data in addition to the chemical structure of the basic target material, thereby reducing time costs and saving experimental testing costs. Attached Figure Description

[0045] Figure 1 A schematic flowchart illustrating the material force field fitting method provided by this invention;

[0046] Figure 2 This is a comparison of the traditional harmonic stretching potential curve and the Morse bond stretching potential curve of the CH3CN cluster model in the embodiment.

[0047] Figure 3 This is a comparison chart of simulation results between the traditional method and the method provided by the present invention in the embodiments;

[0048] Figure 4 This is a top view of a pre-oxidized fiber model under non-equilibrium stretching, generated by the method provided by the present invention in this embodiment.

[0049] Figure 5 This is a top view of the pre-oxidized fiber model generated by the conventional method in the embodiment, under non-equilibrium stretching. Detailed Implementation

[0050] like Figure 1 As shown, a material force field fitting method based on potential energy surface matching includes the following steps:

[0051] S1. Based on the target material, obtain the corresponding chemical structural formula: Taking the fiber pre-oxidized filament model system containing C / H / O / N elements as an example, the chemical structural formulas and corresponding bond types involved in this system are shown in Table 1:

[0052] Table 1 shows the chemical structural formulas and corresponding bond types included in the embodiments.

[0053] Cluster model key type <![CDATA[CH3CN]]> <![CDATA[C sp -C sp3 ]]> HCN <![CDATA[C sp -N sp ]]> HC(O)H <![CDATA[C sp2 -THE sp2 ]]> <![CDATA[CH4]]> <![CDATA[C sp3 -H]]> <![CDATA[C6H5NH2]]> <![CDATA[C ring -N sp3 ]]> <![CDATA[C6H5OH]]> <![CDATA[C ring -THE sp3 ]]> <![CDATA[C6H6]]> <![CDATA[C ring -C ring ]]> <![CDATA[C6H6]]> <![CDATA[C ring -H]]> <![CDATA[C5H5N]]> <![CDATA[C ring -N ring ]]> <![CDATA[NH3]]> <![CDATA[N sp3 -H]]> <![CDATA[CH3OH]]> <![CDATA[O sp3 -H]]> <![CDATA[C2H6]]> <![CDATA[C sp3 -C sp3 ]]>

[0054] Based on the chemical structural formula obtained in S1, and considering the atomic hybridization state of the chemical structural formula, S2 constructs a chemical structural model regarding the connecting bonds in the chemical structural formula:

[0055] With C sp -C sp3 Taking the bond as an example, atoms are supplemented based on the hybrid form, and a chemical structure model is established in Avogadro molecular modeling software according to the cluster model CH3CN. This model includes the cluster model CH3CN and the free radical fragment models ·-CH3 and ·-CN after bond breaking.

[0056] S3 performs preliminary structural optimization on the chemical structure model obtained in S2 to obtain the energy-minimizing structure: Based on the molecular force field of MMFFs, combined with the steepest descent method and the conjugate gradient method, the chemical structure model in S2 is optimized into the energy-minimizing structure, and the inp input file of the Orca quantization calculation software is generated in Extensions.

[0057] Based on the energy-minimizing structure obtained in S3, S4 calculates the chemical structure parameters of the target material using the DFT method: The keywords in the quantitative calculation input file are edited, selecting a functional such as B3LYP containing 20% ​​HF in the exchange-correlation functional; according to the final test accuracy requirements, a basis set type such as def2-SVP split valence bond basis set is selected; the keyword "opt" is used for structure optimization; the keyword "Freq" is used to perform vibrational frequency analysis using the secular equation; simultaneously, the net charge number and spin multiplicity parameters are set, and density theory analysis is performed on the cluster model CH3CN in the Orca software.

[0058] Finally, the structural electronic energy E of the optimal cluster model CH3CN was visualized using the software Chemcraft. el,c l = -8.3194 * 10 4 kcal / mol, equilibrium position bond length The vibrational wavenumber ν of the cluster model CH3CN was determined by combining the mode shape diagram and empirical values ​​from infrared spectroscopy analysis. e =1388.59cm -1 .

[0059] Using the same functional, basis set, and keyword "SP", density functional theory analysis was performed on the post-bonded radical fragment models ·-CH3 and ·-CN using Orca software. The final output files were visualized using Chemcraft software. The optimal structural electron energies of the post-bonded radical fragment models ·-CH3 and ·-CN were -2.4961 × 10⁻⁶. 4 kcal / mol and -5.8102*10 4 kcal / mol.

[0060] Based on the chemical structure parameters obtained in S4, S5 uses the Morse algorithm to obtain the Morse bond stretching potential parameters and saves them in the molecular simulation file:

[0061] Original Morse formula:

[0062]

[0063] Where r is the interatomic distance, r e Let α be the bond length at the equilibrium position, D be the depth of the bond dissociation potential well, and α be the curvature of the potential energy curve at the equilibrium position.

[0064] The discrete expression for vibrational energy is obtained based on multiple parameter transformations, generalized Laguerre polynomials, and fractional-order differential normalized integration algorithms:

[0065]

[0066] Where χ is the anharmonicity constant, νe The wavenumber corresponds to the equilibrium vibration frequency, and the anharmonicity constant is a function of the bond dissociation potential well depth and the vibration wavenumber or frequency.

[0067] Vibrational wavenumbers can be obtained through DFT vibrational analysis and spectral analysis, while the bond dissociation potential well depth can be calculated by solving the electronic energy difference between the cluster structure model and the post-bond breakage free radical fragment model using the DFT method.

[0068]

[0069] Where h is Planck's constant, c is the speed of light, and D ν Bond dissociation potential well expressed in wavelength units. Let i be the molar enthalpy of the i-th free radical fragment. For the molar enthalpy of the cluster model, ΔE c This is the temperature correction term, where R is the ideal gas constant, T is the temperature, and E is the temperature correction term. el,fr,i E represents the electron energy of the free radical fragment after bond breaking. el,cl The electronic energy is for the cluster structure; since the bonds in the diatomic system are calculated repeatedly, the final potential energy parameter D is taken as D0. e Half of the value was used to obtain a bond dissociation potential well depth of 65.80 kcal / mol.

[0070] Calculation of the curvature of the potential energy curve at the equilibrium position:

[0071] α=(8π 2 cμν e χ / h) 1 / 2 =0.2454(Mν) e χ) 1 / 2 =0.1227ν e (M / D ν ) 1 / 2

[0072] Where M = M1M2 / (M1+M2) is the reduced mass of the diatomic system, χ is the anharmonic constant, and v e The wavenumber corresponding to the equilibrium vibration frequency is given; the curvature of the potential energy curve at the equilibrium position is obtained.

[0073] Substitute the parameters obtained above into the fitted Morse-type bond stretching potential expression:

[0074] E Morse =65.80*[1-e -1.945*(r-1.459) ] 2

[0075] A curve was plotted based on the fitted Morse-type bond stretching potential, and compared with a curve obtained by a traditional method:

[0076] The results are as follows Figure 2 As shown, the fitting results of this method have good consistency with the traditional harmonic stretching potential curve at the equilibrium position (i.e., the interatomic distance in the figure is around 1.5). At the same time, bond breakage occurs when the atomic distance is far (i.e., the interatomic distance in the figure is around 2.5). Compared with the traditional harmonic stretching potential curve, it is more in line with the actual situation. The dashed line in the figure is the traditional harmonic stretching potential curve, and the solid line is the Morse bond stretching potential curve.

[0077] S6 repeats S2-S5 until all chemical structural formulas are calculated, and then integrates all molecular simulation files into the force field file:

[0078] Write the parameters into the .in and .data files of LAMMPS.

[0079] The format of the in file is: bond_style morse

[0080] The data file format is: Bond Coeffs${Bond#}65.80 1.945 1.459.

[0081] The bond types and fitting parameters of the fiber pre-oxidized fiber model in this embodiment are shown in Table 2:

[0082]

[0083] Based on the force field file obtained in S6, S7 performs molecular simulation to calculate the mechanical properties of the target material, including yield strength, fracture strength, and elongation at break.

[0084] The simulation results are compared with traditional methods, such as Figure 3 As shown in the figure, the dashed line represents the stretching potential stress field curve of the traditional harmonic bond, and the solid line represents the stretching potential stress field curve of the Morse bond. According to the comparison, the stretching potential stress of the traditional harmonic bond increases linearly with the increase of strain, and no sudden drop was observed at 9% nominal strain (i.e., no bond breakage point appeared); while the stretching potential stress-strain curve of the Morse bond shows a sudden drop at 9% nominal strain (i.e., bond breakage point appears), and its process is more consistent with the actual situation.

[0085] like Figure 4 , 5 As shown in the figure, the thick black bars represent nanopores or microcracks that appear between atoms. In the simulation process, the model generated by the method of this invention can clearly observe the broken bond position when the nominal strain reaches 9%, which is consistent with the actual situation; while the traditional method does not show obvious nanopores or microcracks, which is inconsistent with the actual situation.

[0086] Finally, according to the simulation results, the breaking strength of the pre-oxidized fiber in this embodiment is 7.95 GPa, and the breaking elongation is 8.83%.

Claims

1. A material force field fitting method based on potential energy surface matching, characterized in that, include: S1 obtains the corresponding chemical structural formula based on the target material; Based on the atomic hybridization state of the chemical structure in S1, S2 constructs a chemical structure model of the connecting bonds in the chemical structure. S3 performs preliminary structural optimization on the chemical structure model obtained in S2 to obtain the energy-minimizing structure. S4 is based on the energy-minimizing structure obtained in S3, and the chemical structure parameters of the target material are calculated by the DFT method. S5, based on the chemical structure parameters obtained in S4, uses the Morse algorithm to obtain the Morse bond stretching potential parameters and saves them in the molecular simulation file. The specific process of S5 is as follows: S5.1 calculates the bond dissociation potential well depth using the electronic energy difference between the cluster structure model and the post-bond-break radical fragment model: Where h is Planck's constant, c is the speed of light, and D ν Bond dissociation potential well expressed in wavelength units. Let i be the molar enthalpy of the i-th free radical fragment. For the molar enthalpy of the cluster model, ΔE c This is the temperature correction term, where R is the ideal gas constant, T is the temperature, and E is the temperature correction term. el,fr,i E represents the electron energy of the free radical fragment after bond breaking. el,cl The electronic energy of the cluster structure; S5.2 Calculation of the curvature of the potential energy curve at the control equilibrium position: α=(8π 2 cmn e h / h) 1 / 2 =0.2454(Mν e x) 1 / 2 =0.1227ν e (M / D ν ) 1 / 2 Where M = M1M2 / (M1+M2), M is the reduced mass of the diatomic system, M1 and M2 are the molar masses of each atom in the diatomic system, χ is the anharmonicity constant, and ν e The wavenumber corresponding to the equilibrium vibration frequency; S6 repeats S2-S5 until all chemical structure models have been calculated and all molecular simulation files have been integrated into the force field file; Based on the force field file obtained in S6, S7 performs molecular simulation to calculate the mechanical properties of the target material.

2. The material force field fitting method according to claim 1, characterized in that, The target material in S1 is a fiber material.

3. The material force field fitting method according to claim 1, characterized in that, The chemical structure model in S2 includes a cluster structure model and a free radical fragment model after bond breaking.

4. The material force field fitting method according to claim 1, characterized in that, The initial structural optimization in S3 is based on the molecular force field of MMFFs, and the chemical structure model is optimized by combining the fastest descent method and the conjugate gradient method.

5. The material force field fitting method according to claim 1, characterized in that, The chemical structural parameters in S4 include equilibrium bond length, electronic energy, equilibrium vibrational frequency, and corresponding vibrational wavenumber.

6. The material force field fitting method according to claim 1 or 5, characterized in that, The specific process of S4 is as follows: S4.1 The vibration frequency of the energy-minimizing structure is analyzed by secular equation to obtain the equilibrium vibration frequency; S4.2 combined with S4.1 to obtain the mode shape diagram of the equilibrium vibration frequency, and the corresponding vibration wavenumber is obtained by infrared spectroscopy analysis; S4.3 uses density functional theory to perform wave function analysis on the energy-minimizing structure to obtain the corresponding electronic energy and equilibrium bond length.

7. The material force field fitting method according to claim 1, characterized in that, The Morse bond stretching potential parameters in S5 include the bond dissociation potential well depth, equilibrium bond length, and curvature of the potential energy curve controlling the equilibrium position.

8. The material force field fitting method according to claim 1, characterized in that, The mechanical properties of S7 include yield strength, fracture strength, and elongation at break.

Images