Machine learning molecular dynamics simulations toward exploration of high-temperature properties of nuclear fuel materials: case study of thorium dioxide

Predicting materials properties of nuclear fuel compounds is a challenging task in materials science. Their thermodynamical behaviors around and above the operational temperature are essential for the design of nuclear reactors. However, they are not easy to measure, because the target temperature range is too high to perform various standard experiments safely and accurately. Moreover, theoretical methods such as first-principles calculations also suffer from the computational limitations in calculating thermodynamical properties due to their high calculation-costs and complicated electronic structures stemming from f-orbital occupations of valence electrons in actinide elements. Here, we demonstrate, for the first time, machine-learning molecular-dynamics to theoretically explore high-temperature thermodynamical properties of a nuclear fuel material, thorium dioxide. The target compound satisfies first-principles calculation accuracy because f-electron occupation coincidentally diminishes and the scheme meets sampling sufficiency because it works at the computational cost of classical molecular-dynamics levels. We prepare a set of training data using first-principles molecular dynamics with small number of atoms, which cannot directly evaluate thermodynamical properties but captures essential atomistic dynamics at the high temperature range. Then, we construct a machine-learning molecular-dynamics potential and carry out large-scale molecular-dynamics calculations. Consequently, we successfully access two kinds of thermodynamic phase transitions, namely the melting and the anomalous \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}λ transition induced by large diffusions of oxygen atoms. Furthermore, we quantitatively reproduce various experimental data in the best agreement manner by selecting a density functional scheme known as SCAN. Our results suggest that the present scale-up simulation-scheme using machine-learning techniques opens up a new pathway on theoretical studies of not only nuclear fuel compounds, but also a variety of similar materials that contain both heavy and light elements, like thorium dioxide.


The parameters of symmetry functions and parity plot
In this section, we show the parameters of symmetry functions employed in the main text and the parity plots of the energies and forces given by BPNN as compared to DFT reference values.
We generated symmetry functions with the combination of parameters: η  Table S1, S2, and S3 Type of symmetry function Combination η (2), (4)   Th-Th-Th 0.02367 -1 1.0 6.5 type-4 Th-Th-Th 0.02367 1 1.0 6.5 type-4 Th-Th-Th 0.02367 1 4.0 6.5 type-4 Th-Th-Th 0.05945 1 1.0 6.5  Th-Th-Th 0.02367 1 1.0 6.5 type-4 Th-Th-Th 0.02367 1 4.0 6.5    Using the selected descriptor vector listed in Table S1, S2, and S3, we trained BPNNs using the multistream Kalman filter method 3,4 with 30 epochs. 90% of the reference data was assigned to training data and the remaining 10% to test data. Both root mean square errors (RMSE) decreased rapidly in a few epochs as shown in figure S1. The RMSE of energies continued to show small oscillations even as the epoch progressed, but the RMSE of forces converged to an almost constant value. We selected the BPNN with the smallest RMSE of energies for the test data within 30 epochs. The resulting RMSEs were summarized in table 1 in the main text. The parity plots of the energies and forces obtained by BPNNs are shown in figure S2. We can find that a good linear relationship is realized for both energies and forces obtained by DFTs and BPNNs. Comparing BPNN-LDA, BPNN-PBEsol, and BPNN-SCAN, the ranges of forces included in the reference data of BPNN-LDA and BPNN-SCAN were broader than that of BPNN-PBEsol. The difference is considered due to the method used to generate the reference data set. The reference data set of BPNN-PBEsol was generated by only FPMD simulation, whereas these of BPNN-LDA and BPNN-SCAN were mainly created by adiabatic learning scheme. The structures including large forces acting on atoms do not frequently appear in FPMD simulation. On the other hand, the active learning scheme explores structures that are not sufficiently contained in the reference data. Therefore, the range of forces in the reference data of BPNN-LDA and BPNN-SCAN became broader than that of BPNN-PBEsol.

Validation of BPNNs for dynamical calculation
The validation of BPNNs in the main text focused on static calculations. Therefore, in this section, we show the validation of BPNNs for dynamical calculations. First, we show the stability of MLMD for dynamical calculation. MD simulation with BPNN trained with inappropriate descriptors, hyperparameters, and reference data sometimes be-comes unstable and show structural collapse with a long simulation period, especially at high temperatures 5 . Therefore, we checked the stability of MLMD with a long simulation period at high temperatures. We conducted MLMD NVE simulations at 2000, 3000, and 4000 K. The number of atoms was 2592, and the total simulation time was 1000 ps with 1 fs step size. To verify the energy conservation of MLMD NVE simulations, we show the time evolution of the error for total energy in figure S3, which is defined as where E(t) is the total energy at time t and N is the number of atoms. The errors obtained by MLMDs with BPNN-LDA, BPNN-PBEsol, and BPNN-SCAN were within ±5 × 10 −6 eV per atom, and the total energies computed by MLMDs were well conserved over a long period with no anomalies. These results denote that BPNNs can be smoothly fitted to the potential energy surface of DFTs. In this subsection, we compare the mean square displacement (MSD) for short trajectories obtained by FPMDs and MLMDs. We conducted FPMD and MLMD NVT simulations with 2 × 2 × 2 supercell (96 atoms). The lattice constant was fixed to 5.8Å. Time step size and total simulation time of NVT simulation were 2 fs and 25 ps, respectively. A Nosè-Hoover thermostat were applied for the simulations with relaxation times of 0.08 ps. We calculated ensemble average of MSD from t = 0 to 5 ps time period as

