Melting point calculation method based on first-principle molecular dynamics and bayesian statistics

Abstract

This invention belongs to the field of computational materials and proposes a melting point calculation method based on first-principles molecular dynamics and Bayesian statistics. It employs a small-scale system for calculation, significantly reducing computational complexity compared to large-scale systems. Furthermore, it utilizes an on-the-fly machine learning force field to further accelerate melting point calculation. This method requires no manual intervention, automatically updating the force field during AIMD simulations to maintain the accuracy of first-principles calculations and significantly improve computational efficiency. This method calculated 16 solid-liquid configurations within the 800–1100 K range, consuming approximately 600 cores, which significantly shortens computation time compared to other methods, thus demonstrating great application potential.

Description

Technical Field

[0001] This invention belongs to the field of computational materials, specifically relating to a melting point calculation method based on first-principles molecular dynamics and Bayesian statistics. Background Technology

[0002] Theoretical predictions of material melting points have been around for a long time. In the past two decades, with the improvement of computing power, first-principles molecular dynamics (AIMD) simulations based on density functional theory (DFT) have become accurate and common simulation methods. However, DFT-based melting point predictions are still quite challenging due to limitations in computational accuracy and time. Currently, commonly used methods for melting point calculation include, for example, the free energy-based method [ellers Michael S, Lísal Martin, Brennan John K. Free-energy calculations using classical molecular simulation: application to the determination of the melting point and chemical potential of a flexible RDX model]. The melting point is located at the intersection of the solid and liquid free energy curves. This method requires very high accuracy in Gibbs free energy calculation because the angle between the solid and liquid free energy curves is very small; even a small error in the free energy calculation can cause a large error in the melting temperature. Another mainstream method is the solid-liquid coexistence method [Alfè Dario. Temperature of the inner-core boundary of the Earth: Melting of iron at high pressure from first-principles coexistence simulations]. The melting point is determined by simulating the stable temperature of the solid-liquid coexistence system using AIMD's NVE ensemble. Because there is a solid-liquid interface between the solid and liquid phases, the system requires a large unit cell (typically at least 1000 atoms) to remain stable, and requires a long simulation time (typically at least 20 ps) to obtain results. Summary of the Invention

[0003] To address the shortcomings of existing simulation methods, we propose a melting point calculation method based on first-principles molecular dynamics and Bayesian statistics to solve the problems of excessively long calculation times and high costs. We reduce the size of the simulation system, establishing a small-scale solid-liquid coexistence system (approximately 200 atoms). During the AIMD calculation, we use machine learning force fields combined with Bayesian methods to process the results, thereby reducing computational costs.

[0004] The technical solution of this invention:

[0005] A melting point calculation method based on first-principles molecular dynamics and Bayesian statistics, comprising the following steps:

[0006] S1. Optimize the initial unit cell structure and establish a supercell;

[0007] Unit cell structures were downloaded from the open-source materials database Materials Project or manually constructed according to Perason's Handbook. First-order calculations were performed on the unit cell structures using the DFT-based (Vienna Ab-initio Simulation Package) VASP software. After obtaining the optimized unit cell structures, supercell structures were fabricated, with supercell side lengths greater than [missing value].

[0008] S2. Calculate and estimate the thermal expansion at the melting point temperature;

[0009] The NPT ensemble was supercelled to the experimental temperature or the average temperature of each component and relaxed for 10 ps. The average lattice constant of the final 500–1000 configurations was taken as the lattice constant after thermal expansion. Then, the NVT ensemble was relaxed for 5 ps at the same temperature to bring the system to equilibrium. In the first-principles molecular dynamics (AIMD) calculations, if the size exceeds 100 atoms, the Gamma parameter was used to select the k-grid point, and the Brillouin zone was selected with a 1×1×1 grid point.

[0010] S3. Expand the thermally expanded supercell along the Z direction, fix half of the atoms, heat the other half of the supercell for 10 ps at 1000 K above the estimated melting point temperature to completely melt it, and then continue to relax it at the same temperature for 4 ps. Take a solid-liquid coexistence configuration for every 1 ps, and take a total of 4 solid-liquid coexistence configurations to prepare for the next relaxation step.