Comparison of the mean square displacements obtained by FPMD and MLMD
where N is the number of oxygen atoms and T 0 = 20 ps. Figure S4 shows the MSD obtained by FPMD and MLMD at 2000 and 3500 K. We can find that the results obtained by FPMD and MLMD are in good agreement.  This subsection shows the validation of thermophysical properties obtained by BPNN. Since all reference data about BPNN-PBEsol was generated by FPMD at various temperatures, thermophysical properties obtained by FPMD and MLMD based on PBEsol functional can be compared. Here, we compared the temperature dependence of the enthalpy H(T ) and the lattice constant L(T ) from 300 to 3500 K using FPMD and MLMD NPT. For both methods, we use 3 × 3 × 3 supercells (324 atoms). The time step size and simulation total time of FPMD at each temperature were 2 fs and 16 ps, respectively. MLMD was conducted with 2 fs time step and a total 100 ps simulation time. The results shown in figure S5 reveal that the temperature dependence of enthalpy H(T ) and lattice constant L(T ) obtained by MLMD shows a good agreement with that computed by FPMD.  In this section, we show the detail of smoothing the curve of heat capacity. We performed MLMD NPT simulations with 2592 atoms and total 200 ps run per 10 K temperature step. The molar specific heat capacity was obtained from the numerical differentiation of the enthalpy

Computational efficiency of BPNNs
where n is the amount of substance in moles, H is the enthalpy and dT is temperature step (dT = 10). In order to smooth the curve of specific heat capacity calculated by Eq (3), we averaged specific heat capacity over the interval of ±100 K twice 6 as where M = 10. Figure S6 shows the curves of specific heat with no averaging operation, averaging operation once and averaging operation twice. In figure S7, we also show the curves of the coefficient of linear thermal expansion (CLTE) defined as where L(T ) is the lattice constant at temperature T . The smoothing of the curve was done in the same way as done in the specific heat capacity curve. The ACLTE defined in the main text is equivalent to the average of CLTE from 300 to 1600 K.  Figure S8.(a) shows the size dependence of specific heat capacity by changing the number of atoms in the range from 96 (2 × 2 × 2 cell) to 2592 (6 × 6 × 6 cell). The temperatures of λ-peak became lower with increasing system size, and the shapes of λ-peak were almost unchanged above 1500 atoms. Therefore, the system size with 2592 atoms is considered sufficient to neglect finite size effects for calculating thermal properties using MLMD. The size dependence of CLTE was also small above 768 atoms (4 × 4 × 4 cell) as shown in figure  S8(b). The size dependence for the thermal properties of ThO 2 using MLMD was similar to the previous study about that of UO 2 using empirical force field 6 .
where O i mean the position of an interstitial oxygen atom on 4b site. Finally, the thorium Frenkel pair defect energies (ThFP) were calculated for the second nearest neighbor configuration with a thorium interstitial (4b site) and vacancy (4a site). The defect formation energies obtained by DFTs and BPNNs are listed in table S5. We can find that BPNNs reproduced the defect formation energies with a certain degree of accuracy. Here, we again note that we did not explicitly include the above defect structures in the present reference data. It is considered that the inclusion of liquid structures and the structures with oxygen diffusion at high temperatures made BPNNs possible to describe the various defect structures. However, comparing the errors for the defect formation energies obtained by BPNNs, there were some differences in accuracy against DFT data. Especially, BPNN-PBEsol showed a large error for the ThFP formation energy. The difference in accuracy for DFT data is considered due to differences in the method used to generate the reference data. The reference data set of BPNN-PBEsol was generated by only FPMD simulation, whereas these of BPNN-LDA and BPNN-SCAN were mainly created by the adiabatic learning scheme. The ranges of forces contained in the reference data of BPNN-LDA and BPNN-SCAN were broader than that of BPNN-PBEsol, as shown in figure S2. The reference data of BPNN-LDA and BPNN-SCAN contained more diverse structures than that of BPNN-PBEsol, which enable BPNN-LDA and BPNN-SCAN to predict the ThFP formation energy with a certain degree of accuracy. Of course, the explicit inclusion of the defect structures in the reference data is expected to improve further the accuracy of BPNNs for the defect formation energies.

Born effective charges and dielectric constants
In phonon calculation in the main text, we used the Born effective charge and the dielectric constants to calculate the non-analytical term correction to the dynamical matrix. They were obtained from the response to finite electric fields |E| = 0.01 (eV/Å) using VASP. The phonon dispersion curves computed by (BPNN-)LDA, (BPNN-)PBEsol, and (BPNN-)SCAN are corrected using the dipole-dipole interaction correction 8,9 with the Born effective charge and the dielectric constants obtained by LDA, PBEsol, and SCAN, respectively. The detailed values were summarized in Table S6. In phonon calculations using BD08 and cooper potentials, we used the non-analytical term correction with the Born effective charge and the dielectric constants obtained by LDA calculation.