[0011] S4. Place each of the four solid-liquid coexistence configurations within the estimated melting point temperature range (T). m -200 K)~

[0012] (T m Relaxation was performed at +200K for 10–20 ps, ​​followed by analysis of the solid-liquid distribution within this range. The Steinhardt bond order parameters were used to determine the number of solid atoms in the relaxed solid-liquid mixture. The method is as follows: Given a central atom, its nearest neighbor bonds are projected onto a unit sphere. Based on these projection vectors, a set of local bond order parameters is defined. The bond order parameters of the structure surrounding particle i are defined as follows:

[0013]

[0014] Among them, Y lm It is a spherical harmonic, where N(i) is the number of nearest-neighbor atoms of particle i. It is the vector connecting particles i and j, where l and m are both integers, and m takes values ​​between -l and +l; q is defined in formula (1). lm (i) The value is affected by the choice of coordinate system. A key order parameter that is not affected by the choice of coordinate system is defined as:

[0015]

[0016] This bond order parameter is a quantification of the local order around particle i; the q6 parameter is used to identify the crystal structure and distinguish between solid and liquid atoms; for two adjacent atoms i and j, if the parameter:

[0017]

[0018] It is assumed that the two atoms are "connected"; in the formula: s ij is the scalar product of the correlation between the two atomic structures i and j; threshold is the threshold value, which is set to 0.5;

[0019] Bokeloh proposed using a second-order bond order parameter to distinguish between solid and liquid atoms; this parameter is derived from s ij The average value is defined as follows:

[0020] ij >>average threshold (4)

[0021] In the formula: ij >for s ij The average value; if an atom's Steinhardt bond order parameter satisfies both formulas (3) and (4), then the atom is determined to be a solid atom, otherwise it is a liquid atom; the average threshold value needs to be corrected according to the solid-liquid coexistence configuration. For the solid-liquid coexistence configuration, since half of the solid atoms are fixed and the other half of the solid atoms are heated and melted to make them liquid, the number of solid atoms in the solid-liquid coexistence configuration is 50%. The average threshold parameter needs to be corrected during calculation to ensure the proportion of solid atoms; if the number of liquid atoms in the system accounts for 80% of the total number of atoms, the configuration is considered to be a liquid state, otherwise it is a solid state; finally, the Bayesian method is used to process the statistical data of the solid-liquid state. This method not only provides point estimates, but also makes statistical inferences on the uncertainty of the estimates; thus, the probability curves of the liquid configuration at different temperatures are obtained, and the temperature corresponding to a probability of 50% is the melting point.

[0022] Compared with existing technologies, the beneficial effects of this invention are as follows:​​

[0023] This invention employs a method based on first-principles molecular dynamics (AIMD) to calculate the melting point of a small-scale solid-liquid coexisting system (approximately 200 atoms). By incorporating machine learning force fields during the calculation, the computational workload can be significantly reduced. Bond sequence parameters are used to determine the solid-liquid state, and Bayesian methods are then used to process the solid-liquid distribution data to obtain a curve showing the melting probability as a function of temperature. Attached Figure Description

[0024] Figure 1 This is a schematic diagram of an algorithm for simulating the melting point of a coexisting phase;

[0025] Taking the melting point calculation of pure aluminum as an example, the initial unit cell structure of aluminum is first optimized by expanding the cell in a 3×3×3 manner; the melting point temperature (900K) is calculated and estimated; then the cell is expanded along the Z direction, and half of the supercell is melted to prepare a solid-liquid coexistence configuration. Green represents solid atoms, and brown represents liquid atoms.

[0026] Figure 2 The Steinhardt bond order parameter is used to determine the solid-liquid atom distribution map;

[0027] First, the number of solid atoms in the solid-liquid coexistence configuration is calculated using the Steinhardt bond order parameter. The avgthreshold parameter is initialized so that the number of solid atoms is 108 and the total number of atoms is 216. Then, this parameter can be used to determine the solid-liquid state of the configuration at different temperatures. For example, the experimental melting point of aluminum is 933 K. When the coexistence configuration relaxes for 10 ps at 800 K, the number of solid atoms is 216, which exceeds 80% of the total number of atoms, so it can be determined to be a solid configuration. When it relaxes for 10 ps at 1100 K, the number of solid atoms is 0, so it can be determined to be a liquid configuration.

[0028] Figure 3 The melting point fitting is as follows: (a) shows the distribution of the number of 200 liquid atoms in each configuration within a range of ±100K at different temperatures. When the number of liquid atoms is greater than 80%, the configuration is considered to be liquid after further calculation. (b) shows the Bayesian fitted temperature curve, calculated using Bayesian statistical methods on the processed data. The red curve in the figure corresponds to the temperature at which there is a 50% probability, which is the melting point. Detailed Implementation

[0029] The specific embodiments of the present invention will be further described below with reference to the accompanying drawings and technical solutions.

[0030] The calculation process is as follows, taking the melting point of pure aluminum as an example.

[0031] S1. Optimize the initial unit cell structure and establish a supercell.

[0032] Aluminum unit cell structures were downloaded from the open-source materials database Materals Project. First-order calculations were performed using the Vienna Ab-initio Simulation Package (VASP) software, based on the Discrete Functional Theory (DFT). Electron exchange correlation functionals were calculated using the Localized Density Approximation (LDA). Structural optimization was performed using a plane wave cutoff energy of 400 eV and an energy change convergence threshold of 10 eV. -7 eV / atom; the Monkhorst-Pace method was used to select k grid points. For bulk structures, 15×15×15 grid points were selected for the Brillouin zone; the pseudopotential was selected as the projected augmented plane wave (PAW)

[90] potential, with valence electrons of Al: 3s2 3p1. After structural optimization, the cell was expanded in a 3×3×3 manner to prepare a supercell structure, such as Figure 1 As shown.

[0033] S2. Calculate and estimate the thermal expansion at the melting point temperature.

[0034] After cell expansion, the aluminum supercell was relaxed to the estimated temperature for 10 ps using the NPT ensemble. The average lattice constant of the final 500 configurations was calculated as the lattice constant after thermal expansion. Then, the system was relaxed for another 5 ps using the NVT ensemble to reach equilibrium. The Gamma method was used to select the k-grid points, and 1×1×1 grid points were selected for the Brillouin zone. Machine learning was used to accelerate the computation.

[0035] S3. Expand the cell along the Z direction and melt half of the supercell to prepare multiple solid-liquid coexisting configurations.

[0036] The thermally expanded supercell is expanded along the Z direction, half of the aluminum atoms are fixed, and the other half of the supercell is heated at 2000K for 10 ps to completely melt it. Then, it is relaxed for another 4 ps. A solid-liquid coexistence configuration is taken for every 1 ps, and a total of 4 configurations are taken to prepare for the next step.

[0037] S4. Collect the solid-liquid distribution results and use Bayesian statistical methods to fit the melting point (T). m )

[0038] The solid-liquid coexistence configurations were relaxed for 20 ps within the estimated melting point temperature range (800–1100 K), and the solid-liquid distribution results were collected. The Steinhardt bond order parameters for each configuration were calculated using the Python-based PyScal package to determine the number of solid atoms in each configuration. First, calculations were performed... Figure 2 In the aluminum solid-liquid coexistence model, the number of solid atoms was found to be [number missing] when the average threshold was 0.6. Figure 2 The solid-liquid coexistence model has 10⁸ solid atoms, which is 50%, and the average threshold verification is complete. Then, the number of solid atoms in the relaxed configuration at each temperature is calculated, such as... Figure 2The solid atom count of Al₂O₃ in its 800K configuration is 216; the solid atom count in its 1100K configuration is 0. When the solid atom count is greater than 80%, the configuration is considered to be in a solid state after further calculation; otherwise, it is considered to be in a liquid state. The percentage of liquid atoms for each configuration is as follows: Figure 3 As shown in (a). After determining the solid-liquid state of the configuration, all atoms in the solid configuration are considered solid atoms, and vice versa, corresponding to... Figure 3 (b) shows the black dots with melting probabilities of 0 and 100. Bayesian statistical methods are then used to calculate the processed data, fitting the number of atoms in the processed solid and liquid to obtain the curve of melting probability versus temperature. Figure 3 (b) The curve corresponds to the temperature of 941±3K with a 50% probability, which is the melting point of aluminum, and is basically consistent with the experimental melting point of 933K. This method calculated 16 small solid-liquid configurations in the range of 800–1100K, consuming about 600 cores. Compared with other methods (usually more than 5000 cores), it significantly shortens the calculation time and does not require complex free energy calculations, thus showing great application potential.

Claims

1. A melting point calculation method based on first-principles molecular dynamics and Bayesian statistics, characterized by, The steps are as follows: S1, optimize the initial unit cell structure, and establish the supercell; The unit cell structure was downloaded from the open source material database Materals Project or manually built according to Perason's Handbook; the first principle calculation was performed on the unit cell structure by using the VASP software based on DFT, and the supercell structure was prepared after obtaining the optimized unit cell structure, and the side length of the supercell was greater than S2, calculate the thermal expansion at the estimated melting point temperature; Use the NPT ensemble supercell to relax for 10 ps at the experimental temperature or the average temperature of each component, and take the average lattice constant of the final 500-1000 configurations as the lattice constant after thermal expansion; then use the NVT ensemble to relax for 5 ps at the same temperature, so that the system reaches equilibrium; during the first-principles molecular dynamics calculation, if the size exceeds 100 atoms, use the Gamma parameter to select the k grid point, and select a 1x1x1 grid point in the Brillouin zone; S3, expand the supercell along the Z direction after thermal expansion, fix half of the atoms, heat the other half of the supercell above the estimated melting point temperature by 1000 K for 10 ps, so that it is completely melted, then continue to relax for 4 ps at the same temperature, take one solid-liquid coexistence configuration every 1 ps, a total of 4 solid-liquid coexistence configurations for the next step of relaxation; S4, put the four solid-liquid coexistence configurations into the estimated melting point temperature range (T m -200 K)~(T m +200 K) to relax 10~20 ps, and then analyze the solid-liquid distribution results in the range of 10~20 ps; adopt Steinhardt bond order parameter to judge the number of solid atoms in the relaxation, and the method is as follows: give a central atom, project its near neighbor bond to a unit sphere; based on these projection vectors, define a set of local bond order parameters, and the bond order parameter of the structure around particle i is defined as: where Y lm is the spherical harmonic, N(i) is the number of nearest neighbors of particle i, is the vector connecting particle i and j, and l and m are integers, m taking values between -l and +l; q lm (i) is dependent on the choice of the coordinate system, and a bond order parameter that is independent of the choice of the coordinate system is defined as This bond order parameter is a quantification of the local order around particle i; the q6 parameter is used to identify the crystal structure and distinguish solid and liquid atoms; for two adjacent atoms i and j, if the parameter: two atoms are considered to be "connected"; where: s ij is the scalar product of the correlation between the two atomic structures i and j; threshold is the threshold value, taken to be 0.5; Bokeloh proposed to distinguish solid from liquid atoms by a second order bond order parameter, which is defined by s ij The average value of s is defined as follows: ij >>average threshold (4) Where: ij >for s ij The average value; if an atom's Steinhardt bond order parameter satisfies both formulas (3) and (4), then the atom is determined to be a solid atom, otherwise it is a liquid atom; the average threshold value needs to be corrected according to the solid-liquid coexistence configuration. For the solid-liquid coexistence configuration, since half of the solid atoms are fixed and the other half of the solid atoms are heated and melted to make them liquid, the number of solid atoms in the solid-liquid coexistence configuration is 50%. The average threshold parameter needs to be corrected during calculation to ensure the proportion of solid atoms; if the number of liquid atoms in the system accounts for 80% of the total number of atoms, the configuration is considered to be a liquid state, otherwise it is a solid state; finally, the Bayesian method is used to process the statistical data of the solid-liquid state. This method not only provides point estimates, but also makes statistical inferences on the uncertainty of the estimates; thus, the probability curves of the liquid configuration at different temperatures are obtained, and the temperature corresponding to a probability of 50% is the melting point.​​

Images